Index: trunk/tools/nc_var_tools.py
===================================================================
--- trunk/tools/nc_var_tools.py	(revision 457)
+++ trunk/tools/nc_var_tools.py	(revision 458)
@@ -14864,4 +14864,936 @@
     return tB
 
+def wdismean(pos,val4):
+    """ Function to compute the mean value weighted to its 4 distances
+      pos= x,y position within the 'square'
+      val4= values at the four vertexs (2,2)
+    >>>position = [0.005,0.005]
+    >>>val4 = np.zeros((2,2), dtype=np.float)
+    >>>val4 = np.arange(4).reshape(2,2)
+    0.035707955234
+    """
+    fname = 'wdismean'
+
+    dist = np.zeros((2,2), dtype=np.float)    
+    dist[0,0] = np.sqrt(pos[0]*pos[0]+pos[1]*pos[1])
+    dist[0,1] = np.sqrt(pos[0]*pos[0]+(1.-pos[1])*(1.-pos[1]))
+    dist[1,0] = np.sqrt((1.-pos[0])*(1.-pos[0])+pos[1]*pos[1])
+    dist[1,1] = np.sqrt((1.-pos[0])*(1.-pos[0])+(1.-pos[1])*(1.-pos[1]))
+
+    if np.min(dist) == 0.:
+        pos0 = index_mat(dist, 0.)
+        dist[:,:] = 0.
+        posmean = val4[pos0[0], pos0[1]]
+    else:
+        totdist = np.sum(1./dist)
+        posmean = np.sum(val4.flatten()/dist.flatten())/totdist
+
+    return posmean
+
+def radial_points(angle,dist):
+    """ Function to provide a number of grid point positions for a given angle
+      angle= desired angle (rad)
+      dist= number of grid points
+    >>> radial_points(0.5*np.pi,5)
+    [[  6.12323400e-17   1.00000000e+00]
+     [  1.22464680e-16   2.00000000e+00]
+     [  1.83697020e-16   3.00000000e+00]
+     [  2.44929360e-16   4.00000000e+00]
+     [  3.06161700e-16   5.00000000e+00]]
+    """
+    fname = 'radial_points'
+    radpos = np.zeros((dist,2), dtype=np.float)
+
+    for ip in range(dist):
+        radpos[ip,0] = (ip+1.)*np.cos(angle)
+        radpos[ip,1] = (ip+1.)*np.sin(angle)
+
+    return radpos
+
+def compute_tevolboxtraj_radialsec(values, ncfile, varn):
+    """ Function to compute tevolboxtraj_radialsec: temporal evolution at a given point along 
+          a number of radii at a given frequency following a trajectory
+      values= [trajfile]@[Tbeg],[lonname],[latname],[zname],[timename],[timekind],[Nangle],[radii]
+        [trajfile]: ASCII file with the trajectory ('#' not readed)
+          [time] [xpoint] [ypoint]
+        [Tbeg]: equivalent first time-step of the trajectory within the netCDF file
+        [lonname],[latname],[zname],[timename]: longitude, latitude, z and time variables names
+        [timekind]: kind of time
+          'cf': cf-compilant
+          'wrf': WRF kind
+        [Nangle]: Number of angles
+        [radii]: length in grtid points of the angle
+      ncfile= netCDF file to use
+      varn= ',' list of variables' name ('all', for all variables)
+    """
+    import numpy.ma as ma
+    import subprocess as sub
+
+    fname='compute_tevolboxtraj_radialsec'
+
+    if values == 'h':
+        print fname + '_____________________________________________________________'
+        print compute_tevolboxtraj_radialsec.__doc__
+        quit()
+
+    arguments = '[trajfile]@[Tbeg],[lonname],[latname],[zname],[timename],' +        \
+      '[timekind],[Nangle],[radii]]'
+    check_arguments(fname,len(arguments.split(',')),arguments,',',8)
+
+    trajfile = values.split(',')[0].split('@')[0]
+    Tbeg = int(values.split(',')[0].split('@')[1])
+    lonn = values.split(',')[1]
+    latn = values.split(',')[2]
+    zn = values.split(',')[3]
+    timn = values.split(',')[4]
+    timekind = values.split(',')[5]
+    Nangle = int(values.split(',')[6])
+    radii = int(values.split(',')[7])
+
+    if not os.path.isfile(ncfile):
+        print errormsg
+        print '  ' + fname + " netCDF file: '" + ncfile + "' does not exist !!"
+        quit(-1)
+
+    if not os.path.isfile(trajfile):
+        print errormsg
+        print '  ' + fname + " trajectory file: '" + trajfile + "' does not exist !!"
+        quit(-1)
+
+    objfile = NetCDFFile(ncfile, 'r')
+    lonobj = objfile.variables[lonn]
+    latobj = objfile.variables[latn]
+    timobj = objfile.variables[timn]
+
+    if varn.find(',') != -1:
+        varns = varn.split(',')
+    else:
+        if varn == 'all':
+            varns = objfile.variables
+        else:
+            varns = [varn]
+
+    lonv, latv = lonlat2D(lonobj[:],latobj[:])
+
+    dimx = lonv.shape[1]
+    dimy = lonv.shape[0]
+
+    timv = timobj[:]
+
+# Selecting accordingly a trajectory
+##
+    Ttraj = file_nlines(trajfile,'#')
+    if timekind == 'wrf':
+        dimt = objfile.variables[timn].shape[0] 
+    else:
+        dimt = objfile.variables[timn].shape
+
+    if Tbeg + Ttraj > dimt:
+        print errormsg
+        print '  ' + fname + ': trajectory has ', Ttraj, ' time steps and starts ' + \
+        ' at',Tbeg,' and data',varobj.shape[0], '" !!!!!'
+        quit(-1)
+
+    print '    ' + fname + ': Number of time-steps in trajectory file: ',Ttraj
+
+    trajobj = open(trajfile,'r')
+
+# Trajectory values/slices
+##
+    iline = 0
+    it = 0
+
+    xrangeslice = []
+    yrangeslice = []
+    xrangeslice2D = []
+    yrangeslice2D = []
+
+    cslicev = []
+    cslice2D = []
+    cslicevnoT = []
+
+    gtrajvals = np.zeros((Ttraj,3), dtype=int)
+    trajvals = np.zeros((Ttraj,3), dtype=np.float)
+
+    trajpos = np.zeros(Ttraj,Nangle,raddi)
+
+    it = 0
+    iline = 0
+    for line in trajobj:
+
+## Skipping first line
+##        if not iline == 0:
+##        it = iline - 1
+##        it = iline
+
+# Slicing brings to reduce 1 time-step.... ???
+        if line[0:1] != '#':
+            gtrajvals[it,0] = Tbeg + iline
+            gtrajvals[it,1] = int(line.split(' ')[1])
+            gtrajvals[it,2] = int(line.split(' ')[2])
+#            print it,'t:',gtrajvals[it,0],'x y:', gtrajvals[it,1], gtrajvals[it,2]
+        
+            if timekind == 'wrf':
+                gdate = datetimeStr_conversion(timv[gtrajvals[it,0]],'WRFdatetime',  \
+                  'matYmdHMS')
+                trajvals[it,0] = realdatetime1_CFcompilant(gdate, '19491201000000',  \
+                  'hours')
+                tunits = 'hours since 1949-12-01 00:00:00'
+            elif timekind == 'cf':
+                trajvals[it,0] = timv[gtrajvals[it,0]]
+                tunits = timobj.getncattr('units')
+            else:
+                print errormsg
+                print '  ' + fname + ' time kind: "' + timekind + '" not ready !!'
+                quit(-1)
+
+            trajvals[it,1] = lonv[gtrajvals[it,2],gtrajvals[it,1]]
+            trajvals[it,2] = latv[gtrajvals[it,2],gtrajvals[it,1]]
+
+#            print iline, it,'time:',trajvals[it,0],'lon lat:', trajvals[it,1],       \
+#              trajvals[it,2]
+
+# Assuming time as the first dimension in a fortran like way (t=0 as 1)
+#            trajcut[0] = tstept - 1
+# Assuming time as the first dimension in a C/python like way (t=0)
+
+            varvalst = np.ones((Nangle, radii), dtype=np.float)*fillValue
+            posangle = np.zeros((radii,2), dtype=np.float)
+
+# Angle loop
+            for ia in range(Nangle):
+                angle = 2.*np.pi*ia/Nangle
+
+# Points at that given angle
+                pangle = radial_points(angle,radii)
+                posangle[:,0] = pangle[:,0] + trajvals[it,1]
+                posangle[:,1] = pangle[:,1] + trajvals[it,2]
+
+                for ip in range(radii):
+                    iipos = int(posangle[ip,0])
+                    jjpos = int(posangle[ip,1])
+
+# Creation of a matrix with all the grid points at each time-step
+                    if iipos >= 0 and iipos < dimx and jjpos >= 0 and jjpos < dimy:
+                         print 'Hola'
+                         quit()
+#                        varvalst[ia,ip] = 
+
+
+                if gtrajvals[it,2]-radii < 0: 
+                    yinit = 0
+                    yinit2D = radii - gtrajvals[it,2]
+                else:
+                    yinit = gtrajvals[it,2]-radii
+                    yinit2D = 0
+
+                if gtrajvals[it,2]+radii + 1 > dimy + 1: 
+                    yend = dimy + 1
+                    yend2D = dimy + 1 - gtrajvals[it,2] + radii
+                else:
+                    yend = gtrajvals[it,2]+radii + 1
+                    yend2D = radi
+
+                if gtrajvals[it,1]-box2 < 0: 
+                    xinit = 0
+                    xinit2D = box2 - gtrajvals[it,1]
+                else:
+                    xinit = gtrajvals[it,1]-box2
+                    xinit2D = 0
+
+                if gtrajvals[it,1]+box2 + 1 > dimx + 1: 
+                    xend = dimx + 1
+                    xend2D = dimx + 1 - gtrajvals[it,1] - box2
+                else:
+                    xend = gtrajvals[it,1]+box2 + 1
+                    xend2D = boxs
+
+                yrangeslice.append([yinit, yend])
+                xrangeslice.append([xinit, xend])
+                yrangeslice2D.append([yinit2D, yend2D])
+                xrangeslice2D.append([xinit2D, xend2D])
+            else:
+                yrangeslice.append([gtrajvals[it,2]-box2, gtrajvals[it,2]+box2 + 1])
+                xrangeslice.append([gtrajvals[it,1]-box2, gtrajvals[it,1]+box2 + 1])
+                yrangeslice2D.append([gtrajvals[it,2]-box2, gtrajvals[it,2]+box2 + 1])
+                xrangeslice2D.append([gtrajvals[it,1]-box2, gtrajvals[it,1]+box2 + 1])
+
+# circle values
+            circdist[it,:,:] = radius_dist(dimy,dimx,gtrajvals[it,2],gtrajvals[it,1])
+
+            if gtrajvals[it,2]-Nrad < 0 or gtrajvals[it,2]+Nrad + 1 > dimy + 1 \
+              or gtrajvals[it,1]-Nrad < 0 or gtrajvals[it,1]+Nrad + 1 > dimx + 1:
+
+                if gtrajvals[it,2]-Nrad < 0: 
+                    yinit = 0
+                    yinit2D = Nrad - gtrajvals[it,2]
+                else:
+                    yinit = gtrajvals[it,2]-Nrad
+                    yinit2D = 0
+
+                if gtrajvals[it,2]+Nrad + 1 > dimy + 1: 
+                    yend = dimy + 1
+                    yend2D = dimy + 1 - gtrajvals[it,2] + Nrad
+                else:
+                    yend = gtrajvals[it,2]+Nrad + 1
+                    yend2D = 2*Nrad+1
+
+                if gtrajvals[it,1]-Nrad < 0: 
+                    xinit = 0
+                    xinit2D = Nrad - gtrajvals[it,1]
+                else:
+                    xinit = gtrajvals[it,1]-Nrad
+                    xinit2D = 0
+
+                if gtrajvals[it,1]+Nrad + 1 > dimx + 1: 
+                    xend = dimx + 1
+                    xend2D = dimx + 1 - gtrajvals[it,1] - Nrad
+                else:
+                    xend = gtrajvals[it,1]+Nrad + 1
+                    xend2D = 2*Nrad+1
+
+                cyrangeslice.append([yinit, yend])
+                cxrangeslice.append([xinit, xend])
+                cyrangeslice2D.append([yinit2D, yend2D])
+                cxrangeslice2D.append([xinit2D, xend2D])
+            else:
+                cyrangeslice.append([gtrajvals[it,2]-Nrad, gtrajvals[it,2]+Nrad + 1])
+                cxrangeslice.append([gtrajvals[it,1]-Nrad, gtrajvals[it,1]+Nrad + 1])
+                cyrangeslice2D.append([gtrajvals[it,2]-Nrad, gtrajvals[it,2]+Nrad+1])
+                cxrangeslice2D.append([gtrajvals[it,1]-Nrad, gtrajvals[it,1]+Nrad+1])
+            it = it + 1
+            iline = iline + 1
+
+    trajobj.close()
+# Creation of the netCDF file
+##
+    if varn.find(',') != -1:
+        varnS = 'multi-var'
+    else:
+        varnS = varn.replace(',','-')
+
+    ofile = 'tevolboxtraj_' + varnS + '.nc'
+    objofile = NetCDFFile(ofile, 'w')
+
+# Dimensions
+    newdim = objofile.createDimension('x', boxs)
+    newdim = objofile.createDimension('y', boxs)
+    newdim = objofile.createDimension('xr', Nrad*2+1)
+    newdim = objofile.createDimension('yr', Nrad*2+1)
+    newdim = objofile.createDimension('time', None)
+
+    stsn = ['min', 'max', 'mean', 'mean2', 'stdev', 'ac']
+    vstlname = ['minimum value within', 'maximum value within', 'mean value within', \
+      'squared mean value within', 'standard deviation value within',                \
+      'accumulated value within']
+    Nsts = len(stsn)
+    statnames = []
+    cstatnames = []
+    for i in range(Nsts):
+        statnames.append(stsn[i] + 'box')
+        cstatnames.append(stsn[i] + 'circle')
+
+# Getting values
+##
+    ivar = 0
+
+    for vn in varns:
+        vnst = variables_values(vn)[0]
+        if not searchInlist(objfile.variables, vn):
+            print errormsg
+            print '  ' + fname + ": variable name '" + vn + "' is not in file " +    \
+              ncfile + '" !!!!!'
+            quit(-1)
+
+        varobj = objfile.variables[vn]
+        Nvardims = len(varobj.shape)
+
+        print '  ' + fname + ": getting '" + vn + "' ... .. ." 
+
+        slicev = []
+        slice2D = []
+        slicevnoT = []
+        cslicev = []
+        cslice2D = []
+        cslicevnoT = []
+
+        if Nvardims == 4:
+# Too optimistic, but at this stage...
+            if not objofile.dimensions.has_key('z'):
+                varzobj = objfile.variables[zn]
+                if len(varzobj.shape) == 1:
+                    dimz = varzobj.shape
+                    zvals = varzobj[:]
+                elif len(varzobj.shape) == 2:
+                    dimz = varzobj.shape[1]
+                    zvals = varzobj[0,:]
+                elif len(varzobj.shape) == 3:
+                    dimz = varzobj.shape[2]
+                    zvals = varzobj[0,0,:]
+
+                objofile.createDimension('z', dimz)
+                newvar = objofile.createVariable(zn, 'f4', ('z'),fill_value=fillValue)
+                if searchInlist(varzobj.ncattrs(),'standard_name'):
+                    vsname = varzobj.getncattr('standard_name')
+                else:
+                    vsname = variables_values(zn)[1]
+                if searchInlist(varzobj.ncattrs(),'long_name'):
+                    vlname = varzobj.getncattr('long_name')
+                else:
+                    vlname = variables_values(zn)[4].replace('|',' ')
+                if searchInlist(varzobj.ncattrs(),'units'):
+                    vunits = varzobj.getncattr('units')
+                else:
+                    vunits = variables_values(zn)[5].replace('|',' ')
+    
+                newattr = basicvardef(newvar, vsname, vlname, vunits)
+                newvar[:] = zvals
+ 
+            varvals = np.ones(tuple([Ttraj,dimz,boxs,boxs]), dtype=np.float)
+            lonvals = np.ones(tuple([Ttraj,dimz,boxs,boxs]), dtype=np.float)
+            latvals = np.ones(tuple([Ttraj,boxs,boxs]), dtype=np.float)
+            rvarvals = np.ones(tuple([Ttraj,dimz,Nrad*2+1,Nrad*2+1]), dtype=np.float)
+            rlonvals = np.ones(tuple([Ttraj,dimz,Nrad*2+1,Nrad*2+1]), dtype=np.float)
+            rlatvals = np.ones(tuple([Ttraj,dimz,Nrad*2+1,Nrad*2+1]), dtype=np.float)
+
+# One dimension plus for the values at the center of the trajectory
+            statvarvals = np.ones(tuple([Ttraj,dimz,Nsts+1]), dtype=np.float)
+            rstatvarvals = np.ones(tuple([Ttraj,dimz,Nsts+1]), dtype=np.float)
+
+            for it in range(Ttraj):
+                it0 = Tbeg + it
+
+                slicev = []
+                slice2D = []
+                slicevnoT = []
+                cslicev = []
+                cslice2D = []
+                cslice2Dhor = []
+                cslicevnoT = []
+                cslicevnoThor = []
+
+                slicev.append(gtrajvals[it,0])
+                if gtrajvals[it,2]-box2 < 0 or gtrajvals[it,2]+box2 + 1 > dimy + 1   \
+                  or gtrajvals[it,1]-box2 < 0 or gtrajvals[it,1]+box2 + 1 > dimx +1:
+# box values
+                    slicev.append(slice(0,dimz))
+                    slicev.append(slice(yrangeslice[it][0],yrangeslice[it][1]))
+                    slicev.append(slice(xrangeslice[it][0],xrangeslice[it][1]))
+
+                    slicevnoT.append(slice(0,dimz))
+                    slicevnoT.append(slice(yrangeslice[it][0],yrangeslice[it][1]))
+                    slicevnoT.append(slice(xrangeslice[it][0],xrangeslice[it][1]))
+
+                    slice2D.append(slice(0,dimz))
+                    slice2D.append(slice(0,yrangeslice[it][1]-yrangeslice[it][0]))
+                    slice2D.append(slice(0,xrangeslice[it][1]-xrangeslice[it][0]))
+
+                    rvarvalst = np.ones((dimz, Nrad*2+1, Nrad*2+1),dtype=np.float)*  \
+                      fillValue
+
+                    varvalst[tuple(slice2D)] = varobj[tuple(slicev)]
+                    varvals[it,:,:,:] = varvalst
+    
+# box stats values
+                    maskedvals = ma.masked_values (varvalst, fillValue)
+                    maskedvals2 = maskedvals*maskedvals
+                    for iz in range(dimz):
+                        statvarvals[it,iz,0] = varvalst[iz,box2,box2]
+                        statvarvals[it,iz,1] = np.min(varvalst[iz,:,:])
+                        statvarvals[it,iz,2] = np.max(varvalst[iz,:,:])
+                        statvarvals[it,iz,3] = np.mean(varvalst[iz,:,:])
+                        statvarvals[it,iz,4] = maskedvals2[iz,:,:].mean()
+                        statvarvals[it,iz,5] = np.sqrt(statvarvals[it,iz,4] -        \
+                          statvarvals[it,iz,3]*statvarvals[it,iz,3])
+                        statvarvals[it,iz,6] = np.sum(varvalst[iz,:,:])
+
+                else:
+                    slicev.append(slice(0,dimz))
+                    slicev.append(slice(gtrajvals[it,2]-box2, gtrajvals[it,2]+box2+1))
+                    slicev.append(slice(gtrajvals[it,1]-box2, gtrajvals[it,1]+box2+1))
+                    slicevnoT.append(slice(0,dimz))
+                    slicevnoT.append(slice(gtrajvals[it,2]-box2, gtrajvals[it,2]+    \
+                      box2+1))
+                    slicevnoT.append(slice(gtrajvals[it,1]-box2, gtrajvals[it,1]+    \
+                      box2+1))
+                    slice2D.append(slice(0,dimz))
+                    slice2D.append(slice(gtrajvals[it,2]-box2, gtrajvals[it,2] +     \
+                      box2 + 1))
+                    slice2D.append(slice(gtrajvals[it,1]-box2, gtrajvals[it,1] +     \
+                      box2 + 1))
+
+                    varvalst = varobj[tuple(slicev)]
+# box values
+
+                    varvals[it,:,:,:] = varvalst
+#                    print 'values at time t______'
+#                    print varvalst
+
+# box stats values
+                    for iz in range(dimz):
+                        statvarvals[it,:,0] = varvalst[:,box2,box2]
+                        statvarvals[it,iz,1] = np.min(varvalst[iz,:,:])
+                        statvarvals[it,iz,2] = np.max(varvalst[iz,:,:])
+                        statvarvals[it,iz,3] = np.mean(varvalst[iz,:,:])
+                        statvarvals[it,iz,4] = np.mean(varvalst[iz,:,:]*             \
+                          varvalst[iz,:,:])
+                        statvarvals[it,iz,5] = np.sqrt(statvarvals[it,iz,4] -        \
+                          statvarvals[it,iz,3]*statvarvals[it,iz,3])
+                        statvarvals[it,iz,6] = np.sum(varvalst[iz,:,:])
+
+# Circle values
+                cslicev.append(gtrajvals[it,0])
+                if gtrajvals[it,2]-Nrad < 0 or gtrajvals[it,2]+Nrad + 1 >= dimy       \
+                  or gtrajvals[it,1]-Nrad < 0 or gtrajvals[it,1]+Nrad + 1 >= dimx:
+
+                    maxx = np.min([cxrangeslice[it][1], dimx])
+                    maxy = np.min([cyrangeslice[it][1], dimy])
+                    cslicev.append(slice(0,dimz))
+                    cslicev.append(slice(cyrangeslice[it][0], maxy))
+                    cslicev.append(slice(cxrangeslice[it][0], maxx)) 
+
+                    cslicevnoT.append(slice(0,dimz))
+                    cslicevnoT.append(slice(cyrangeslice[it][0], cyrangeslice[it][1]))
+                    cslicevnoT.append(slice(cxrangeslice[it][0], cxrangeslice[it][1]))
+
+                    cslice2D.append(slice(0,dimz))
+                    cslice2D.append(slice(0,maxy-cyrangeslice[it][0]))
+                    cslice2D.append(slice(0,maxx-cxrangeslice[it][0]))
+                    cslice2Dhor.append(slice(0, maxy - cyrangeslice[it][0]))
+                    cslice2Dhor.append(slice(0, maxx -cxrangeslice[it][0]))
+
+                    rvarvalst = np.ones((dimz,Nrad*2+1,Nrad*2+1),dtype=np.float)*    \
+                      fillValue
+                    rvarvalst[tuple(cslice2D)] = varobj[tuple(cslicev)]
+                    for iz in range(dimz):
+                        tslice = [slice(it,it)]+cslice2Dhor
+                        zslice = [slice(iz,iz)]+cslice2Dhor
+                        rvarvalst[tuple(zslice)] = np.where(circdist[tuple(tslice)] > \
+                          np.float(Nrad), fillValue, rvarvalst[tuple(zslice)]) 
+
+                    rvarvals[it,:,:,:] = rvarvalst
+
+# circle stats values
+                    maskedvals = ma.masked_values (rvarvalst, fillValue)
+                    maskedvals2 = maskedvals*maskedvals
+                    for iz in range(dimz):
+                        rstatvarvals[it,iz,0] = varvalst[iz,box2,box2]
+                        rstatvarvals[it,iz,1] = np.min(varvalst[iz,:,:])
+                        rstatvarvals[it,iz,2] = np.max(varvalst[iz,:,:])
+                        rstatvarvals[it,iz,3] = np.mean(varvalst[iz,:,:])
+                        rstatvarvals[it,iz,4] = maskedvals2[iz,:,:].mean()
+                        rstatvarvals[it,iz,5] = np.sqrt(rstatvarvals[it,iz,4] -      \
+                          rstatvarvals[it,iz,3]*rstatvarvals[it,iz,3])
+                        rstatvarvals[it,iz,6] = np.sum(varvalst[iz,:,:])
+
+                else:
+                    cslicev.append(slice(0,dimz))
+                    cslicev.append(slice(gtrajvals[it,2]-Nrad,gtrajvals[it,2]+Nrad+1))
+                    cslicev.append(slice(gtrajvals[it,1]-Nrad,gtrajvals[it,1]+Nrad+1))
+                    cslicevnoT.append(slice(0,dimz))
+                    cslicevnoT.append(slice(gtrajvals[it,2]-Nrad, gtrajvals[it,2]+   \
+                      Nrad+1))
+                    cslicevnoT.append(slice(gtrajvals[it,1]-Nrad, gtrajvals[it,1]+   \
+                      Nrad+1))
+                    cslicevnoThor.append(slice(gtrajvals[it,2]-Nrad,                 \
+                      gtrajvals[it,2] + Nrad+1))
+                    cslicevnoThor.append(slice(gtrajvals[it,1]-Nrad,                 \
+                      gtrajvals[it,1] + Nrad+1))
+                    cslice2D.append(slice(0,dimz))
+                    cslice2D.append(slice(gtrajvals[it,2]-Nrad,gtrajvals[it,2] +     \
+                      Nrad+1))
+                    cslice2D.append(slice(gtrajvals[it,1]-Nrad,gtrajvals[it,1] +     \
+                      Nrad+1))
+
+                    rvarvalst = varobj[tuple(cslicev)]
+# circle values
+                    for iz in range(dimz):
+                        tslice = [it]+cslicevnoThor
+                        rvarvalst[iz,:,:] = np.where(circdist[tuple(tslice)] >       \
+                          np.float(Nrad), fillValue, rvarvalst[iz,:,:]) 
+
+                    rvarvals[it,:,:,:] = rvarvalst
+
+# circle stats values
+                    maskedvals = ma.masked_values (rvarvalst, fillValue)
+                    maskedvals2 = maskedvals*maskedvals
+                    for iz in range(dimz):
+                        rstatvarvals[it,iz,0] = rvarvalst[iz,Nrad,Nrad]
+                        rstatvarvals[it,iz,1] = maskedvals[iz,:,:].min()
+                        rstatvarvals[it,iz,2] = maskedvals[iz,:,:].max()
+                        rstatvarvals[it,iz,3] = maskedvals[iz,:,:].mean()
+                        rstatvarvals[it,iz,4] = maskedvals2[iz,:,:].mean()
+                        rstatvarvals[it,iz,5] = np.sqrt(rstatvarvals[it,iz,4] -      \
+                          rstatvarvals[it,iz,3]*rstatvarvals[it,iz,3])
+                        rstatvarvals[it,iz,6] = maskedvals[iz,:,:].sum()
+
+#            print 'statistics:',rstatvarvals[it,:] 
+
+# variable box values
+            newvar = objofile.createVariable(vnst + 'box', 'f4', ('time','z','y','x'),\
+              fill_value=fillValue)
+            if searchInlist(varobj.ncattrs(),'standard_name'):
+                vsname = varobj.getncattr('standard_name')
+            else:
+                vsname = variables_values(vn)[1]
+            if searchInlist(varobj.ncattrs(),'long_name'):
+                vlname = varobj.getncattr('long_name')
+            else:
+                vlname = variables_values(vn)[4].replace('|',' ')
+            if searchInlist(varobj.ncattrs(),'units'):
+                vunits = varobj.getncattr('units')
+            else:
+                vunits = variables_values(vn)[5].replace('|',' ')
+
+            newattr = basicvardef(newvar, vsname, vlname, vunits)
+            newattr = newvar.setncattr('projection','lon lat')
+            newvar[:] = varvals
+
+# center of the trajectory
+            newvar = objofile.createVariable('trj_' + vnst, 'f', ('time','z'),     \
+              fill_value=fillValue)
+            newattr = basicvardef(newvar, 'trj_' + vsname, 'value along the ' +    \
+              'trajectory of '+ vn, vunits)
+            newvar[:] = statvarvals[:,:,0]
+
+# variable box statistics
+            ist = 0
+            for statn in statnames:
+                newvar = objofile.createVariable(statn + '_' + vnst,'f',('time','z'),\
+                  fill_value=fillValue)
+                newattr = basicvardef(newvar, statn + '_' + vsname, vstlname[ist] +  \
+                  ' the box ('+str(boxs)+'x'+str(boxs)+') of ' + vnst, vunits)
+                newvar[:] = statvarvals[:,:,ist+1]
+#                newattr = set_attributek(newvar,'_FillValue',fillValue,'npfloat')
+                ist = ist + 1
+        
+# variable circle values
+            newvar = objofile.createVariable(vnst + 'circle','f4',('time','z','yr',    \
+              'xr'), fill_value=fillValue)
+            if searchInlist(varobj.ncattrs(),'standard_name'):
+                vsname = varobj.getncattr('standard_name')
+            else:
+                vsname = variables_values(vn)[1]
+            if searchInlist(varobj.ncattrs(),'long_name'):
+                vlname = varobj.getncattr('long_name')
+            else:
+                vlname = variables_values(vn)[4].replace('|',' ')
+            if searchInlist(varobj.ncattrs(),'units'):
+                vunits = varobj.getncattr('units')
+            else:
+                vunits = variables_values(vn)[5].replace('|',' ')
+
+            newattr = basicvardef(newvar, vsname, vlname, vunits)
+            newattr = newvar.setncattr('projection','lon lat')
+            newvar[:] = rvarvals
+
+# variable circle statistics
+            ist = 0
+            for statn in cstatnames:
+                newvar = objofile.createVariable(statn + '_' + vnst,'f',('time','z'),\
+                  fill_value=fillValue)
+                newattr = basicvardef(newvar, statn + '_' + vsname, vstlname[ist] +  \
+                  ' the circle of radius ('+ str(circler)+') of ' + vnst, vunits)
+                newvar[:] = rstatvarvals[:,:,ist+1]
+#                newattr = set_attributek(newvar,'_FillValue',fillValue,'npfloat')
+                ist = ist + 1
+
+        elif Nvardims == 3:
+# Too optimistic, but at this stage...
+            dimt = varobj.shape[0]
+ 
+            varvals = np.ones(tuple([Ttraj,boxs,boxs]), dtype=np.float)
+            lonvals = np.ones(tuple([Ttraj,boxs,boxs]), dtype=np.float)
+            latvals = np.ones(tuple([Ttraj,boxs,boxs]), dtype=np.float)
+            rvarvals = np.ones(tuple([Ttraj,Nrad*2+1,Nrad*2+1]), dtype=np.float)
+            rlonvals = np.ones(tuple([Ttraj,Nrad*2+1,Nrad*2+1]), dtype=np.float)
+            rlatvals = np.ones(tuple([Ttraj,Nrad*2+1,Nrad*2+1]), dtype=np.float)
+
+# One dimension plus for the values at the center of the trajectory
+            statvarvals = np.ones(tuple([Ttraj,Nsts+1]), dtype=np.float)
+            rstatvarvals = np.ones(tuple([Ttraj,Nsts+1]), dtype=np.float)
+
+            for it in range(Ttraj):
+                it0 = Tbeg + it
+
+                slicev = []
+                slice2D = []
+                slicevnoT = []
+                cslicev = []
+                cslice2D = []
+                cslicevnoT = []
+                cslicevC = []
+
+                slicev.append(gtrajvals[it,0])
+                if gtrajvals[it,2]-box2 < 0 or gtrajvals[it,2]+box2 + 1 > dimy + 1   \
+                  or gtrajvals[it,1]-box2 < 0 or gtrajvals[it,1]+box2 + 1 > dimx + 1:
+# box values
+                    slicev.append(slice(yrangeslice[it][0],yrangeslice[it][1]))
+                    slicev.append(slice(xrangeslice[it][0],xrangeslice[it][1]))
+
+                    slicevnoT.append(slice(yrangeslice[it][0],yrangeslice[it][1]))
+                    slicevnoT.append(slice(xrangeslice[it][0],xrangeslice[it][1]))
+
+                    slice2D.append(slice(0,yrangeslice[it][1]-yrangeslice[it][0]))
+                    slice2D.append(slice(0,xrangeslice[it][1]-xrangeslice[it][0]))
+
+                    rvarvalst = np.ones((Nrad*2+1, Nrad*2+1),dtype=np.float)*fillValue
+
+                    varvalst[tuple(slice2D)] = varobj[tuple(slicev)]
+                    varvals[it,:,:] = varvalst
+    
+# box stats values
+                    maskedvals = ma.masked_values (varvalst, fillValue)
+                    statvarvals[it,0] = varvalst[box2,box2]
+                    statvarvals[it,1] = maskedvals.min()
+                    statvarvals[it,2] = maskedvals.max()
+                    statvarvals[it,3] = maskedvals.mean()
+                    maskedvals2 = maskedvals*maskedvals
+                    statvarvals[it,4] = maskedvals2.mean()
+                    statvarvals[it,5] = np.sqrt(statvarvals[it,4] -                   \
+                      statvarvals[it,3]*statvarvals[it,3])
+                    statvarvals[it,6] = maskedvals.sum()
+
+                    varvalst[tuple(slice2D)] = lonv[tuple(slicevnoT)]
+                    lonvals[it,:,:] = varvalst
+                    varvalst[tuple(slice2D)] = latv[tuple(slicevnoT)]           
+                    latvals[it,:,:] = varvalst
+
+                else:
+                    slicev.append(slice(gtrajvals[it,2]-box2, gtrajvals[it,2]+box2+1))
+                    slicev.append(slice(gtrajvals[it,1]-box2, gtrajvals[it,1]+box2+1))
+                    slicevnoT.append(slice(gtrajvals[it,2]-box2, gtrajvals[it,2]+    \
+                      box2+1))
+                    slicevnoT.append(slice(gtrajvals[it,1]-box2, gtrajvals[it,1]+    \
+                      box2+1))
+
+                    slice2D.append(slice(gtrajvals[it,2]-box2, gtrajvals[it,2]+box2 +\
+                      1))
+                    slice2D.append(slice(gtrajvals[it,1]-box2, gtrajvals[it,1]+box2 +\
+                      1))
+
+                    varvalst = varobj[tuple(slicev)]
+# box values
+
+                    varvals[it,:,:] = varvalst
+#                    print 'values at time t______'
+#                    print varvalst
+
+# box stats values
+                    statvarvals[it,0] = varvalst[box2,box2]
+                    statvarvals[it,1] = np.min(varvalst)
+                    statvarvals[it,2] = np.max(varvalst)
+                    statvarvals[it,3] = np.mean(varvalst)
+                    statvarvals[it,4] = np.mean(varvalst*varvalst)
+                    statvarvals[it,5] = np.sqrt(statvarvals[it,4] -                  \
+                      statvarvals[it,3]*statvarvals[it,3])
+                    statvarvals[it,6] = np.sum(varvalst)
+
+                    varvalst = lonv[tuple(slicevnoT)]
+                    lonvals[it,:,:] = varvalst
+                    varvalst = latv[tuple(slicevnoT)]           
+                    latvals[it,:,:] = varvalst
+
+# Circle values
+                cslicev.append(gtrajvals[it,0])
+                if gtrajvals[it,2]-Nrad < 0 or gtrajvals[it,2]+Nrad + 1 >= dimy      \
+                  or gtrajvals[it,1]-Nrad < 0 or gtrajvals[it,1]+Nrad + 1 >= dimx:
+                    maxx = np.min([cxrangeslice[it][1], dimx])
+                    maxy = np.min([cyrangeslice[it][1], dimy])
+                    cslicev.append(slice(cyrangeslice[it][0],maxy))
+                    cslicev.append(slice(cxrangeslice[it][0],maxx))
+
+                    cslicevnoT.append(slice(cyrangeslice[it][0],cyrangeslice[it][1]))
+                    cslicevnoT.append(slice(cxrangeslice[it][0],cxrangeslice[it][1]))
+
+                    cslice2D.append(slice(0,maxy - cyrangeslice[it][0]))
+                    cslice2D.append(slice(0,maxx - cxrangeslice[it][0]))
+
+                    rvarvalst = np.ones((Nrad*2+1,Nrad*2+1),dtype=np.float)*fillValue
+                    rvarvalst[tuple(cslice2D)] = varobj[tuple(cslicev)]
+                    rvarvalst[tuple(cslice2D)] = np.where(circdist[tuple(cslicev)] >\
+                      np.float(Nrad), fillValue, rvarvalst[tuple(cslice2D)])
+
+                    rvarvals[it,:,:] = rvarvalst
+
+# circle stats values
+                    maskedvals = ma.masked_values (rvarvalst, fillValue)
+                    rstatvarvals[it,0] = rvarvalst[Nrad,Nrad]
+                    rstatvarvals[it,1] = maskedvals.min()
+                    rstatvarvals[it,2] = maskedvals.max()
+                    rstatvarvals[it,3] = maskedvals.mean()
+                    maskedvals2 = maskedvals*maskedvals
+                    rstatvarvals[it,4] = maskedvals2.mean()
+                    rstatvarvals[it,5] = np.sqrt(rstatvarvals[it,4] -                \
+                      rstatvarvals[it,3]*rstatvarvals[it,3])
+                    rstatvarvals[it,6] = maskedvals2.sum()
+
+                else:
+                    cslicev.append(slice(gtrajvals[it,2]-Nrad,gtrajvals[it,2]+Nrad+1))
+                    cslicev.append(slice(gtrajvals[it,1]-Nrad,gtrajvals[it,1]+Nrad+1))
+                    cslicevnoT.append(slice(gtrajvals[it,2]-Nrad, gtrajvals[it,2]+    \
+                      Nrad+1))
+                    cslicevnoT.append(slice(gtrajvals[it,1]-Nrad, gtrajvals[it,1]+    \
+                      Nrad+1))
+
+                    cslice2D.append(slice(gtrajvals[it,2]-Nrad,gtrajvals[it,2]+Nrad+1))
+                    cslice2D.append(slice(gtrajvals[it,1]-Nrad,gtrajvals[it,1]+Nrad+1))
+
+                    cslicevC.append(it)
+                    cslicevC.append(slice(gtrajvals[it,2]-Nrad,gtrajvals[it,2]+Nrad+1))
+                    cslicevC.append(slice(gtrajvals[it,1]-Nrad,gtrajvals[it,1]+Nrad+1))
+
+                    rvarvalst = varobj[tuple(cslicev)]
+                    cdist = circdist[tuple(cslicevC)]
+# circle values
+                    rvarvalst = np.where(cdist > np.float(Nrad),fillValue,rvarvalst)
+                    rvarvals[it,:,:] = rvarvalst
+
+# circle stats values
+                    maskedvals = ma.masked_values (rvarvalst, fillValue)
+                    rstatvarvals[it,0] = rvarvalst[Nrad,Nrad]
+                    rstatvarvals[it,1] = maskedvals.min()
+                    rstatvarvals[it,2] = maskedvals.max()
+                    rstatvarvals[it,3] = maskedvals.mean()
+                    maskedvals2 = maskedvals*maskedvals
+                    rstatvarvals[it,4] = maskedvals2.mean()
+                    rstatvarvals[it,5] = np.sqrt(rstatvarvals[it,4] -                \
+                      rstatvarvals[it,3]*rstatvarvals[it,3])
+                    rstatvarvals[it,6] = maskedvals.sum()
+
+#            print 'statistics:',rstatvarvals[it,:] 
+
+# variable box values
+            newvar = objofile.createVariable(vnst + 'box', 'f4', ('time', 'y', 'x'), \
+              fill_value=fillValue)
+            if searchInlist(varobj.ncattrs(),'standard_name'):
+                vsname = varobj.getncattr('standard_name')
+            else:
+                vsname = variables_values(vn)[1]
+            if searchInlist(varobj.ncattrs(),'long_name'):
+                vlname = varobj.getncattr('long_name')
+            else:
+                vlname = variables_values(vn)[4].replace('|',' ')
+            if searchInlist(varobj.ncattrs(),'units'):
+                vunits = varobj.getncattr('units')
+            else:
+                vunits = variables_values(vn)[5].replace('|',' ')
+
+            newattr = basicvardef(newvar, vsname, vlname, vunits)
+            newattr = newvar.setncattr('projection','lon lat')
+            newvar[:] = varvals
+
+# center of the trajectory
+            newvar = objofile.createVariable('trj_' + vnst, 'f', ('time'),           \
+              fill_value=fillValue)
+            newattr = basicvardef(newvar, 'trj_' + vsname, 'value along the ' +      \
+              'trajectory of '+ vnst, vunits)
+            newvar[:] = statvarvals[:,0]
+
+# variable box statistics
+            ist = 0
+            for statn in statnames:
+                newvar = objofile.createVariable(statn + '_' + vnst, 'f', ('time'),  \
+                  fill_value=fillValue)
+                newattr = basicvardef(newvar, statn + '_' + vsname, vstlname[ist] +  \
+                  ' the box ('+str(boxs)+'x'+str(boxs)+') of ' + vnst, vunits)
+                newvar[:] = statvarvals[:,ist+1]
+#                newattr = set_attributek(newvar,'_FillValue',fillValue,'npfloat')
+                ist = ist + 1
+        
+# variable circle values
+            newvar = objofile.createVariable(vnst + 'circle','f4',('time','yr','xr'),\
+              fill_value=fillValue)
+            if searchInlist(varobj.ncattrs(),'standard_name'):
+                vsname = varobj.getncattr('standard_name')
+            else:
+                vsname = variables_values(vn)[1]
+            if searchInlist(varobj.ncattrs(),'long_name'):
+                vlname = varobj.getncattr('long_name')
+            else:
+                vlname = variables_values(vn)[4].replace('|',' ')
+            if searchInlist(varobj.ncattrs(),'units'):
+                vunits = varobj.getncattr('units')
+            else:
+                vunits = variables_values(vn)[5].replace('|',' ')
+
+            newattr = basicvardef(newvar, vsname, vlname, vunits)
+            newattr = newvar.setncattr('projection','lon lat')
+            newvar[:] = rvarvals
+
+# variable circle statistics
+            ist = 0
+            for statn in cstatnames:
+                newvar = objofile.createVariable(statn + '_' + vnst, 'f', ('time'),  \
+                  fill_value=fillValue)
+                newattr = basicvardef(newvar, statn + '_' + vsname, vstlname[ist] +  \
+                  ' the circle of radius ('+ str(circler)+') of ' + vnst, vunits)
+                newvar[:] = rstatvarvals[:,ist+1]
+#                newattr = set_attributek(newvar,'_FillValue',fillValue,'npfloat')
+                ist = ist + 1
+        else:
+# Other variables
+            print warnmsg
+            print '  ' + fname + ": variable '" + vn + "' shape:",varobj.shape,' not'\
+              + ' ready!!'
+            print '    skipping variable'
+            if len(varns) == 1:
+                objofile.close()
+                sub.call(['rm', ofile])
+                print '    uniq variable! removing file and finishing program'
+                quit()
+
+        if not objofile.variables.has_key('trlon') and Nvardims == 3:
+# var dimensions
+            newvar = objofile.createVariable('trlon', 'f8', ('time'))
+            newattr = basicvardef(newvar,'trlon','trajectory longitude',             \
+              'degrees west_east')
+            newvar[:] = trajvals[:,1]
+
+            newvar = objofile.createVariable('trlat', 'f8', ('time'))
+            newattr = basicvardef(newvar,'trlat','trajectory latitude',              \
+              'degrees north_south')
+            newvar[:] = trajvals[:,2]
+
+            newvar = objofile.createVariable('time', 'f8', ('time'))
+            newattr = basicvardef(newvar, 'time', 'time', tunits)
+            newvar[:] = trajvals[:,0]
+
+            newvar = objofile.createVariable('lon', 'f8', ('time', 'y', 'x'),        \
+              fill_value=fillValue)
+            newattr = basicvardef(newvar, 'longitude', 'longitude',                  \
+              'degrees west_east')
+            newvar[:] = lonvals
+
+            newvar = objofile.createVariable('lat', 'f8', ('time', 'y', 'x'),        \
+              fill_value=fillValue)
+            newattr = basicvardef(newvar, 'latitude', 'latitude',                    \
+              'degrees north_south')
+            newvar[:] = latvals
+           
+        ivar = ivar + 1
+
+
+# global attributes
+    objofile.setncattr('author', 'L. Fita')
+    objofile.setncattr('institution', 'Laboratire de Meteorologie Dynamique')
+    objofile.setncattr('university', 'Pierre Marie Curie - Jussieu')
+    objofile.setncattr('center', 'Centre National de Recherches Scientifiques')
+    objofile.setncattr('city', 'Paris')
+    objofile.setncattr('country', 'France')
+    objofile.setncattr('script', 'nc_var_tools.py')
+    objofile.setncattr('function', 'compute_tevolboxtraj')
+    objofile.setncattr('version', '1.0')
+    objofile.setncattr('data_file',ncfile)
+
+    objfile.close()
+
+    objofile.sync()
+    objofile.close()
+
+    print '  ' + fname + ': successful creation of file "' + ofile + '" !!!'
+
+    return 
+
 
 #quit()
Index: trunk/tools/read_ISD.py
===================================================================
--- trunk/tools/read_ISD.py	(revision 458)
+++ trunk/tools/read_ISD.py	(revision 458)
@@ -0,0 +1,718 @@
+# -*- coding: iso-8859-15 -*-
+# Python reader of ISD files
+## e.g. # read_ISD.py -c '#' -f 722317-13970-2005 -d description.dat -g true -t 19491201000000,seconds -s 'Baton Rouge METRO R,30.537,-91.147,23.2'
+import numpy as np
+from optparse import OptionParser
+import os
+from netCDF4 import Dataset as NetCDFFile
+
+main = 'read_ISD.py'
+errormsg = 'ERROR -- error -- ERROR -- error'
+warnmsg = 'WARNING -- warning -- WARNING -- warning'
+
+# version
+version=1.1
+
+# Filling values for floats, integer and string
+fillValueF = 1.e20
+fillValueI = -99999
+fillValueS = '---'
+
+# Length of the string variables
+StringLength = 200
+
+####### Different fixed tabuled values
+
+quality_vals = ['Passed gross limits check', 'Passed all quality control checks',    \
+  'Suspect', 'Erroneous',                                                            \
+  'Passed gross limits check , data originate from an NCDC data source',             \
+  'Passed all quality control checks, data originate from an NCDC data source',      \
+  'Suspect, data originate from an NCDC data source',                                \
+  'Erroneous, data originate from an NCDC data source',                              \
+  'Passed gross limits check if element is present' ]
+
+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 remove_NONascii(string):
+    """ Function to remove that characters which are not in the standard 127 ASCII
+      string= string to transform
+    >>> remove_NONascii('LluÃ­s')
+    Lluis
+    """
+    fname = 'remove_NONascii'
+
+    newstring = string
+
+    RTFchar= ['Ã¡', 'Ã©', 'Ã­', 'Ã³', 'Ãº', 'Ã ', 'Ãš', 'Ã¬', 'Ã²', 'Ã¹', 'Ã¢', 'Ãª', 'Ã®', 'ÃŽ',  \
+      'Ã»', 'Ã€', 'Ã«', 'Ã¯', 'Ã¶', 'ÃŒ', 'Ã§', 'Ã±','ÃŠ', 'Å', 'Ã', 'Ã', 'Ã', 'Ã', 'Ã', 'Ã', \
+      'Ã', 'Ã', 'Ã', 'Ã', 'Ã', 'Ã', 'Ã', 'Ã', 'Ã', 'Ã', 'Ã', 'Ã', 'Ã', 'Ã', 'Ã', 'Ã',\
+      'Ã', 'Å', '\n', '\t']
+    ASCchar= ['a', 'e', 'i', 'o', 'u', 'a', 'e', 'i', 'o', 'u', 'a', 'e', 'i', 'o',  \
+      'u', 'a', 'e', 'i', 'o', 'u', 'c', 'n','ae', 'oe', 'A', 'E', 'I', 'O', 'U', 'A', \
+      'E', 'I', 'O', 'U', 'A', 'E', 'I', 'O', 'U', 'A', 'E', 'I', 'O', 'U', 'C', 'N',\
+      'AE', 'OE', '', ' ']
+
+    Nchars = len(RTFchar)
+    for ichar in range(Nchars):
+        foundchar = string.find(RTFchar[ichar])
+        if foundchar != -1:
+            newstring = newstring.replace(RTFchar[ichar], ASCchar[ichar])
+
+    return newstring
+
+def read_description(fdobs, dbg):
+    """ reads the description file of the observational data-set 
+    fdobs= descriptive observational data-set
+    dbg= boolean argument for debugging
+    * Each station should have a 'description.dat' file with:
+      institution=Institution who creates the data
+      department=Department within the institution
+      scientists=names of the data producers
+      contact=contact of the data producers
+      description=description of the observations
+      acknowledgement=sentence of acknowlegement
+      comment=comment for the measurements
+
+      MissingValue='|' list of ASCII values for missing values within the data 
+        (as they appear!)
+      comment=comments
+
+      varN='|' list of variable names
+      varLN='|' list of long variable names
+      varU='|' list units of the variables
+      varBUFR='|' list BUFR code of the variables
+      varTYPE='|' list of variable types ('D', 'F', 'I', 'I64', 'S')
+      varL='|' list of length of the sections for each variable within a given line (Fortran-like)
+
+      NAMElon=name of the variable with the longitude (x position)
+      NAMElat=name of the variable with the latitude (y position)
+      NAMEheight=ind_alt
+      NAMEtime=name of the varibale with the time
+      FMTtime=format of the time (as in 'C', 'CFtime' for already CF-like time)
+
+    """
+    fname = 'read_description'
+
+    descobj = open(fdobs, 'r')
+    desc = {}
+
+    namevalues = []
+
+    for line in descobj:
+        if line[0:1] != '#' and len(line) > 1 :
+            descn = remove_NONascii(line.split('=')[0])
+            descv = remove_NONascii(line.split('=')[1])
+            namevalues.append(descn)
+            if descn[0:3] != 'var':
+                if descn != 'MissingValue':
+                    desc[descn] = descv
+                else:
+                    desc[descn] = []
+                    for dn in descv.split('|'): desc[descn].append(dn)
+                    print '  ' + fname + ': missing values found:',desc[descn]
+            elif descn[0:3] == 'var':
+                desc[descn] = descv.split('|')
+            elif descn[0:4] == 'NAME':
+                desc[descn] = descv
+            elif descn[0:3] == 'FMT':
+                desc[descn] = descv
+
+    if dbg:
+        Nvars = len(desc['varN'])
+        print '  ' + fname + ": description content of '" + fdobs + "'______________"
+        for varn in namevalues:
+            if varn[0:3] != 'var':
+                print '    ' + varn + ':',desc[varn]
+            elif varn == 'varN':
+                print '    * Variables:'
+                for ivar in range(Nvars):
+                    varname = desc['varN'][ivar]
+                    varLname = desc['varLN'][ivar]
+                    varunits = desc['varU'][ivar]
+                    varsec = int(desc['varL'][ivar])
+                    if desc.has_key('varBUFR'): 
+                        varbufr=[ivar]
+                        print '      ', ivar, varname + ':',varLname,'[',varunits,   \
+                          ']','bufr code:',varbufr, 'L:', varsec
+                    else:
+                        print '      ', ivar, varname + ':',varLname,'[',varunits,   \
+                          '] L:',varsec
+
+    descobj.close()
+
+    return desc
+
+def value_fmt(val, miss, fmt):
+    """ Function to transform an ASCII value to a given format
+      val= value to transform
+      miss= list of possible missing values
+      fmt= format to take:
+        'D': float double precission
+        'F': float
+        'I': integer
+        'I64': 64-bits integer
+        'S': string
+    >>> value_fmt('9876.12', '-999', 'F')
+    9876.12
+    """
+
+    fname = 'value_fmt'
+
+    fmts = ['D', 'F', 'I', 'I64', 'S']
+    Nfmts = len(fmts)
+
+    if not searchInlist(miss,val):
+        if not searchInlist(fmts, fmt):
+            print errormsg
+            print '  ' + fname + ": format '" + fmt + "' not ready !!"
+            quit(-1)
+        else:
+            if fmt == 'D':
+                newval = np.float32(val)
+            elif fmt == 'F':
+                newval = np.float(val)
+            elif fmt == 'I':
+                newval = int(val)
+            elif fmt == 'I64':
+                newval = np.int64(val)
+            elif fmt == 'S':
+                newval = val
+    else:
+        newval = None
+
+    return newval
+
+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 Stringtimes_CF(tvals, fmt, Srefdate, tunits, dbg):
+    """ Function to transform a given data in String formats to a CF date
+      tvals= string temporal values
+      fmt= format of the the time values
+      Srefdate= reference date in [YYYY][MM][DD][HH][MI][SS] format
+      tunits= units to use ('weeks', 'days', 'hours', 'minutes', 'seconds') 
+      dbg= debug
+    >>> Stringtimes_CF(['19760217082712','20150213101837'], '%Y%m%d%H%M%S', 
+      '19491201000000', 'hours', False)
+    [229784.45333333  571570.31027778]
+    """
+    import datetime as dt
+
+    fname = 'Stringtimes'
+
+    dimt = len(tvals)
+
+    yrref = int(Srefdate[0:4])
+    monref = int(Srefdate[4:6])
+    dayref = int(Srefdate[6:8])
+    horref = int(Srefdate[8:10])
+    minref = int(Srefdate[10:12])
+    secref = int(Srefdate[12:14])
+    refdate = dt.datetime( yrref, monref, dayref, horref, minref, secref)
+
+    cftimes = np.zeros((dimt), dtype=np.float)
+
+    Nfmt=len(fmt.split('%'))
+
+    if dbg: print '  ' + fname + ': fmt=',fmt,'refdate:',Srefdate,'uits:',tunits,    \
+      'date dt_days dt_time deltaseconds CFtime _______'
+    for it in range(dimt):
+
+# Removing excess of mili-seconds (up to 6 decimals)
+        if fmt.split('%')[Nfmt-1] == 'f':
+            tpoints = tvals[it].split('.')
+            if len(tpoints[len(tpoints)-1]) > 6:
+                milisec = '{0:.6f}'.format(np.float('0.'+tpoints[len(tpoints)-1]))[0:7]
+                newtval = ''
+                for ipt in range(len(tpoints)-1):
+                    newtval = newtval + tpoints[ipt] + '.'
+                newtval = newtval + str(milisec)[2:len(milisec)+1]
+            else:
+                newtval = tvals[it]
+            tval = dt.datetime.strptime(newtval, fmt)
+        else:
+            tval = dt.datetime.strptime(tvals[it], fmt)
+
+        deltatime = tval - refdate
+        deltaseconds = deltatime.days*24.*3600. + deltatime.seconds +                \
+          deltatime.microseconds/100000.
+        if tunits == 'weeks':
+            deltat = 7.*24.*3600.
+        elif tunits == 'days':
+            deltat = 24.*3600.
+        elif tunits == 'hours':
+            deltat = 3600.
+        elif tunits == 'minutes':
+            deltat = 60.
+        elif tunits == 'seconds':
+            deltat = 1.
+        else:
+            print errormsg
+            print '  ' + fname + ": time units '" + tunits + "' not ready !!"
+            quit(-1)
+
+        cftimes[it] = deltaseconds / deltat
+        if dbg:
+            print '  ' + tvals[it], deltatime, deltaseconds, cftimes[it]
+
+    return cftimes
+
+def additional_vars(string):
+    """ Function for the additional variables
+      string: sectino of the file line which has to be evaluated
+    """
+    fname = 'additional_vars'
+
+    Lstring=len(string)
+
+# Generic list of additional variables '-' for range of values
+    addvars0 = ['AA1-4', 'AB1', 'AC1', 'AD1', 'AE1', 'AG1', 'AH1-6', 'AI1-6', 'AJ1', \
+      'AK1', 'AL1-4', 'AM1', 'AN1', 'AO1-4', 'AP1-4', 'HPD', 'AU1-9', 'AW1-4',       \
+      'AX1-6', 'AY1-2', 'AZ1-2', 'CB1-2', 'CF1-3', 'CG1-3', 'CH1-2', 'CI1', 'CN1-4', \
+      'CO1-9', 'CR1', 'CT1-3', 'CU1-3', 'CV1-3', 'CW1-3', 'CX1-3', 'ED1', 'GA1-6',   \
+      'GE1', 'GF1', 'GG1-6', 'GH1', 'GJ1', 'GK1', 'GM1', 'GN1', 'GO1', 'GP1', 'GQ1', \
+      'GR1', 'HL1', 'IA1-2', 'IB1-2', 'IC1', 'KA1-4', 'KB1-3', 'KC1-2', 'KD1-2',     \
+      'KE1', 'KF1', 'KG1-2', 'MA1', 'MD1', 'ME1', 'MF1', 'MG1', 'MH1', 'MK1',        \
+      'MV1-7', 'MW1-7', 'OA1-3', 'OB1-2', 'OC1', 'OD1-3', 'OE1-3', 'RH1-3', 'SA1',   \
+      'ST1', 'UA1', 'UG1-2', 'WA1', 'WD1', 'WG1', 'WJ1', 'REM', 'EQD', 'N01-99',     \
+      'QNN']
+
+    addvars = []
+    for avar in addvars0:
+        if avar.find('-') == -1:
+            addvars.append(avar)
+        else:
+            totrange=avar.split('-')
+            erange = totrange[1]
+
+            if len(erange) > 1:
+                irange = totrange[0][1:3]
+            else:
+                irange = totrange[0][2:3]
+          
+            for ir in range(int(irange), int(erange)+1):
+                if len(erange) > 1:
+                    addvars.append(avar[0:1] + str(ir).zfill(2))
+                else:
+                    addvars.append(avar[0:2] + str(ir))
+
+#    print addvars
+# Looking for the variables
+    ichar=0
+
+    additionals = []
+    while ichar <= Lstring:
+        sec3 = string[ichar:ichar+3]
+        if searchInlist(addvars, sec3):
+            additionals.append(sec3)
+#            print '  ' , sec3, 'Additional!!'
+            ichar = ichar + 2
+
+        ichar = ichar + 1
+
+    return additionals
+
+def read_datavalues_conline(dataf, comchar, fmt, Lv, miss, varns, dbg):
+    """ Function to read from an ASCII file values in ISH format (continuos line)
+      dataf= data file
+      comchar= list of the characters indicating comment in the file
+      dbg= debug mode or not
+      fmt= list of kind of values to be found
+      Lv= length of the value section
+      miss= missing value
+      varns= list of name of the variables to find
+    """
+    fname = 'read_datavalues_conline'
+
+    ofile = open(dataf, 'r')
+    Nvals = len(fmt)
+
+    finalvalues = {}
+
+    iline = 0
+    for line in ofile:
+        line = line.replace('\n','').replace(chr(13),'')
+        Lline = len(line)
+        if not searchInlist(comchar,line[0:1]) and len(line) > 1:
+# Removing no-value columns
+            values = []
+
+            for ivar in range(Nvals):
+                if ivar == 0:
+                    isec = 0
+                else:
+                    isec = int(Lv[ivar-1])
+
+                esec = int(Lv[ivar])
+                val = line[isec:esec]
+
+                values.append(val)
+                if dbg: 
+                    print iline, varns[ivar],'value:',values[ivar],miss,fmt[ivar]
+
+                if iline == 0:
+                    listvals = []
+                    listvals.append(value_fmt(values[ivar], miss, fmt[ivar]))
+                    finalvalues[varns[ivar]] = listvals
+                else:
+                    listvals = finalvalues[varns[ivar]]
+                    listvals.append(value_fmt(values[ivar], miss, fmt[ivar]))
+                    finalvalues[varns[ivar]] = listvals
+# Additional variables
+            if Lline > esec and line[esec:esec+3] == 'ADD':
+                if iline == 0: print '  ' + fname + ': Additional variables!!'
+                ichar=esec+3
+                print '    ', additional_vars(line[esec+4:Lline+1])
+
+        else:
+# First line without values
+            if iline == 0: iline = -1
+    
+        iline = iline + 1
+
+    ofile.close()
+
+    return finalvalues
+
+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 = objfile.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 = objfile.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 = objfile.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 = objfile.createVariable( heightname, 'f4', ('nst'))
+    basicvardef(newvar, heightname, 'height above sea level', 'm' )
+    newvar[:] = stdesc[3]
+
+    return
+
+####### ###### ##### #### ### ## #
+
+strCFt="Refdate,tunits (CF reference date [YYYY][MM][DD][HH][MI][SS] format and " +  \
+  " and time units: 'weeks', 'days', 'hours', 'miuntes', 'seconds')"
+
+parser = OptionParser()
+parser.add_option("-c", "--comments", dest="charcom",
+  help="':', list of characters used for comments", metavar="VALUES")
+parser.add_option("-d", "--descriptionfile", dest="fdesc",
+  help="description file to use", metavar="FILE")
+parser.add_option("-f", "--file", dest="obsfile",
+  help="observational file to use", metavar="FILE")
+parser.add_option("-g", "--debug", dest="debug",
+  help="whther debug is required ('false', 'true')", metavar="VALUE")
+parser.add_option("-s", "--stationLocation", dest="stloc",  
+  help="',' list with Name, longitude, latitude and height of the station", metavar="VALUES")
+parser.add_option("-t", "--CFtime", dest="CFtime", help=strCFt, metavar="VALUE")
+(opts, args) = parser.parse_args()
+
+#######    #######
+## MAIN
+    #######
+
+ofile='ISD.nc'
+
+if opts.charcom is None:
+    print warnmsg
+    print '  ' + main + ': No list of comment characters provided!!'
+    print '    assuming no need!'
+    charcomments = []
+else:
+    charcomments = opts.charcom.split(':')    
+
+if opts.fdesc is None:
+    print errormsg
+    print '  ' + main + ': No description file for the observtional data provided!!'
+    quit(-1)
+
+if opts.obsfile is None:
+    print errormsg
+    print '  ' + main + ': ISD file not provided !!'
+    quit(-1)
+
+if opts.debug is None:
+    print warnmsg
+    print '  ' + main + ': No debug provided!!'
+    print "    assuming 'False'"
+    debug = False
+else:
+    if opts.debug == 'true':
+        debug = True
+    else:
+        debug = False
+
+if not os.path.isfile(opts.fdesc):
+    print errormsg
+    print '   ' + main + ": description file '" + opts.fdesc + "' does not exist !!"
+    quit(-1)
+
+if not os.path.isfile(opts.obsfile):
+    print errormsg
+    print '   ' + main + ": observational file '" + opts.obsfile + "' does not exist !!"
+    quit(-1)
+
+if opts.CFtime is None:
+    print warnmsg
+    print '  ' + main + ': No CFtime criteria are provided !!'
+    print "    either time is already in CF-format ('timeFMT=CFtime') in '" +        \
+      opts.fdesc + "'"
+    print "    or assuming refdate: '19491201000000' and time units: 'hours'"
+    referencedate = '19491201000000'
+    timeunits = 'hours'
+else:
+    referencedate = opts.CFtime.split(',')[0]
+    timeunits = opts.CFtime.split(',')[1]
+
+if opts.stloc is None:
+    print errornmsg
+    print '  ' + main + ': No station location provided !!'
+    quit(-1)
+else:
+    stationdesc = [opts.stloc.split(',')[0], np.float(opts.stloc.split(',')[1]),     \
+      np.float(opts.stloc.split(',')[2]), np.float(opts.stloc.split(',')[3])]
+
+# Reading description file
+##
+description = read_description(opts.fdesc, debug)
+
+Nvariables=len(description['varN'])
+formats = description['varTYPE']
+
+if len(formats) != Nvariables: 
+    print errormsg
+    print '  ' + main + ': number of formats:',len(formats),' and number of ' +     \
+      'variables', Nvariables,' does not coincide!!'
+    print '    * what is found is _______'
+    if Nvariables > len(formats):
+        Nshow = len(formats)
+        for ivar in range(Nshow):
+            print '      ',description['varN'][ivar],':', description['varLN'][ivar],\
+               '[', description['varU'][ivar], '] fmt:', formats[ivar]
+        print '      missing values for:', description['varN'][Nshow:Nvariables+1]
+    else:
+        Nshow = Nvariables
+        for ivar in range(Nshow):
+            print '      ',description['varN'][ivar],':', description['varLN'][ivar],\
+               '[', description['varU'][ivar], '] fmt:', formats[ivar]
+        print '      excess of formats for:', formats[Nshow:len(formats)+1]
+
+    quit(-1)
+
+# Reading values
+##
+datavalues = read_datavalues_conline(opts.obsfile, charcomments, formats,            \
+  description['varL'], description['MissingValue'], description['varN'], debug)
+
+# Total number of values
+Ntvalues = len(datavalues[description['varN'][0]])
+print main + ': total temporal values found:',Ntvalues
+
+objfile = NetCDFFile(ofile, 'w')
+
+# Creation of dimensions
+##
+objfile.createDimension('time',None)
+objfile.createDimension('StrLength',StringLength)
+
+# Creation of variables
+##
+for ivar in range(Nvariables):
+    varn = description['varN'][ivar]
+    print "  including: '" + varn + "' ... .. ."
+
+    if formats[ivar] == 'D':
+        newvar = objfile.createVariable(varn, 'f32', ('time'), fill_value=fillValueF)
+        basicvardef(newvar, varn, description['varLN'][ivar],                        \
+          description['varU'][ivar])
+        newvar[:] = np.where(datavalues[varn] is None, fillValueF, datavalues[varn])
+    elif formats[ivar] == 'F':
+        newvar = objfile.createVariable(varn, 'f', ('time'), fill_value=fillValueF)
+        basicvardef(newvar, varn, description['varLN'][ivar],                        \
+          description['varU'][ivar])
+        newvar[:] = np.where(datavalues[varn] is None, fillValueF, datavalues[varn])
+    elif formats[ivar] == 'I':
+        newvar = objfile.createVariable(varn, 'i', ('time'), fill_value=fillValueI)
+        basicvardef(newvar, varn, description['varLN'][ivar],                        \
+          description['varU'][ivar])
+# Why is not wotking with integers?
+        vals = np.array(datavalues[varn])
+        for iv in range(Ntvalues):
+            if vals[iv] is None: vals[iv] = fillValueI
+        newvar[:] = vals
+    elif formats[ivar] == 'S':
+        newvar = objfile.createVariable(varn, 'c', ('time','StrLength'))
+        basicvardef(newvar, varn, description['varLN'][ivar],                        \
+          description['varU'][ivar])
+        writing_str_nc(newvar, datavalues[varn], StringLength)
+
+# Time variable in CF format
+##
+if description['FMTtime'] == 'CFtime':
+    timevals = datavalues[description['NAMEtime']]
+    iv = 0
+    for ivn in description['varN']:
+        if ivn == description['NAMEtime']:
+            tunits = description['varU'][iv]
+            break
+        iv = iv + 1
+else:
+    # Time as a composition of different columns
+    tcomposite = description['NAMEtime'].find('@')
+    if tcomposite != -1:
+        timevars = description['NAMEtime'].split('@')
+
+        print warnmsg
+        print '  ' + main + ': time values as combination of different columns!'
+        print '    combining:',timevars,' with a final format: ',description['FMTtime']
+
+        timeSvals = []
+        if debug: print '  ' + main + ': creating times _______'
+        for it in range(Ntvalues):
+            tSvals = ''
+            for tvar in timevars:
+                tSvals = tSvals + datavalues[tvar][it] + ' '
+
+            timeSvals.append(tSvals[0:len(tSvals)-1])
+            if debug: print it, '*' + timeSvals[it] + '*'
+
+        timevals = Stringtimes_CF(timeSvals, description['FMTtime'], referencedate,      \
+          timeunits, debug)
+    else:
+        timevals = Stringtimes_CF(datavalues[description['NAMEtime']],                   \
+          description['FMTtime'], referencedate, timeunits, debug)
+
+    CFtimeRef = referencedate[0:4] +'-'+ referencedate[4:6] +'-'+ referencedate[6:8] +   \
+      ' ' + referencedate[8:10] +':'+ referencedate[10:12] +':'+ referencedate[12:14]
+    tunits = timeunits + ' since ' + CFtimeRef
+
+if objfile.variables.has_key('time'):
+    print warnmsg
+    print '  ' + main + ": variable 'time' already exist !!"
+    print "    renaming it as 'CFtime'"
+    timeCFname = 'CFtime'
+    newdim = objfile.renameDimension('time','CFtime')
+    newvar = objfile.createVariable( timeCFname, 'f8', ('CFtime'))
+    basicvardef(newvar, timeCFname, 'time', tunits )
+else:
+    timeCFname = 'time'
+    newvar = objfile.createVariable( timeCFname, 'f8', ('time'))
+    basicvardef(newvar, timeCFname, 'time', tunits )
+
+set_attribute(newvar, 'calendar', 'standard')
+newvar[:] = timevals
+
+# Global attributes
+##
+for descn in description.keys():
+    if descn[0:3] != 'var' and descn[0:4] != 'NAME' and descn[0:3] != 'FMT':
+        set_attribute(objfile, descn, description[descn])
+
+set_attribute(objfile,'author_nc','Lluis Fita')
+set_attribute(objfile,'institution_nc','Laboratoire de Meteorology Dynamique, ' +    \
+  'LMD-Jussieu, UPMC, Paris')
+set_attribute(objfile,'country_nc','France')
+set_attribute(objfile,'script_nc',main)
+set_attribute(objfile,'version_script',version)
+
+# Adding three variables with the station location, longitude, latitude and height
+adding_station_desc(objfile,stationdesc)
+
+objfile.sync()
+objfile.close()
+
+print main + ": Successfull generation of netcdf observational file '" + ofile + "' !!"
