Context Navigation

← Previous Changeset
Next Changeset →

Changeset 747 in lmdz_wrf

Timestamp:

May 4, 2016, 11:19:50 AM (9 years ago)

Author:

lfita

Message:

Ready version to test and suffer with the use of `generic_tools'.

Now each function (when used) will demand the readiness with gen

File:

: 1 edited

trunk/tools/nc_var_tools.py (modified) (37 diffs)

Legend:

: Unmodified
: Added
: Removed

trunk/tools/nc_var_tools.py

-                      r746
+                      r747
                   ncf = NetCDFFile(ncfile,'a')
     return
+  return
 def chvarname(values, ncfile, varn):
 …
                 self.dimx=self.dims[3]
     return
+        return
 def subyrs(values, ncfile, varn):
 …
     return loadvar
-AQUI AQUI AQUI AQUI AQUI AQUI AQUI
 def fill_ncvariable_lastdims(ncf, varn, prevdims, Nlastdims, fillvals):
     """ Function to fill the last [Nlastdims] dimensions of a variable [varn] from a netCDF object at the other dimensions at [prevdims]
 …
         quit(-1)
+    return
 def maskvar(values, filename, varn):
     """ Function to mask a variable using a mask. It is assumed that mask[...,dimM,dimJ,dimK] and
 …
     print 'masked file "' + ofile + '" generated !!!'
+def lonlat2D(ncfobj, lonn, latn):
+    """ Function to create 2D matricies of longitudes and latitudes from the file
+    ncfobj= netCDF object file
+    lonn= longitude name
+    latn= latitude name
+    """
+    fname='lonlat2D'
+    lonvar = ncfobj.variables[lonn]
+    latvar = ncfobj.variables[latn]
+    loninf = variable_inf(lonvar)
+    latinf = variable_inf(latvar)
+    if loninf.Ndims == 1:
+        lonv = lonvar[:]
+        latv = latvar[:]
+        dx=loninf.dims[0]
+        dy=latinf.dims[0]
+        lonval = np.zeros((dy, dx), dtype=float)
+        latval = np.zeros((dy, dx), dtype=float)
+        for iy in range(dy):
+            lonval[iy,:] = lonv
+        for ix in range(dx):
+            latval[:,ix] = latv
+    elif loninf.Ndims == 2:
+        lonval = lonvar[:]
+        latval = latvar[:]
+    elif loninf.Ndims == 3:
+        lonval = lonvar[0,:,:]
+        latval = latvar[0,:,:]
+    else:
+        print errormsg
+        print '  ' + fname + ' dimension ', lonvar.Ndims, ' of the lon/lat matrices not ready !!!!'
+    return lonval, latval
+def indices_mask(mask):
+    """ Function to retrieve the indices from a 2D mask
+    mask= 2D mask matrix
+    mat=np.array(range(100), dtype=float).reshape((10,10))
+    >>> mask1= mat >= 45.
+    >>> mask2= mat <= 55.
+    >>> mask = mask1 * mask2
+    indices_mask(mask)
+    [[4 5]
+    [4 6]
+    [4 7]
+    [4 8]
+    [4 9]
+    [5 0]
+    [5 1]
+    [5 2]
+    [5 3]
+    [5 4]
+    [5 5]]
+    """
+    dy=mask.shape[0]
+    dx=mask.shape[1]
+    Nptmask = np.sum(mask)
+    maskindices = np.zeros((Nptmask,2), dtype=int)
+    ij=0
+    for j in range(dy):
+        for i in range(dx):
+            if mask[j,i] == True:
+                maskindices[ij,0] = j
+                maskindices[ij,1] = i
+                ij = ij + 1
+    return maskindices
+    return
 def chgtimestep(values, origfile, varN):
 …
     ncofile.close()
+    return
 def uploading_timestep(ncfile, varN, time_step):
     """ Function to upload a given time-step of a variable inside a netCDF
 …
     fileobj.close()
+def typemod(value, typeval):
+    """ Function to give back a value in a given dtype
+    >>> print(typemod(8.2223, 'np.float64'))
+    <type 'numpy.float64'>
+    >>> print(typemod(8.2223, 'tuple'))
+    <type 'tuple'>
+    """
+    fname='typemod'
+    if typeval == 'int':
+        return int(value)
+    elif typeval == 'long':
+        return long(value)
+    elif typeval == 'float':
+        return float(value)
+    elif typeval == 'complex':
+        return complex(value)
+    elif typeval == 'str':
+        return str(value)
+    elif typeval == 'bool':
+        return bool(value)
+    elif typeval == 'list':
+        newval = []
+        newval.append(value)
+        return newval
+    elif typeval == 'dic':
+        newval = {}
+        newval[value] = value
+        return newval
+    elif typeval == 'tuple':
+        newv = []
+        newv.append(value)
+        newval = tuple(newv)
+        return newval
+    elif typeval == 'np.int8':
+        return np.int8(value)
+    elif typeval == 'np.int16':
+        return np.int16(value)
+    elif typeval == 'np.int32':
+        return np.int32(value)
+    elif typeval == 'np.int64':
+        return np.int64(value)
+    elif typeval == 'np.uint8':
+        return np.uint8(value)
+    elif typeval == 'np.uint16':
+        return np.uint16(value)
+    elif typeval == 'np.np.uint32':
+        return np.uint32(value)
+    elif typeval == 'np.uint64':
+        return np.uint64(value)
+    elif typeval == 'np.float':
+        return np.float(value)
+    elif typeval == 'np.float16':
+        return np.float16(value)
+    elif typeval == 'np.float32':
+        return np.float32(value)
+    elif typeval == 'float32':
+        return np.float32(value)
+    elif typeval == 'np.float64':
+        return np.float64(value)
+    elif typeval == 'float64':
+        return np.float64(value)
+    elif typeval == 'np.complex':
+        return np.complex(value)
+    elif typeval == 'np.complex64':
+        return np.complex64(value)
+    elif typeval == 'np.complex128':
+        return np.complex128(value)
+    else:
+        print errormsg
+        print fname + ':  data type "' + typeval + '" is not ready !'
+        print errormsg
+        quit(-1)
+    return
 def addvals(values, filename, varN):
 …
         quit(-1)
+    return
 def TimeInf(filename, tname):
     """ Function to print all the information from the variable time
 …
     fmtprinting_class(timeinf)
     ncfobj.close()
+    return
 def sorttimesmat(filename, tname):
 …
     ncft.close()
+def diffdate_units(diffdate, units):
+    """ Function to transform a difference of dates to a certain units
+    diffdate = difference of dates (as timedelta)
+    units = units to transform: 'weeks', 'days', 'hours', 'minutes', 'seconds', 'miliseconds', 'nanoseconds'
+    >>> diffdate_units(dt.datetime(1976, 2, 17, 8, 32, 5) - dt.datetime(1949, 12, 1, 0, 0, 0), 'seconds')
+    827224325.0
+    """
+    fname = 'diffdate_units'
+    if units == 'weeks':
+        diffunits = diffdate.days/7. + diffdate.seconds/(3600.*24.*7.)
+    elif units == 'days':
+        diffunits = diffdate.days + diffdate.seconds/(3600.*24.)
+    elif units == 'hours':
+        diffunits = diffdate.days*24. + diffdate.seconds/(3600.)
+    elif units == 'minutes':
+        diffunits = diffdate.days*24.*60. + diffdate.seconds/(60.)
+    elif units == 'seconds':
+        diffunits = diffdate.days*24.*3600. + diffdate.seconds
+    elif units == 'miliseconds':
+        diffunits = diffdate.days*24.*3600.*1000. + diffdate.seconds*1000.
+    elif units == 'nanoseconds':
+        diffunits = diffdate.days*24.*3600.*1000000. + diffdate.seconds*1000000.
+    else:
+        print errormsg
+        print fname + ': time units "' + units + '" not ready!!!'
+        print errormsg
+        quit(-1)
+    return diffunits
+def days_month(year,month):
+    """ Function to give the number of days of a month of a given year
+    >>> days_month(1976,2)
+    """
+    import datetime as dt
+    date1=dt.date(year,month,1)
+    date2=dt.date(year,month+1,1)
+    diffdate = date2-date1
+    return diffdate.days
+def mid_period(yr1,mon1,day1,hour1,min1,sec1,yr2,mon2,day2,hour2,min2,sec2):
+    """ Function to give the mid poiint of a period
+    >>> mid_period(1976,2,1,0,0,0,1976,3,1,0,0,0)
+    [1976    2   15   12    0    0]
+    """
+    import datetime as dt
+    date1=dt.datetime(yr1,mon1,day1,hour1,min1,sec1)
+    date2=dt.datetime(yr2,mon2,day2,hour2,min2,sec2)
+    diffdate = date2-date1
+    diffseconds = diffdate.days*24*3600 + diffdate.seconds
+    diff2 = dt.timedelta(seconds=diffseconds/2)
+    datenew = date1 + diff2
+    datenewvals = np.array([datenew.year, datenew.month, datenew.day, datenew.hour,  \
+      datenew.minute, datenew.second])
+    return datenewvals
+def days_year(year):
+    """ Function to give the number of days of a year
+    >>> days_year(1976)
+    """
+    import datetime as dt
+    date1=dt.date(year,1,1)
+    date2=dt.date(year+1,1,1)
+    diffdate = date2-date1
+    return diffdate.days
+def days_period(yeari,yearf):
+    """ Function to give the number of days for a period
+    >>> days_period(1976,1980)
+    """
+    ndays=0
+    for iyr in range(yearf-yeari+1):
+        ndays=ndays+days_year(yeari+iyr)
+    return ndays
+    return
 def check_times_file(values, filename, timen):
 …
                 print matdates[it,:]
+    return
 def writing_str_nc(varo, values, Lchar):
     """ Function to write string values in a netCDF variable as a chain of 1char values
 …
     return strings
-class ysubsetstats(object):
-    """ Class to compute multi-year statistics of a given subset of values
-    values= values to use
-    dates=dates of the values in matrix format ([year], [month], [day], [hour], [minute], [second])
-    sec=section to use for the period='year', 'month', 'day', 'hour', 'minute', 'second'
-    vali=initial value of the period
-    valf=final value of the period
-    Nq= number of quantiles (20 for 5% bins, it is going to produce 21 values)
-    missval= missing value
-      self.Nvalues = Number of non masked values
-      self.min = minimum non masked value
-      self.max = maximum non masked value
-      self.mean = mean non masked value
-      self.mean2 = quandratic mean non masked value
-      self.stdev = standard deviation non masked value
-      self.quantiles = quantiles non masked value
-    """
-    def __init__(self, values, dates, sec, vali, valf, Nq, missVal):
-        import numpy.ma as ma
-        fname = 'ysubsetstats'
-        if values is None:
-            self.Nvalues = None
-            self.min = None
-            self.max = None
-            self.mean = None
-            self.mean2 = None
-            self.stdev = None
-            self.quantiles = None
-        else:
-            Nvalues = len(values)
-            Ndates = len(dates)
-            if Nvalues != Ndates:
-                print errormsg
-                print fname + ': number of values ', Nvalues,' does not coincide with the number of dates ',Ndates, '!!!!!!'
-                print errormsg
-                quit(-1)
-            initialsub = np.zeros(Nvalues, dtype=bool)
-            endsub = initialsub.copy()
-            if missVal > 0.:
-                valuesmask = ma.masked_greater_equal(values, missVal*0.9)
-            else:
-                valuesmask = ma.masked_less_equal(values, missVal*0.9)
-            if sec == 'year':
-                isec = 0
-            elif sec == 'month':
-                isec = 1
-            elif sec == 'day':
-                isec = 2
-            elif sec == 'hour':
-                isec = 3
-            elif sec == 'minute':
-                isec = 4
-            elif sec == 'second':
-                isec = 5
-            else:
-                print errormsg
-                print fname + ': sction "' + sec + '" not ready!!!!'
-                print errormsg
-                quit(-1)
-            if vali < valf:
-                timesmaskperi = ma.masked_less(dates[:,isec], vali)
-                timesmaskperf = ma.masked_greater(dates[:,isec], valf)
-            else:
-                timesmaskperi = ma.masked_less(dates[:,isec], valf)
-                timesmaskperf = ma.masked_greater(dates[:,isec], vali)
-            finalmask = valuesmask.mask+timesmaskperi.mask+timesmaskperf.mask
-            finalvalues = ma.array(values, mask=finalmask)
-            finalvalues2 = finalvalues*finalvalues
-            self.Nvalues = finalvalues.count()
-            self.min = finalvalues.min()
-            self.max = finalvalues.max()
-            self.mean = finalvalues.mean()
-            self.mean2 = finalvalues2.mean()
-            self.stdev = finalvalues.std()
-            self.quantiles = mask_quantiles(finalvalues, Nq)
 def remapnn(values, filename, varname):
 …
     print 'remapnn: new file "' + ofile + '" with remapped variable written !!'
+    return
 #remapnn('met_em.d01.1979-01-01_00:00:00.nc|XLONG_M,XLAT_M|east_weast,south_north|longitude,latitude|longitude,latitude', \
 …
     print 'remapnn: new file "' + ofile + '" with remapped variable written !!'
+    return
 #remapnn('data/R1_CCRC_NARCliM_MOM_1950-2009_pracc.nc:lon,lat:lon,lat', 'data/AWAP_pre_1950-2009_pre_mask.nc', 'pre', 1. , 1.)
 #quit()
 …
     ncobj.close()
+def give_fix1dim_values(mat,iddim,dimval):
+    """ Function to return matrix values when one dimension has been removed at a
+      given value
+    mat = matrix
+    iddim = dimesion to remove (starting from 0 tacking the right most)
+    dimval = value of the dimension
+    >>> matrix = np.array( (range(27)), dtype=np.float).reshape(3,3,3)
+    give_fix1dim_values(matrix,2,2)
+    [[[  0.   1.   2.]
+      [  3.   4.   5.]
+      [  6.   7.   8.]]
+     [[  9.  10.  11.]
+      [ 12.  13.  14.]
+      [ 15.  16.  17.]]]
+    return
+def valmod_dim(values, ncfile, varn):
+    """ Function to modify the value of a variable at a given dimension and value
+    values = dimn,dimval,modins,modval1,[modval2,...]
+      dimn = name of the dimension
+      dimval = value of the dimension at which change the values
+      modins = instruction:
+        'sumc', add [modval1]
+        'subc', substraction [modval1]
+        'mulc', multiply by [modval1]
+        'divc', divide by [modval1]
+        'lowthres': modify all values below [modval1] to [modval2]
+        'upthres': modify all values above [modval1] to [modval2]
+    ncfile = netCDF file name
+    varn = name of the variable
+    >>> valmod_dim('num_metgrid_levels,10,mulc,-1','aqua_met-em_control.nc','PRES')
     """
+    fname='give_fix1dim_values'
+# Not working ma.squeeze
+##    import numpy.ma as ma
+##    print ma.squeeze(mat, axis = (idim,))
+    matshape = mat.shape
+    mattype = mat.dtype
+    matmask = np.ones( tuple(matshape), dtype=bool)
+    fname = 'valmod_dim'
+    if values == '-h':
+        print """ Function to modify the value of a variable at a given dimension and value
+          values = dimn,dimval,modins,modval1,[modval2,...]
+          dimn = name of the dimension
+          dimval = value of the dimension at which change the values
+          modins = instruction:
+            'sumc', add [modval1]
+            'subc', substraction [modval1]
+            'mulc', multiply by [modval1]
+            'divc', divide by [modval1]
+            'lowthres': modify all values below [modval1] to [modval2]
+            'upthres': modify all values above [modval1] to [modval2]
+          ncfile = netCDF file name
+          varn = name of the variable
+          >>> valmod_dim('num_metgrid_levels,10,mulc,-1','aqua_met-em_control.nc','PRES')
+          """
+        quit()
+##  mat = np.array( (range(27)), dtype=np.float).reshape(3,3,3)
+    vals = values.split(',')
+    Lvals = len(vals)
+    dimn = vals[0]
+    dimval = np.int(vals[1])
+    modins = vals[2]
+    modval = np.float(vals[3])
+    newvalues = ','.join(vals[2:Lvals])
+    if not os.path.isfile(ncfile):
+      print errormsg
+      print '   valmod: File "' + ncfile + '" does not exist !!'
+      print errormsg
+      quit(-1)
+    ncf = NetCDFFile(ncfile,'a')
+    if ncf.dimensions.has_key('plev'):
+      # removing pressure level from variable name
+      varn = re.sub("\d+", "", varn)
+    if not ncf.dimensions.has_key(dimn):
+      print errormsg
+      print '   ' + fname + ': File "' + ncfile + '" does not have dimension "' +    \
+        dimn + '" !!!!'
+      ncf.close()
+      quit(-1)
+    if not ncf.variables.has_key(varn):
+      print errormsg
+      print '   ' + fname + ': File "' + ncfile + '" does not have variable "' +     \
+        varn + '" !!!!'
+      ncf.close()
+      quit(-1)
+    varobj = ncf.variables[varn]
+    varinf = variable_inf(varobj)
+    varvals = varobj[:]
+# Dimension id
+    iddim=varinf.Ndims - 1
+    for idimn in varinf.dimns:
+#        print iddim, idimn
+        if dimn == idimn:
+            break
+        iddim = iddim - 1
+    matshape = varvals.shape
+    mattype = varvals.dtype
+    matmask = np.zeros( tuple(matshape), dtype=bool)
     zerodims=np.ones( (6), dtype=int)
     Ndims = len(matshape)
     if Ndims > 6:
         print errmsg
 …
     zerodims[5-Ndims+1:6] = matshape
-    newdims = list(matshape)
-    newdims[Ndims-1-iddim] = matshape[iddim]-1
-    mask = np.ones(tuple(zerodims), dtype=bool)
-    vals = np.zeros(tuple(zerodims), dtype=mattype)
-    vals[0:zerodims[0],0:zerodims[1],0:zerodims[2],0:zerodims[3],0:zerodims[4],
-:zerodims[5]] = mat
-    for i5 in range(zerodims[5]):
-        for i4 in range(zerodims[4]):
-            for i3 in range(zerodims[3]):
-                for i2 in range(zerodims[2]):
-                    for i1 in range(zerodims[1]):
-                        for i0 in range(zerodims[0]):
-#                            print 'iddim:',iddim,'dimval:',dimval,'-->',i0,i1,i2,i3,i4,i5
-                            if iddim == 5 and dimval == i0:
-                                mask[i0,i1,i2,i3,i4,i5] = False
-                            elif iddim == 4 and dimval == i1:
-                                mask[i0,i1,i2,i3,i4,i5] = False
-                            elif iddim == 3 and dimval == i2:
-                                mask[i0,i1,i2,i3,i4,i5] = False
-                            elif iddim == 2 and dimval == i3:
-                                mask[i0,i1,i2,i3,i4,i5] = False
-                            elif iddim == 1 and dimval == i4:
-                                mask[i0,i1,i2,i3,i4,i5] = False
-                            elif iddim == 0 and dimval == i5:
-                                mask[i0,i1,i2,i3,i4,i5] = False
-    newmat = vals[mask].reshape(tuple(newdims))
-#    print 'newmat dim:',iddim,'val:',dimval,':',newdims,'______'
-#    print newmat
-    return newmat
-def portion_fix1dim_values(mat,iddim,dimval):
-    """ Function to return that portion of the matrix values when one dimension has
-      been removed at a given value
-    mat = matrix
-    iddim = dimesion to remove (starting from 0 tacking the right most)
-    dimval = value of the dimension
-    >>> matrix = np.array( (range(27)), dtype=np.float).reshape(3,3,3)
-    portion_fix1dim_values(matrix,2,2)
-    [[ 18.  19.  20.]
-     [ 21.  22.  23.]
-     [ 24.  25.  26.]]
-    """
-    fname='portion_fix1dim_values'
-# Not working ma.squeeze
-##    import numpy.ma as ma
-##    print ma.squeeze(mat, axis = (idim,))
-    matshape = mat.shape
-    mattype = mat.dtype
-    matmask = np.zeros( tuple(matshape), dtype=bool)
-    zerodims=np.ones( (6), dtype=int)
-    Ndims = len(matshape)
-    if Ndims > 6:
-        print errmsg
-        print '  ' + fname + ' number of dimensions: ',Ndims,' not ready!!!!!'
-        quit(-1)
-    zerodims[5-Ndims+1:6] = matshape
     newdims = []
     for idim in range(Ndims):
+#        print idim,'matshape[idim]:',matshape[Ndims-1-idim],'iddim:',iddim
         if idim != iddim:
             newdims.append(matshape[idim])
+            newdims.append(matshape[Ndims-1-idim])
     newdims.reverse()
 …
     mask = np.zeros(tuple(zerodims), dtype=bool)
     vals = np.zeros(tuple(zerodims), dtype=mattype)
+    vals[0:zerodims[0],0:zerodims[1],0:zerodims[2],0:zerodims[3],0:zerodims[4],
+:zerodims[5]] = mat
+    newvals = np.zeros(tuple(zerodims), dtype=mattype)
+    vals[0:zerodims[0],0:zerodims[1],0:zerodims[2],0:zerodims[3],0:zerodims[4],      \
+:zerodims[5]] = varvals
+#    print 'vals shape:',vals.shape,':',newdims
     for i5 in range(zerodims[5]):
 …
 #    print 'newmat dim:',iddim,'val:',dimval,':',newdims,'______'
     newmat = vals[mask].reshape(tuple(newdims))
-#    print newmat
-    return newmat
-def valmod_dim(values, ncfile, varn):
-    """ Function to modify the value of a variable at a given dimension and value
-    values = dimn,dimval,modins,modval1,[modval2,...]
-      dimn = name of the dimension
-      dimval = value of the dimension at which change the values
-      modins = instruction:
-        'sumc', add [modval1]
-        'subc', substraction [modval1]
-        'mulc', multiply by [modval1]
-        'divc', divide by [modval1]
-        'lowthres': modify all values below [modval1] to [modval2]
-        'upthres': modify all values above [modval1] to [modval2]
-    ncfile = netCDF file name
-    varn = name of the variable
-    >>> valmod_dim('num_metgrid_levels,10,mulc,-1','aqua_met-em_control.nc','PRES')
-    """
-    fname = 'valmod_dim'
-    if values == '-h':
-        print """ Function to modify the value of a variable at a given dimension and value
-          values = dimn,dimval,modins,modval1,[modval2,...]
-          dimn = name of the dimension
-          dimval = value of the dimension at which change the values
-          modins = instruction:
-            'sumc', add [modval1]
-            'subc', substraction [modval1]
-            'mulc', multiply by [modval1]
-            'divc', divide by [modval1]
-            'lowthres': modify all values below [modval1] to [modval2]
-            'upthres': modify all values above [modval1] to [modval2]
-          ncfile = netCDF file name
-          varn = name of the variable
-          >>> valmod_dim('num_metgrid_levels,10,mulc,-1','aqua_met-em_control.nc','PRES')
-          """
-        quit()
-##  mat = np.array( (range(27)), dtype=np.float).reshape(3,3,3)
-    vals = values.split(',')
-    Lvals = len(vals)
-    dimn = vals[0]
-    dimval = np.int(vals[1])
-    modins = vals[2]
-    modval = np.float(vals[3])
-    newvalues = ','.join(vals[2:Lvals])
-    if not os.path.isfile(ncfile):
-      print errormsg
-      print '   valmod: File "' + ncfile + '" does not exist !!'
-      print errormsg
-      quit(-1)
-    ncf = NetCDFFile(ncfile,'a')
-    if ncf.dimensions.has_key('plev'):
-      # removing pressure level from variable name
-      varn = re.sub("\d+", "", varn)
-    if not ncf.dimensions.has_key(dimn):
-      print errormsg
-      print '   ' + fname + ': File "' + ncfile + '" does not have dimension "' +    \
-        dimn + '" !!!!'
-      ncf.close()
-      quit(-1)
-    if not ncf.variables.has_key(varn):
-      print errormsg
-      print '   ' + fname + ': File "' + ncfile + '" does not have variable "' +     \
-        varn + '" !!!!'
-      ncf.close()
-      quit(-1)
-    varobj = ncf.variables[varn]
-    varinf = variable_inf(varobj)
-    varvals = varobj[:]
-# Dimension id
-    iddim=varinf.Ndims - 1
-    for idimn in varinf.dimns:
-#        print iddim, idimn
-        if dimn == idimn:
-            break
-        iddim = iddim - 1
-    matshape = varvals.shape
-    mattype = varvals.dtype
-    matmask = np.zeros( tuple(matshape), dtype=bool)
-    zerodims=np.ones( (6), dtype=int)
-    Ndims = len(matshape)
-    if Ndims > 6:
-        print errmsg
-        print '  ' + fname + ' number of dimensions: ',Ndims,' not ready!!!!!'
-        quit(-1)
-    zerodims[5-Ndims+1:6] = matshape
-    newdims = []
-    for idim in range(Ndims):
-#        print idim,'matshape[idim]:',matshape[Ndims-1-idim],'iddim:',iddim
-        if idim != iddim:
-            newdims.append(matshape[Ndims-1-idim])
-    newdims.reverse()
-    mask = np.zeros(tuple(zerodims), dtype=bool)
-    vals = np.zeros(tuple(zerodims), dtype=mattype)
-    newvals = np.zeros(tuple(zerodims), dtype=mattype)
-    vals[0:zerodims[0],0:zerodims[1],0:zerodims[2],0:zerodims[3],0:zerodims[4],      \
-:zerodims[5]] = varvals
-#    print 'vals shape:',vals.shape,':',newdims
-    for i5 in range(zerodims[5]):
-        for i4 in range(zerodims[4]):
-            for i3 in range(zerodims[3]):
-                for i2 in range(zerodims[2]):
-                    for i1 in range(zerodims[1]):
-                        for i0 in range(zerodims[0]):
-#                            print 'iddim:',iddim,'dimval:',dimval,'-->',i0,i1,i2,i3,i4,i5
-                            if iddim == 5 and dimval == i0:
-                                mask[i0,i1,i2,i3,i4,i5] = True
-                            elif iddim == 4 and dimval == i1:
-                                mask[i0,i1,i2,i3,i4,i5] = True
-                            elif iddim == 3 and dimval == i2:
-                                mask[i0,i1,i2,i3,i4,i5] = True
-                            elif iddim == 2 and dimval == i3:
-                                mask[i0,i1,i2,i3,i4,i5] = True
-                            elif iddim == 1 and dimval == i4:
-                                mask[i0,i1,i2,i3,i4,i5] = True
-                            elif iddim == 0 and dimval == i5:
-                                mask[i0,i1,i2,i3,i4,i5] = True
-#    print 'newmat dim:',iddim,'val:',dimval,':',newdims,'______'
-    newmat = vals[mask].reshape(tuple(newdims))
     varmod = np.zeros(tuple(newdims), dtype=mattype)
 …
     ncf.sync()
     ncf.close()
+    return
 def checkNaNs(values, filen):
 …
     ncobj.close()
+    return
 ##checkNaNs(10.e10, '/d4/lflmd/etudes/WRF_LMDZ/WaquaL/WRF/control/wrfout_d01_1979-12-01_00:00:00')
 ##checkNaNs(10.e10, '/d4/lflmd/etudes/FF/tests/em_quarter_ss/run/wrfout_d01_0001-01-01_00:00:00')
 …
     ncobj.close()
+    return
 ##checkAllValues(0., '/d4/lflmd/etudes/WRF_LMDZ/WaquaL/WRF/control/wrfout_d01_1979-12-01_00:00:00')
 …
     return
 def sellonlatboxold(values, ncfile, varn):
 …
 #compute_opervaralltime('add|time', 'obs/PRCP.nc', 'product')
-def retype(val, vtype):
-    """ Function to transform a value to a given type
-    retype(val, vtype)
-      [val]= value
-      [vtype]= type to transform
-    >>> retype(0, type(True))
-    False
-    """
-    fname = 'retype'
-    if val == 'h':
-        print fname + '_____________________________________________________________'
-        print retype.__doc__
-        quit()
-    if vtype == type(int(1)):
-        newval = int(val)
-    elif vtype == type(float(1)):
-        newval = float(val)
-#    elif vtype == type(float32(1)):
-#        newval = float32(val)
-    elif vtype == type(np.int(1)):
-        typeinf = np.iinfo(np.int)
-        if (abs(np.float64(val)) > typeinf.max):
-            print warnmsg
-            print '  ' + fname + ': value to transform ', val,' overpasses type:',   \
-               vtype,' limits!!'
-            print '  Changing value to kind largest allowed value:', typeinf.max
-            newval = abs(np.float64(val))/np.float64(val) * np.int(typeinf.max)
-        else:
-            newval = np.int(val)
-    elif vtype == type(np.int16(1)):
-        typeinf = np.iinfo(np.int16)
-        if (abs(np.float64(val)) > typeinf.max):
-            print warnmsg
-            print '  ' + fname + ': value to transform ', val,' overpasses type:',   \
-               vtype,' limits!!'
-            print '  Changing value to kind largest allowed value:', typeinf.max
-            newval = abs(np.float64(val))/np.float64(val) * np.int16(typeinf.max)
-        else:
-            newval = np.int16(val)
-    elif vtype == type(np.int32(1)):
-        typeinf = np.iinfo(np.int32)
-        if (np.abs(np.float64(val)) > typeinf.max):
-            print warnmsg
-            print '  ' + fname + ': value to transform ', val,' overpasses type:',   \
-               vtype,' limits!!'
-            print '  Changing value to kind largest allowed value:', typeinf.max
-            newval = abs(np.float64(val))/np.float64(val) * np.int32(typeinf.max)
-        else:
-            newval = np.int32(np.float64(val))
-    elif vtype == type(np.int64(1)):
-        newval = np.int64(val)
-    elif vtype == type(np.float(1)):
-        newval = np.float(val)
-#    elif vtype == type(np.float16(1)):
-#        newval = np.float16(val)
-    elif vtype == type(np.float32(1)):
-        newval = np.float32(val)
-    elif vtype == type(np.float64(1)):
-        newval = np.float64(val)
-    elif vtype == type(True):
-        if val == 0:
-            newval = False
-        else:
-            newval = True
-    elif vtype == '|S1':
-        newval = str(val)
-    else:
-        print errormsg
-        print '  ' + fname + ': variable type "', vtype, '" not ready!!!'
-        quit(-1)
-    return newval
 def setvar_asciivalues(values, ncfile, varname):
     """ Function to set de values of a variable with an ASCII file (common
 …
     objnc.close()
+    return
 #setvar_asciivalues('obstimes.inf', 'obs/re_project_accumulated-1h_PRCP.nc', 'time')
-class stats_space2D(object):
-    """spatial statistics for a 2D file:
-      vals= variable ro check (assuming var[t,dy,dx])
-        self.minv[t]: spatial minimum at each time-step
-        self.maxv[t]: spatial maximum at each time-step
-        self.meanv[t]: spatial mean at each time-step
-        self.mean2v[t]: spatial quadratic mean at each time-step
-        self.stdv[t]: spatial standard deviation at each time-step
-        self.anomv[t]: spatial anomaly from the whole mean at each time-step
-    """
-    def __init__(self, vals):
-        from scipy import stats as sts
-        fname = 'stats_space2D'
-        if vals1 == 'h':
-            print fname + '_____________________________________________________________'
-            print stats_space2Dvars.__doc__
-            quit()
-        if vals1 is None:
-            self.minv = None
-            self.maxv = None
-            self.meanv = None
-            self.mean2v = None
-            self.stdv = None
-        else:
-            dimt = vals.shape[0]
-            stats=np.zeros((dimt,5), dtype=np.float)
-            absmean = np.mean(vals)
-            for it in range(dimt):
-                stats[it,0]=np.min(vals[it,:,:])
-                stats[it,1]=np.max(vals[it,:,:])
-                stats[it,2]=np.mean(vals[it,:,:])
-                stats[it,3]=np.mean(vals[it,:,:]*vals[it,:,:])
-                stats[it,4]=absmean - stats[it,2]
-            stats = np.where(stats > 0.1*fillValue, fillValue, stats)
-            stats = np.where(stats is nan, fillValue, stats)
-            stats = np.where(stats is np.inf, fillValue, stats)
-            stats = np.where(stats is None, fillValue, stats)
-            self.minv=stats[:,0]
-            self.maxv=stats[:,1]
-            self.meanv=stats[:,2]
-            self.mean2v=stats[:,3]
-            self.stdv=np.sqrt(stats[:,3]-stats[:,2]*stats[:,2])
-            self.anomv=stats[:,4]
-class stats_space2Dvars(object):
-    """spatial statistics beween 2 2D files:
-    valsA = variable 1 (assuming var[t,dy,dx])
-    valsB = variable 2 (assuming var[t,dy,dx])
-    self.min[t,var], self.max[t,var], self.mean[t,var], self.mean2[t,var],
-      self.std[t,var]
-      [var] = var1+var2[v1Av2], var1-var2[v1Sv2], var1/var2[v1Dv2], var1*var2[v1Pv2]
-    self.mae=mean[t,abs(var1-var2)]
-    self.correlation=correlation[var1,var2] (and p-value)
-    self.bias=[t,mean(var1)-mean(var2)]
-    self.meancorrelation=correlation[mean(var1),mean(var2)] (and p-value)
-    """
-    def __init__(self, vals1, vals2):
-        from scipy import stats as sts
-        fname = 'stats_space2Dvars'
-        if vals1 == 'h':
-            print fname + '_____________________________________________________________'
-            print stats_space2Dvars.__doc__
-            quit()
-        if vals1 is None:
-            self.minv1Av2 = None
-            self.maxv1Av2 = None
-            self.meanv1Av2 = None
-            self.mean2v1Av2 = None
-            self.stdv1Av2 = None
-            self.minv1Sv2 = None
-            self.maxv1Sv2 = None
-            self.meanv1Sv2 = None
-            self.mean2v1Sv2 = None
-            self.stdv1Sv2 = None
-            self.minv1Dv2 = None
-            self.maxv1Dv2 = None
-            self.meanv1Dv2 = None
-            self.mean2v1Dv2 = None
-            self.stdv1Dv2 = None
-            self.minv1Pv2 = None
-            self.maxv1Pv2 = None
-            self.meanv1Pv2 = None
-            self.mean2v1Pv2 = None
-            self.stdv1Pv2 = None
-            self.mae = None
-            self.corr = None
-        else:
-            dimt = vals1.shape[0]
-            stats=np.zeros((dimt,26), dtype=np.float)
-            meanvals1 = np.zeros((dimt), dtype=np.float)
-            meanvals2 = np.zeros((dimt), dtype=np.float)
-            for it in range(dimt):
-                v1 = vals1[it,:,:].flatten()
-                v2 = vals2[it,:,:].flatten()
-# add
-                vs = v1 + v2
-                stats[it,0] = np.min(vs)
-                stats[it,1] = np.max(vs)
-                stats[it,2] = np.mean(vs)
-                stats[it,3] = np.mean(vs*vs)
-                stats[it,4] = np.sqrt(stats[it,3] - stats[it,2]*stats[it,2])
-# sub
-                stats[it,20] = np.mean(np.abs(v1-v2))
-                vs = v1 - v2
-                stats[it,5] = np.min(vs)
-                stats[it,6] = np.max(vs)
-                stats[it,7] = np.mean(vs)
-                stats[it,8] = np.mean(vs*vs)
-                stats[it,9] = np.sqrt(stats[it,8] - stats[it,7]*stats[it,7])
-# mul
-                vs = v1 * v2
-                stats[it,10] = np.min(vs)
-                stats[it,11] = np.max(vs)
-                stats[it,12] = np.mean(vs)
-                stats[it,13] = np.mean(vs*vs)
-                stats[it,14] = np.sqrt(stats[it,13] - stats[it,12]*stats[it,12])
-# div
-                vs = v1 / v2
-                stats[it,15] = np.min(vs)
-                stats[it,16] = np.max(vs)
-                stats[it,17] = np.mean(vs)
-                stats[it,18] = np.mean(vs*vs)
-                stats[it,19] = np.sqrt(stats[it,18] - stats[it,17]*stats[it,17])
-# corr
-                stats[it,21], stats[it,22] = mask_pearsonr(v1, v2)
-# Mean values
-                meanvals1[it] = np.mean(v1)
-                meanvals2[it] = np.mean(v2)
-                stats[it,23] =  meanvals1[it] - meanvals2[it]
-            stats = np.where(stats > 0.1*fillValue, fillValue, stats)
-            stats = np.where(stats is np.nan, fillValue, stats)
-            stats = np.where(stats is np.inf, fillValue, stats)
-            stats = np.where(stats is None, fillValue, stats)
-            self.minv1Av2 = stats[:,0]
-            self.maxv1Av2 = stats[:,1]
-            self.meanv1Av2 = stats[:,2]
-            self.meanv21Av2 = stats[:,3]
-            self.stdv1Av2 = stats[:,4]
-            self.minv1Sv2 = stats[:,5]
-            self.maxv1Sv2 = stats[:,6]
-            self.meanv1Sv2 = stats[:,7]
-            self.meanv21Sv2 = stats[:,8]
-            self.stdv1Sv2 = stats[:,9]
-            self.minv1Pv2 = stats[:,10]
-            self.maxv1Pv2 = stats[:,11]
-            self.meanv1Pv2 = stats[:,12]
-            self.meanv21Pv2 = stats[:,13]
-            self.stdv1Pv2 = stats[:,14]
-            self.minv1Dv2 = stats[:,15]
-            self.maxv1Dv2 = stats[:,16]
-            self.meanv1Dv2 = stats[:,17]
-            self.meanv21Dv2 = stats[:,18]
-            self.stdv1Dv2 = stats[:,19]
-            self.mae = stats[:,20]
-            self.corr = stats[:,21]
-            self.p_value = stats[:,22]
-            self.bias = stats[:,23]
-            self.meancorr, self.meanp_value = sts.pearsonr(meanvals1, meanvals2)
-class stats_time2D(object):
-    """temporal statistics for a 2D file:
-      vals= variable ro check (assuming var[t,dy,dx])
-        self.minv[t]: temporal minimum at each grid point
-        self.maxv[t]: temporal maximum at each grid point
-        self.meanv[t]: temporal mean at each grid point
-        self.mean2v[t]: temporal quadratic mean at each grid point
-        self.stdv[t]: temporal standard deviation at each grid point
-        self.anomv[t]: temporal anomaly from the whole mean at each grid point
-    """
-    def __init__(self, vals):
-        from scipy import stats as sts
-        fname = 'stats_time2D'
-        if vals1 == 'h':
-            print fname + '_____________________________________________________________'
-            print stats_time2Dvars.__doc__
-            quit()
-        if vals1 is None:
-            self.minv = None
-            self.maxv = None
-            self.meanv = None
-            self.mean2v = None
-            self.stdv = None
-        else:
-            dimx = vals.shape[2]
-            dimy = vals.shape[1]
-            stats=np.zeros((dimy,dimx,5), dtype=np.float)
-            absmean = np.mean(vals,axis=0)
-            stats[:,:,0]=np.min(vals, axis=0)
-            stats[:,:,1]=np.max(vals, axis=0)
-            stats[:,:,2]=np.mean(vals, axis=0)
-            stats[:,:,3]=np.mean(vals*vals, axis=0)
-            stats[:,:,4]=absmean - stats[:,:,2]
-            stats = np.where(stats > 0.1*fillValue, fillValue, stats)
-            stats = np.where(stats is np.nan, fillValue, stats)
-            stats = np.where(stats is np.inf, fillValue, stats)
-            stats = np.where(stats is None, fillValue, stats)
-            self.minv=stats[:,:,0]
-            self.maxv=stats[:,:,1]
-            self.meanv=stats[:,:,2]
-            self.mean2v=stats[:,:,3]
-            self.stdv=np.sqrt(stats[:,:,3]-stats[:,:,2]*stats[:,:,2])
-            self.anomv=stats[:,:,4]
-class stats_time2Dvars(object):
-    """temporal statistics beween 2 2D files:
-    valsA = variable 1 (assuming var[t,dy,dx])
-    valsB = variable 2 (assuming var[t,dy,dx])
-    self.min[dy,dx,var], self.max[dy,dx,var], self.mean[dy,dx,var],
-      self.mean2[dy,dx,var], self.std[dy,dx,var]
-      [var] = var1+var2[v1Av2], var1-var2[v1Sv2], var1/var2[v1Dv2], var1*var2[v1Pv2]
-    self.mae=mean[dy,dx,abs(var1-var2)]
-    self.correlation=correlation[var1,var2] (and p-value)
-    self.bias=[dy,dx,mean(var1)-mean(var2)]
-    self.meancorrelation=correlation[mean(var1),mean(var2)] (and p-value)
-    """
-    def __init__(self, vals1, vals2):
-        from scipy import stats as sts
-        fname = 'stats_time2Dvars'
-        if vals1 == 'h':
-            print fname + '_____________________________________________________________'
-            print stats_time2Dvars.__doc__
-            quit()
-        if vals1 is None:
-            self.minv1Av2 = None
-            self.maxv1Av2 = None
-            self.meanv1Av2 = None
-            self.mean2v1Av2 = None
-            self.stdv1Av2 = None
-            self.minv1Sv2 = None
-            self.maxv1Sv2 = None
-            self.meanv1Sv2 = None
-            self.mean2v1Sv2 = None
-            self.stdv1Sv2 = None
-            self.minv1Dv2 = None
-            self.maxv1Dv2 = None
-            self.meanv1Dv2 = None
-            self.mean2v1Dv2 = None
-            self.stdv1Dv2 = None
-            self.minv1Pv2 = None
-            self.maxv1Pv2 = None
-            self.meanv1Pv2 = None
-            self.mean2v1Pv2 = None
-            self.stdv1Pv2 = None
-            self.mae = None
-            self.corr = None
-        else:
-            dimx = vals1.shape[1]
-            dimy = vals1.shape[1]
-            stats=np.zeros((dimy,dimx,24), dtype=np.float)
-            meanvals1 = np.zeros((dimy,dimx), dtype=np.float)
-            meanvals2 = np.zeros((dimy,dimx), dtype=np.float)
-# add
-            vs = vals1 + vals2
-            stats[:,:,0] = np.min(vs,axis=0)
-            stats[:,:,1] = np.max(vs,axis=0)
-            stats[:,:,2] = np.mean(vs,axis=0)
-            stats[:,:,3] = np.mean(vs*vs,axis=0)
-            stats[:,:,4] = np.sqrt(stats[:,:,3] - stats[:,:,2]*stats[:,:,2])
-# sub
-            stats[:,:,20] = np.mean(np.abs(vals1-vals2), axis=0)
-            vs = vals1 - vals2
-            stats[:,:,5] = np.min(vs,axis=0)
-            stats[:,:,6] = np.max(vs,axis=0)
-            stats[:,:,7] = np.mean(vs,axis=0)
-            stats[:,:,8] = np.mean(vs*vs,axis=0)
-            stats[:,:,9] = np.sqrt(stats[:,:,8] - stats[:,:,7]*stats[:,:,7])
-# mul
-            vs = vals1 * vals2
-            stats[:,:,10] = np.min(vs,axis=0)
-            stats[:,:,11] = np.max(vs,axis=0)
-            stats[:,:,12] = np.mean(vs,axis=0)
-            stats[:,:,13] = np.mean(vs*vs,axis=0)
-            stats[:,:,14] = np.sqrt(stats[:,:,13] - stats[:,:,12]*stats[:,:,12])
-# div
-            vs = vals1 / vals2
-            stats[:,:,15] = np.min(vs,axis=0)
-            stats[:,:,16] = np.max(vs,axis=0)
-            stats[:,:,17] = np.mean(vs,axis=0)
-            stats[:,:,18] = np.mean(vs*vs,axis=0)
-            stats[:,:,19] = np.sqrt(stats[:,:,18] - stats[:,:,17]*stats[:,:,17])
-# Mean values
-            meanvals1[:,:] = np.mean(vals1,axis=0)
-            meanvals2[:,:] = np.mean(vals2,axis=0)
-            stats[:,:,23] =  meanvals1[:,:] - meanvals2[:,:]
-# corr
-            self.meancorr, self.meanp_value = sts.pearsonr(meanvals1.flatten(), meanvals2.flatten())
-            for j in range(dimy):
-                for i in range(dimx):
-                    stats[j,i,21], stats[j,i,22] = sts.pearsonr(vals1[:,j,i], vals2[:,j,i])
-            stats = np.where(stats > 0.1*fillValue, fillValue, stats)
-            stats = np.where(stats is np.nan, fillValue, stats)
-            stats = np.where(stats is np.inf, fillValue, stats)
-            stats = np.where(stats is None, fillValue, stats)
-            self.minv1Av2 = stats[:,:,0]
-            self.maxv1Av2 = stats[:,:,1]
-            self.meanv1Av2 = stats[:,:,2]
-            self.meanv21Av2 = stats[:,:,3]
-            self.stdv1Av2 = stats[:,:,4]
-            self.minv1Sv2 = stats[:,:,5]
-            self.maxv1Sv2 = stats[:,:,6]
-            self.meanv1Sv2 = stats[:,:,7]
-            self.meanv21Sv2 = stats[:,:,8]
-            self.stdv1Sv2 = stats[:,:,9]
-            self.minv1Pv2 = stats[:,:,10]
-            self.maxv1Pv2 = stats[:,:,11]
-            self.meanv1Pv2 = stats[:,:,12]
-            self.meanv21Pv2 = stats[:,:,13]
-            self.stdv1Pv2 = stats[:,:,14]
-            self.minv1Dv2 = stats[:,:,15]
-            self.maxv1Dv2 = stats[:,:,16]
-            self.meanv1Dv2 = stats[:,:,17]
-            self.meanv21Dv2 = stats[:,:,18]
-            self.stdv1Dv2 = stats[:,:,19]
-            self.mae = stats[:,:,20]
-            self.corr = stats[:,:,21]
-            self.p_value = stats[:,:,22]
-            self.bias = stats[:,:,23]
-#vals1 = np.arange(27).reshape(3,3,3)*1.
-#vals2 = np.arange(1,28).reshape(3,3,3)*1.
-#printing_class(stats_time2Dvars(vals1,vals2))
 def statcompare_files(values):
 …
 #statcompare_files('control/sellonlatbox_wss_17-20.nc:wss,obs/re_project_Vu_17-20.nc:wss')
 def sellonlatlevbox(values, ncfile, varn):
 …
 #sellonlatlevbox('XLONG,XLAT,ZNU,-12.4,25.35,32.4,52.65,0,0', 'control/wrfout_d01_2001-11-09_00:00:00', 'P')
-def file_nlines(filen,char):
-    """ Function to provide the number of lines of a file
-    filen= name of the file
-    char= character as no line
-    >>> file_nlines('trajectory.dat','#')
-    """
-    fname = 'file_nlines'
-    if not os.path.isfile(filen):
-        print errormsg
-        print '  ' + fname + ' file: "' + filen + '" does not exist !!'
-        quit(-1)
-    fo = open(filen,'r')
-    nlines=0
-    for line in fo:
-        if line[0:1] != char: nlines = nlines + 1
-    fo.close()
-    return nlines
-def datetimeStr_conversion(StringDT,typeSi,typeSo):
-    """ Function to transform a string date to an another date object
-    StringDT= string with the date and time
-    typeSi= type of datetime string input
-    typeSo= type of datetime string output
-      [typeSi/o]
-        'cfTime': [time],[units]; ]time in CF-convention format [units] = [tunits] since [refdate]
-        'matYmdHMS': numerical vector with [[YYYY], [MM], [DD], [HH], [MI], [SS]]
-        'YmdHMS': [YYYY][MM][DD][HH][MI][SS] format
-        'Y-m-d_H:M:S': [YYYY]-[MM]-[DD]_[HH]:[MI]:[SS] format
-        'Y-m-d H:M:S': [YYYY]-[MM]-[DD] [HH]:[MI]:[SS] format
-        'Y/m/d H-M-S': [YYYY]/[MM]/[DD] [HH]-[MI]-[SS] format
-        'WRFdatetime': [Y], [Y], [Y], [Y], '-', [M], [M], '-', [D], [D], '_', [H],
-          [H], ':', [M], [M], ':', [S], [S]
-    >>> datetimeStr_conversion('1976-02-17_08:32:05','Y-m-d_H:M:S','matYmdHMS')
-    [1976    2   17    8   32    5]
-    >>> datetimeStr_conversion(str(137880)+',minutes since 1979-12-01_00:00:00','cfTime','Y/m/d H-M-S')
-/03/05 18-00-00
-    """
-    import datetime as dt
-    fname = 'datetimeStr_conversion'
-    if StringDT[0:1] == 'h':
-        print fname + '_____________________________________________________________'
-        print datetimeStr_conversion.__doc__
-        quit()
-    if typeSi == 'cfTime':
-        timeval = np.float(StringDT.split(',')[0])
-        tunits = StringDT.split(',')[1].split(' ')[0]
-        Srefdate = StringDT.split(',')[1].split(' ')[2]
-# Does reference date contain a time value [YYYY]-[MM]-[DD] [HH]:[MI]:[SS]
-##
-        yrref=Srefdate[0:4]
-        monref=Srefdate[5:7]
-        dayref=Srefdate[8:10]
-        trefT = Srefdate.find(':')
-        if not trefT == -1:
-#            print '  ' + fname + ': refdate with time!'
-            horref=Srefdate[11:13]
-            minref=Srefdate[14:16]
-            secref=Srefdate[17:19]
-            refdate = datetimeStr_datetime( yrref + '-' + monref + '-' + dayref +    \
-              '_' + horref + ':' + minref + ':' + secref)
-        else:
-            refdate = datetimeStr_datetime( yrref + '-' + monref + '-' + dayref +    \
-              + '_00:00:00')
-        if tunits == 'weeks':
-            newdate = refdate + dt.timedelta(weeks=float(timeval))
-        elif tunits == 'days':
-            newdate = refdate + dt.timedelta(days=float(timeval))
-        elif tunits == 'hours':
-            newdate = refdate + dt.timedelta(hours=float(timeval))
-        elif tunits == 'minutes':
-            newdate = refdate + dt.timedelta(minutes=float(timeval))
-        elif tunits == 'seconds':
-            newdate = refdate + dt.timedelta(seconds=float(timeval))
-        elif tunits == 'milliseconds':
-            newdate = refdate + dt.timedelta(milliseconds=float(timeval))
-        else:
-              print errormsg
-              print '    timeref_datetime: time units "' + tunits + '" not ready!!!!'
-              quit(-1)
-        yr = newdate.year
-        mo = newdate.month
-        da = newdate.day
-        ho = newdate.hour
-        mi = newdate.minute
-        se = newdate.second
-    elif typeSi == 'matYmdHMS':
-        yr = StringDT[0]
-        mo = StringDT[1]
-        da = StringDT[2]
-        ho = StringDT[3]
-        mi = StringDT[4]
-        se = StringDT[5]
-    elif typeSi == 'YmdHMS':
-        yr = int(StringDT[0:4])
-        mo = int(StringDT[4:6])
-        da = int(StringDT[6:8])
-        ho = int(StringDT[8:10])
-        mi = int(StringDT[10:12])
-        se = int(StringDT[12:14])
-    elif typeSi == 'Y-m-d_H:M:S':
-        dateDT = StringDT.split('_')
-        dateD = dateDT[0].split('-')
-        timeT = dateDT[1].split(':')
-        yr = int(dateD[0])
-        mo = int(dateD[1])
-        da = int(dateD[2])
-        ho = int(timeT[0])
-        mi = int(timeT[1])
-        se = int(timeT[2])
-    elif typeSi == 'Y-m-d H:M:S':
-        dateDT = StringDT.split(' ')
-        dateD = dateDT[0].split('-')
-        timeT = dateDT[1].split(':')
-        yr = int(dateD[0])
-        mo = int(dateD[1])
-        da = int(dateD[2])
-        ho = int(timeT[0])
-        mi = int(timeT[1])
-        se = int(timeT[2])
-    elif typeSi == 'Y/m/d H-M-S':
-        dateDT = StringDT.split(' ')
-        dateD = dateDT[0].split('/')
-        timeT = dateDT[1].split('-')
-        yr = int(dateD[0])
-        mo = int(dateD[1])
-        da = int(dateD[2])
-        ho = int(timeT[0])
-        mi = int(timeT[1])
-        se = int(timeT[2])
-    elif typeSi == 'WRFdatetime':
-        yr = int(StringDT[0])*1000 + int(StringDT[1])*100 + int(StringDT[2])*10 +    \
-          int(StringDT[3])
-        mo = int(StringDT[5])*10 + int(StringDT[6])
-        da = int(StringDT[8])*10 + int(StringDT[9])
-        ho = int(StringDT[11])*10 + int(StringDT[12])
-        mi = int(StringDT[14])*10 + int(StringDT[15])
-        se = int(StringDT[17])*10 + int(StringDT[18])
-    else:
-        print errormsg
-        print '  ' + fname + ': type of String input date "' + typeSi +              \
-          '" not ready !!!!'
-        quit(-1)
-    if typeSo == 'matYmdHMS':
-        dateYmdHMS = np.zeros((6), dtype=int)
-        dateYmdHMS[0] =  yr
-        dateYmdHMS[1] =  mo
-        dateYmdHMS[2] =  da
-        dateYmdHMS[3] =  ho
-        dateYmdHMS[4] =  mi
-        dateYmdHMS[5] =  se
-    elif typeSo == 'YmdHMS':
-        dateYmdHMS = str(yr).zfill(4) + str(mo).zfill(2) + str(da).zfill(2) +        \
-          str(ho).zfill(2) + str(mi).zfill(2) + str(se).zfill(2)
-    elif typeSo == 'Y-m-d_H:M:S':
-        dateYmdHMS = str(yr).zfill(4) + '-' + str(mo).zfill(2) + '-' +               \
-          str(da).zfill(2) + '_' + str(ho).zfill(2) + ':' + str(mi).zfill(2) + ':' + \
-          str(se).zfill(2)
-    elif typeSo == 'Y-m-d H:M:S':
-        dateYmdHMS = str(yr).zfill(4) + '-' + str(mo).zfill(2) + '-' +               \
-          str(da).zfill(2) + ' ' + str(ho).zfill(2) + ':' + str(mi).zfill(2) + ':' + \
-          str(se).zfill(2)
-    elif typeSo == 'Y/m/d H-M-S':
-        dateYmdHMS = str(yr).zfill(4) + '/' + str(mo).zfill(2) + '/' +               \
-          str(da).zfill(2) + ' ' + str(ho).zfill(2) + '-' + str(mi).zfill(2) + '-' + \
-          str(se).zfill(2)
-    elif typeSo == 'WRFdatetime':
-        dateYmdHMS = []
-        yM = yr/1000
-        yC = (yr-yM*1000)/100
-        yD = (yr-yM*1000-yC*100)/10
-        yU = yr-yM*1000-yC*100-yD*10
-        mD = mo/10
-        mU = mo-mD*10
-        dD = da/10
-        dU = da-dD*10
-        hD = ho/10
-        hU = ho-hD*10
-        miD = mi/10
-        miU = mi-miD*10
-        sD = se/10
-        sU = se-sD*10
-        dateYmdHMS.append(str(yM))
-        dateYmdHMS.append(str(yC))
-        dateYmdHMS.append(str(yD))
-        dateYmdHMS.append(str(yU))
-        dateYmdHMS.append('-')
-        dateYmdHMS.append(str(mD))
-        dateYmdHMS.append(str(mU))
-        dateYmdHMS.append('-')
-        dateYmdHMS.append(str(dD))
-        dateYmdHMS.append(str(dU))
-        dateYmdHMS.append('_')
-        dateYmdHMS.append(str(hD))
-        dateYmdHMS.append(str(hU))
-        dateYmdHMS.append(':')
-        dateYmdHMS.append(str(miD))
-        dateYmdHMS.append(str(miU))
-        dateYmdHMS.append(':')
-        dateYmdHMS.append(str(sD))
-        dateYmdHMS.append(str(sU))
-    else:
-        print errormsg
-        print '  ' + fname + ': type of output date "' + typeSo + '" not ready !!!!'
-        quit(-1)
-    return dateYmdHMS
-def table_tex(tablevals, colnames, rownames, of):
-    """ Function to write into a LaTeX tabukar from a table of values
-      tablevals = (ncol nrow) of values
-      colnames = list with ncol labels for the columns
-      rownames = list with nrow labels for the rows
-      of= object ASCII file to write the table
-    """
-    errormsg = 'ERROR -- error -- ERROR -- error'
-    fname = 'table_tex'
-    Ncol = tablevals.shape[0]
-    Nrow = tablevals.shape[1]
-    if Ncol + 1 != len(colnames):
-        print errormsg
-        print '  ' + fname + ': wrong number of column names!!'
-        print '    data has:', Ncol, 'and you provided:', len(colnames)
-        print '    remember to provide one more for the label of the rows!'
-        print '    you provide:',colnames
-        quit(-1)
-    if Nrow != len(rownames):
-        print errormsg
-        print '  ' + fname + ': wrong number of row names!!'
-        print '    data has:', Nrow, 'and you provided:', len(rownames)
-        print '    you provide:',rownames
-        quit(-1)
-    colcs = ''
-    colns = ''
-    for icol in colnames:
-        colcs = colcs + 'c'
-        if icol == colnames[0]:
-            colns = ' {\\bfseries{' + icol + '}}'
-        else:
-            colns = colns + ' & {\\bfseries{' + icol + '}}'
-    of.write('\n')
-    of.write("%Automatically written file from function '" + fname + "'\n")
-    of.write('\\begin{tabular}{l'+colcs+'}\n')
-    of.write(colns + ' \\\\ \\hline\n')
-    ir = 0
-    for irow in rownames:
-        rowns = '{\\bfseries{' + irow + '}}'
-        for ic in range(Ncol):
-            rowns = rowns + ' & ' + str(tablevals[ic,ir])
-        rowns = rowns + ' \\\\'
-        of.write(rowns + '\n')
-        ir = ir + 1
-    of.write('\\end{tabular}\n')
-    return
-def radius_dist(dx,dy,ptx,pty):
-    """ Function to generate a matrix with the distance at a given point
-    radius_dist(dx,dy,ptx,pty)
-      [dx/y]: dimension of the matrix
-      [ptx/y]: grid point coordinates of the point
-    >>> radius_dist(3,5,2,2)
-    [[ 1.41421356  1.          1.41421356  2.23606798  3.16227766]
-     [ 1.          0.          1.          2.          3.        ]
-     [ 1.41421356  1.          1.41421356  2.23606798  3.16227766]]
-    """
-    fname = 'radius_dist'
-    if ptx < 0 or ptx > dx-1 or pty < 0 or pty > dy-1:
-        print errormsg
-        print '  ' + fname + ': wrong point coordinates:',dx,',',dy,'for matrix;',dx \
-         ,'x',dy
-        quit(-1)
-    xdist =  np.zeros((dx,dy), dtype=np.float)
-    ydist =  np.zeros((dx,dy), dtype=np.float)
-    dist =  np.zeros((dx,dy), dtype=np.float)
-    for ix in range(dx):
-        xdist[ix,:] = np.float(ix-ptx)
-    for iy in range(dy):
-        ydist[:,iy] = np.float(iy-pty)
-    dist = np.sqrt(xdist*xdist + ydist*ydist)
-    return dist
-def lonlat2D(lon,lat):
-    """ Function to return lon, lat 2D matrices from any lon,lat matrix
-      lon= matrix with longitude values
-      lat= matrix with latitude values
-    """
-    fname = 'lonlat2D'
-    if len(lon.shape) != len(lat.shape):
-        print errormsg
-        print '  ' + fname + ': longitude values with shape:', lon.shape,            \
-          'is different that latitude values with shape:', lat.shape, '(dif. size) !!'
-        quit(-1)
-    if len(lon.shape) == 3:
-        lonvv = lon[0,:,:]
-        latvv = lat[0,:,:]
-    elif len(lon.shape) == 2:
-        lonvv = lon[:]
-        latvv = lat[:]
-    elif len(lon.shape) == 1:
-        lonlatv = np.meshgrid(lon[:],lat[:])
-        lonvv = lonlatv[0]
-        latvv = lonlatv[1]
-    return lonvv, latvv
 def compute_tevolboxtraj(values, ncfile, varn):
 …
 #compute_tevolboxtraj('001/trajectory.dat@0,XLONG,XLAT,ZNU,Times,wrf,3,3',       \
 #  '001/full_concatenated.nc', 'PSFC')
-def interpolate_locs(locs,coords,kinterp):
-    """ Function to provide interpolate locations on a given axis
-    interpolate_locs(locs,axis,kinterp)
-      locs= locations to interpolate
-      coords= axis values with the reference of coordinates
-      kinterp: kind of interpolation
-        'lin': linear
-    >>> coordinates = np.arange((10), dtype=np.float)
-    >>> values = np.array([-1.2, 2.4, 5.6, 7.8, 12.0])
-    >>> interpolate_locs(values,coordinates,'lin')
-    [ -1.2   2.4   5.6   7.8  13. ]
-    >>> coordinates[0] = 0.5
-    >>> coordinates[2] = 2.5
-    >>> interpolate_locs(values,coordinates,'lin')
-    [ -3.4          1.93333333   5.6          7.8         13.        ]
-    """
-    fname = 'interpolate_locs'
-    if locs == 'h':
-        print fname + '_____________________________________________________________'
-        print interpolate_locs.__doc__
-        quit()
-    Nlocs = locs.shape[0]
-    Ncoords = coords.shape[0]
-    dcoords = coords[Ncoords-1] - coords[0]
-    intlocs = np.zeros((Nlocs), dtype=np.float)
-    minc = np.min(coords)
-    maxc = np.max(coords)
-    for iloc in range(Nlocs):
-        for icor in range(Ncoords-1):
-            if locs[iloc] < minc and dcoords > 0.:
-                a = 0.
-                b = 1. / (coords[1] - coords[0])
-                c = coords[0]
-            elif locs[iloc] > maxc and dcoords > 0.:
-                a = (Ncoords-1)*1.
-                b = 1. / (coords[Ncoords-1] - coords[Ncoords-2])
-                c = coords[Ncoords-2]
-            elif locs[iloc] < minc and dcoords < 0.:
-                a = (Ncoords-1)*1.
-                b = 1. / (coords[Ncoords-1] - coords[Ncoords-2])
-                c = coords[Ncoords-2]
-            elif locs[iloc] > maxc and dcoords < 0.:
-                a = 0.
-                b = 1. / (coords[1] - coords[0])
-                c = coords[0]
-            elif locs[iloc] >= coords[icor] and locs[iloc] < coords[icor+1] and dcoords > 0.:
-                a = icor*1.
-                b = 1. / (coords[icor+1] - coords[icor])
-                c = coords[icor]
-                print coords[icor], locs[iloc], coords[icor+1], ':', icor, '->', a, b
-            elif locs[iloc] <= coords[icor] and locs[iloc] > coords[icor+1] and dcoords < 0.:
-                a = icor*1.
-                b = 1. / (coords[icor+1] - coords[icor])
-                c = coords[icor]
-        if kinterp == 'lin':
-            intlocs[iloc] = a + (locs[iloc] - c)*b
-        else:
-            print errormsg
-            print '  ' + fname + ": interpolation kind '" + kinterp + "' not ready !!!!!"
-            quit(-1)
-    return intlocs
-def vertical_interpolation2D(varb,vart,zorigvals,znewval,kinterp):
-    """ Function to vertically integrate a 3D variable
-    vertical_interpolation2D(varb,vart,zorigvals,znewval)
-      varb= values at the base of the interval of interpolation (2D field)
-      vart= values at the top of the interval of interpolation (2D field)
-      zorigvals= pair of original values (2, 2D field)
-      znewval= new value to interpolate
-      Possible cases:
-        zorigvals[0,:,:] <= znewval < zorigvals[1,:,:]
-        znewval <= zorigvals[0,:,:] < zorigvals[1,:,:]
-        zorigvals[0,:,:] < zorigvals[1,:,:] <= znewval
-      kinterp: kind of interpolation
-        'lin': linear
-    >>> dx=5
-    >>> dy=7
-    >>> vals1 = np.ones((dy,dx), dtype=np.float)
-    >>> vals2 = np.ones((dy,dx), dtype=np.float)*2.
-    >>> zbase = np.zeros((2,dy,dx), dtype=np.float)
-    >>> zbase[0,:,:] = 0.5
-    >>> zbase[1,:,:] = 1.
-    >>> vertical_interpolation2D(vals1,vals2, zbase, newz,'lin')
-    [[ 1.5  1.5  1.5  1.5  1.5]
-     [ 1.5  1.5  1.5  1.5  1.5]
-     [ 1.5  1.5  1.5  1.5  1.5]
-     [ 1.5  1.5  1.5  1.5  1.5]
-     [ 1.5  1.5  1.5  1.5  1.5]
-     [ 1.5  1.5  1.5  1.5  1.5]
-     [ 1.5  1.5  1.5  1.5  1.5]]
-    """
-    fname = 'vertical_interpolation2D'
-    if varb == 'h':
-        print fname + '_____________________________________________________________'
-        print vertical_interpolation2D.__doc__
-        quit()
-    newvar = np.zeros((varb.shape), dtype=np.float)
-    if kinterp == 'lin':
-##        print '  ' + fname + ' Vertical linear interpolation at',znewval
-# Standard linear interpolation (y = a + b*incz)
-#   a = zorig[0], b = (vart-varb)/(zorig[1]-zorig[0])), incz = znewval - zorig[0]
-        a = varb
-        b = np.where(zorigvals[1,:,:] == zorigvals[0,:,:], 0.,                       \
-          (vart - varb)/(zorigvals[1,:,:] - zorigvals[0,:,:]))
-        incz = np.ones((varb.shape), dtype=np.float)*znewval - zorigvals[0,:,:]
-        newvar = a + b*incz
-# Too code for not be used... but maybe?
-#        dx=varb.shape[1]
-#        dy=varb.shape[0]
-#        for j in range(dy):
-#            for i in range(dx):
-#                  if zorigvals[1,j,i] == zorigvals[0,j,i]: print 'equals!',          \
-#                    zorigvals[1,j,i], zorigvals[0,j,i],':',newvar[j,i],'@',vart[j,i]
-#                  if zorigvals[0,j,i] != zorigvals[0,j,i]: print '0 Nan',            \
-#                    zorigvals[0,j,i],':',newvar[j,i],'@',vart[j,i]
-#                  if zorigvals[1,j,i] != zorigvals[1,j,i]: print '1 Nan',            \
-#                    zorigvals[1,j,i],':',newvar[j,i],'@',vart[j,i]
-#                  if zorigvals[0,j,i] is None: print '0 None', zorigvals[0,j,i],':', \
-#                    newvar[j,i],'@',vart[j,i]
-#                  if zorigvals[1,j,i] is None: print '1 None', zorigvals[1,j,i],':', \
-#                    newvar[j,i],'@',vart[j,i]
-    else:
-        print errormsg
-        print '  ' + fname + ": interpolation kind '" + kinterp + "' not ready !!!!!"
-        quit(-1)
-    return newvar
-#dx=5
-#dy=7
-#vals1 = np.ones((dy,dx), dtype=np.float)
-#vals2 = np.ones((dy,dx), dtype=np.float)*2.
-#zbase = np.zeros((2,dy,dx), dtype=np.float)
-#zbase[0,:,:] = 0.5
-#for i in range(dx):
-#    for j in range(dy):
-#        zbase[1,j,i] = np.sqrt((j-dy/2.)**2. + (i-dx/2.)**2.) /                      \
-#          np.sqrt((dy/2.)**2.+(dy/2.)**2.) + 1.
-#zbase[1,:,:] = 1.
-#newz = 0.75
-#print vertical_interpolation2D(vals1,vals2, zbase, newz,'lin')
-def vertical_interpolation(varb,vart,zorigvals,znewval,kinterp):
-    """ Function to vertically integrate a 1D variable
-    vertical_interpolation(varb,vart,zorigvals,znewval)
-      varb= values at the base of the interval of interpolation
-      vart= values at the top of the interval of interpolation
-      zorigvals= pair of original values (2)
-      znewval= new value to interpolate
-      Possible cases:
-        zorigvals[0,:] <= znewval < zorigvals[1,:]
-        znewval <= zorigvals[0,:] < zorigvals[1,:]
-        zorigvals[0,:] < zorigvals[1,:] <= znewval
-      kinterp: kind of interpolation
-        'lin': linear
-    >>> dx=5
-    >>> vals1 = np.ones((dx), dtype=np.float)
-    >>> vals2 = np.ones((dx), dtype=np.float)*2.
-    >>> zbase = np.zeros((2,dx), dtype=np.float)
-    >>> zbase[0,:] = 0.5
-    >>> for i in range(dx):
-    >>>     zbase[1,i] = i + 1.
-    >>> newz = 0.75
-    >>> vertical_interpolation2D(vals1,vals2, zbase, newz,'lin')
-    [ 1.5         1.16666667  1.1         1.07142857  1.05555556]
-    """
-    fname = 'vertical_interpolation'
-    if varb == 'h':
-        print fname + '_____________________________________________________________'
-        print vertical_interpolation.__doc__
-        quit()
-    newvar = np.zeros((varb.shape), dtype=np.float)
-    if kinterp == 'lin':
-##        print '  ' + fname + ' Vertical linear interpolation at',znewval
-# Standard linear interpolation (y = a + b*incz)
-#   a = zorig[0], b = (vart-varb)/(zorig[1]-zorig[0])), incz = znewval - zorig[0]
-        a = varb
-        b = np.where(zorigvals[1,:] == zorigvals[0,:], 0.,                           \
-          (vart - varb)/(zorigvals[1,:] - zorigvals[0,:]))
-        incz = np.ones((varb.shape), dtype=np.float)*znewval - zorigvals[0,:]
-        newvar = a + b*incz
-    else:
-        print errormsg
-        print '  ' + fname + ": interpolation kind '" + kinterp + "' not ready !!!!!"
-        quit(-1)
-    return newvar
-def interpolate_3D(zorigvs, varv, zintvs, kint):
-    """ Function to interpolate a 3D variable
-    interpolate_3D(zintvs, varv, zorigvals)
-      zorigvs= original 3D z values
-      varv= 3D variable to interpolate
-      zintvs= new series of z values to interpolate
-      kint= kind of interpolation
-        'lin': linear
-    """
-    fname = 'interpolate_3D'
-    if zorigvs == 'h':
-        print fname + '_____________________________________________________________'
-        print interpolate_3D.__doc__
-        quit()
-    Ninv = len(zintvs)
-    dimv = varv.shape
-# Sense of the original z-variable
-    Lzorig = zorigvs.shape[0]
-    incz = zorigvs[Lzorig-1,0,0] - zorigvs[0,0,0]
-    varin = np.zeros((Ninv,dimv[1],dimv[2]), dtype=np.float)
-    for ilev in range(Ninv):
-#        print ' vertical interpolate level: ',zintvs[ilev]
-        valinf = np.zeros((dimv[1], dimv[2]), dtype=np.float)
-        valsup = np.zeros((dimv[1], dimv[2]), dtype=np.float)
-        zorigv = np.zeros((2,dimv[1], dimv[2]), dtype=np.float)
-        for j in range(valinf.shape[0]):
-            for i in range(valinf.shape[1]):
-                if np.min(zorigvs[:,j,i]) > zintvs[ilev] and incz > 0.:
-                    valinf[j,i] = varv[0,j,i]
-                    valsup[j,i] = varv[1,j,i]
-                    zorigv[0,j,i] = zorigvs[0,j,i]
-                    zorigv[1,j,i] = zorigvs[1,j,i]
-                elif np.max(zorigvs[:,j,i]) < zintvs[ilev] and incz < 0.:
-                    valinf[j,i] = varv[0,j,i]
-                    valsup[j,i] = varv[1,j,i]
-                    zorigv[0,j,i] = zorigvs[0,j,i]
-                    zorigv[1,j,i] = zorigvs[1,j,i]
-                elif np.max(zorigvs[:,j,i]) < zintvs[ilev] and incz > 0.:
-                    valinf[j,i] = varv[Lzorig-2,j,i]
-                    valsup[j,i] = varv[Lzorig-1,j,i]
-                    zorigv[0,j,i] = zorigvs[Lzorig-2,j,i]
-                    zorigv[1,j,i] = zorigvs[Lzorig-1,j,i]
-                elif np.min(zorigvs[:,j,i]) > zintvs[ilev] and incz < 0.:
-                    valinf[j,i] = varv[Lzorig-2,j,i]
-                    valsup[j,i] = varv[Lzorig-1,j,i]
-                    zorigv[0,j,i] = zorigvs[Lzorig-2,j,i]
-                    zorigv[1,j,i] = zorigvs[Lzorig-1,j,i]
-#                    print 'top: ',i,j,':',zorigv[0,j,i], zintvs[ilev], zorigv[1,j,i]
-                else:
-                    for z in range(Lzorig-1):
-                        if (zorigvs[z,j,i]-zintvs[ilev])*(zorigvs[z+1,j,i]-          \
-                          zintvs[ilev]) <= 0.:
-                            valinf[j,i] = varv[z,j,i]
-                            valsup[j,i] = varv[z+1,j,i]
-                            zorigv[0,j,i] = zorigvs[z,j,i]
-                            zorigv[1,j,i] = zorigvs[z+1,j,i]
-                            break
-#        print 'zorigvs: ',zorigvs[:,0,0]
-#        print '  valinf:', valinf[0,0]
-#        print '  valsup:', valsup[0,0]
-#        print '  zorigv 0:',zorigv[0,0,0]
-#        print '  zorigv 1:',zorigv[1,0,0]
-        varin[ilev,:,:] = vertical_interpolation2D(valinf,valsup,zorigv,             \
-          zintvs[ilev], kint)
-#        print '  varin:',varin[ilev,0,0]
-#        quit()
-#    quit()
-    return varin
-def interpolate_2D(zorigvs0, varv0, zdim, zintvs, kint):
-    """ Function to interpolate a 2D variable
-    interpolate_2D(zintvs, varv, zorigvals)
-      zorigvs= original 2D z values
-      varv= 2D variable to interpolate
-      zdim= number of the dimension with the z-values
-      zintvs= new series of z values to interpolate
-      kint= kind of interpolation
-        'lin': linear
-    """
-    fname = 'interpolate_2D'
-    if zorigvs0 == 'h':
-        print fname + '_____________________________________________________________'
-        print interpolate_2D.__doc__
-        quit()
-    Ninv = len(zintvs)
-    if zdim == 0:
-        varv = varv0
-        zorigvs = zorigvs0
-    else:
-        varv = np.transpose(varv0)
-        zorigvs = np.transpose(zorigvs0)
-    dimv = varv.shape
-# Sense of the original z-variable
-    Lzorig = zorigvs.shape[0]
-    incz = zorigvs[Lzorig-1,0] - zorigvs[0,0]
-    varin = np.zeros((Ninv,dimv[1]), dtype=np.float)
-    for ilev in range(Ninv):
-        valinf = np.zeros((dimv[1]), dtype=np.float)
-        valsup = np.zeros((dimv[1]), dtype=np.float)
-        zorigv = np.zeros((2,dimv[1]), dtype=np.float)
-        for i in range(valinf.shape[0]):
-            if np.min(zorigvs[:,i]) > zintvs[ilev] and incz > 0.:
-                valinf[i] = varv[0,i]
-                valsup[i] = varv[1,i]
-                zorigv[0,i] = zorigvs[0,i]
-                zorigv[1,i] = zorigvs[1,i]
-            elif np.max(zorigvs[:,i]) < zintvs[ilev] and incz < 0.:
-                valinf[i] = varv[0,i]
-                valsup[i] = varv[1,i]
-                zorigv[0,i] = zorigvs[0,i]
-                zorigv[1,i] = zorigvs[1,i]
-            elif np.max(zorigvs[:,i]) < zintvs[ilev] and incz > 0.:
-                valinf[i] = varv[Lzorig-2,i]
-                valsup[i] = varv[Lzorig-1,i]
-                zorigv[0,i] = zorigvs[Lzorig-2,i]
-                zorigv[1,i] = zorigvs[Lzorig-1,i]
-            elif np.min(zorigvs[:,i]) > zintvs[ilev] and incz < 0.:
-                valinf[i] = varv[Lzorig-2,i]
-                valsup[i] = varv[Lzorig-1,i]
-                zorigv[0,i] = zorigvs[Lzorig-2,i]
-                zorigv[1,i] = zorigvs[Lzorig-1,i]
-            else:
-                for z in range(Lzorig-1):
-                    if (zorigvs[z,i]-zintvs[ilev])*(zorigvs[z+1,i]-                  \
-                      zintvs[ilev]) <= 0.:
-                        valinf[i] = varv[z,i]
-                        valsup[i] = varv[z+1,i]
-                        zorigv[0,i] = zorigvs[z,i]
-                        zorigv[1,i] = zorigvs[z+1,i]
-                        break
-        varin[ilev,:] = vertical_interpolation(valinf,valsup,zorigv,                 \
-          zintvs[ilev], kint)
-    if zdim == 0:
-        varin0 = varin
-    else:
-        varin0 = np.transpose(varin)
-    return varin0
-def list_toVec(listv, kind):
-    """ Function to transform from a list of values to a vector
-    list_toVec(listv, kind)
-      listv= list with the values
-      kind= kind of values
-      >>> list_toVec(['1', 2, '3', 4, 5, '6.9', 7, 9], 'npfloat')
-      [ 1.   2.   3.   4.   5.   6.9  7.   9. ]
-    """
-    fname = 'list_toVec'
-    if listv == 'h':
-        print fname + '_____________________________________________________________'
-        print list_toVec.__doc__
-        quit()
-    Nvals = len(listv)
-    if kind == 'npfloat':
-        vecv = np.zeros((Nvals), dtype=np.float)
-        for iv in range(Nvals):
-            vecv[iv] = np.float(listv[iv])
-    else:
-        print errormsg
-        print '  ' + fname + ': type of value "' + kind + '" not ready!!'
-        quit(-1)
-    return vecv
-def running_mean(values, run):
-    """ Function to compute a running mean of a series of values
-      running_mean(vals, int)
-      [vals]: series of values
-      [run]: interval to use to compute the running mean
-        NOTE: resultant size would be the same but with None except at [run/2,size(vals)-run/2]
-      >>> running_mean(np.arange(100),10)
-      [None None None None None 4.5 5.5 6.5 7.5 8.5 9.5 10.5 11.5 12.5 13.5 14.5
-.5 16.5 17.5 18.5 19.5 20.5 21.5 22.5 23.5 24.5 25.5 26.5 27.5 28.5 29.5
-.5 31.5 32.5 33.5 34.5 35.5 36.5 37.5 38.5 39.5 40.5 41.5 42.5 43.5 44.5
-.5 46.5 47.5 48.5 49.5 50.5 51.5 52.5 53.5 54.5 55.5 56.5 57.5 58.5 59.5
-.5 61.5 62.5 63.5 64.5 65.5 66.5 67.5 68.5 69.5 70.5 71.5 72.5 73.5 74.5
-.5 76.5 77.5 78.5 79.5 80.5 81.5 82.5 83.5 84.5 85.5 86.5 87.5 88.5 89.5
-.5 91.5 92.5 93.5 None None None None None]
-    """
-    fname = 'running_mean'
-    if values == 'h':
-        print fname + '_____________________________________________________________'
-        print running_mean.__doc__
-        quit()
-    runmean = np.ones((len(values)), dtype=np.float)*fillValue
-    for ip in range(run/2,len(values)-run/2):
-        runmean[ip] = np.mean(values[ip-run/2:ip+run/2])
-    runmean = np.where(runmean == fillValue, None, runmean)
-    return runmean
 def slice_variable(varobj, dimslice):
 …
 #WRF_CFlon_creation('longitude', 'wrfout_d01_1980-03-01_00:00:00_Time_B0-E48-I1.nc', 'XLONG')
 def dimVar_creation(values, ncfile):
     """ Function to add a 1D variable with the size of a given dimension in a file
 …
 #dimVar_creation('bottom_top', 'wrfout_d01_1979-12-01_00:00:00_south_north_B3-E3-I1_west_east_B26-E26-I1.nc')
-def files_folder(folder,headfile):
-    """ Function to retrieve a list of files from a folder [fold] and head [headfile]
-    >>> files_folder('/media/data2/etudes/WRF_LMDZ/WL_HyMeX_HighRes/medic950116/',
-      'wrfout_d03')
-    ['wrfout_d03_1995-01-13_02:00:00', 'wrfout_d03_1995-01-13_04:00:00', 'wrfout_d03_1995-01-13_06:00:00',
-     'wrfout_d03_1995-01-13_06:15:00', 'wrfout_d03_1995-01-13_08:15:00', 'wrfout_d03_1995-01-13_10:15:00',
-     'wrfout_d03_1995-01-13_12:15:00', 'wrfout_d03_1995-01-13_14:15:00', 'wrfout_d03_1995-01-13_16:15:00',
-     (...)
-     'wrfout_d03_1995-01-17_18:15:00', 'wrfout_d03_1995-01-17_20:15:00', 'wrfout_d03_1995-01-17_22:15:00']
-    """
-    import subprocess as sub
-    fname = 'files_folder'
-    if folder == 'h':
-        print fname + '_____________________________________________________________'
-        print files_folder.__doc__
-        quit()
-#    ins = folder + "/" + headfile + "*"
-    ins = folder + "/"
-    print 'ins:',ins
-    files = sub.Popen(["/bin/ls","-1",ins], stdout=sub.PIPE)
-    fileslist = files.communicate()
-    listfiles0 = str(fileslist).split('\\n')
-# Filtering output
-    Lheader=len(headfile)
-    listfiles = []
-    for ifile in listfiles0:
-        if ifile[0:Lheader] == headfile:
-            listfiles.append(ifile)
-    return listfiles
 def netcdf_fold_concatenation(values, ncfile, varn):
 …
 #netcdf_concatenation('mean_pluc.nc@all|mean_dtcon.nc@all')
-def count_cond(var, val, con):
-    """ Function to count values of a variable which attain a condition
-      [var]= variable
-      [val]= value
-      [con]= statistics/condition
-        'eq': equal than [val]
-        'neq': not equal than [val]
-        'lt': less than [val]
-        'le': less and equal than [val]
-        'gt': greater than [val]
-        'ge': greater and equal than [val]
-        'in': inside the range [val[0]] <= [var] <= [val[1]] (where [val]=[val1, val2]
-        'out': inside the range [val[0]] > [var]; [val[1]] < [var] (where [val]=[val1, val2]
-      >>> vals = np.arange(100)
-      >>> count_cond(vals,50,'eq')
-      >>> count_cond(vals,50,'neq')
-      >>> count_cond(vals,50,'lt')
-      >>> count_cond(vals,50,'le')
-      >>> count_cond(vals,50,'gt')
-      >>> count_cond(vals,50,'ge')
-      >>> count_cond(vals,[25,75],'in')
-      >>> count_cond(vals,[25,75],'out')
-    """
-    fname = 'count_cond'
-    cons = ['eq', 'neq', 'lt', 'le', 'gt', 'ge', 'in', 'out']
-    if con == 'eq':
-        Nvals = np.sum(var == val)
-    elif con == 'neq':
-        Nvals = np.sum(var != val)
-    elif con == 'lt':
-        Nvals = np.sum(var < val)
-    elif con == 'le':
-        Nvals = np.sum(var <= val)
-    elif con == 'gt':
-        Nvals = np.sum(var > val)
-    elif con == 'ge':
-        Nvals = np.sum(var >= val)
-    elif con == 'in':
-        if len(val) != 2:
-            print errormsg
-            print '  ' + fname + ": for condition '" + con + "' we must have two " + \
-              " values!!"
-            print '  we have:',  val
-            quit(-1)
-        Nvals = np.sum((var >= val[0])*(var <= val[1]))
-    elif con == 'out':
-        if len(val) != 2:
-            print errormsg
-            print '  ' + fname + ": for condition '" + con + "' we must have two " + \
-              " values!!"
-            print '  we have:',  val
-            quit(-1)
-        Nvals = np.sum((var < val[0])+(var > val[1]))
-    else:
-        print errormsg
-        print '  ' + fname + ": condition '" + con + "' not ready!!"
-        print '  only ready:', cons
-        quit(-1)
-    return Nvals
 def field_stats(values, ncfile, varn):
     """ Function to retrieve statistics from a field
 …
 #field_stats('full', 'geo_em.d01.nc', 'HGT_M')
-def get_namelist_vars(values, namelist):
-    """ Function to get namelist-like  values ([varname] = [value])
-    get_namelist_vars(namelist)
-      values= [sectionname],[kout]
-        [sectionname]: name of the section from which one want the values ('all', for all)
-        [kout]: kind of output
-          'tex3': printed output as LaTeX table of three columns \verb+namelist_name+ & value
-          'column': printed output as namelist_name value
-          'dict': as two python dictionary object (namelistname, value and namelistname, sectionname)
-      namelist= namelist_like file to retrieve values
-    >>> get_namelist_vars('geogrid,dic','/home/lluis/etudes/domains/medic950116/namelist.wps')
-    {'e_sn': '32, 97,', 'stand_lon': '0.', 'opt_geogrid_tbl_path': "'./'", 'geog_data_path': "'/home/lluis/DATA/WRF/geog'", 'pole_lat': '90.0', 'ref_lat': '35.0,', 'map_proj': "'lat-lon',", 'parent_id': '1, 1,', 'geog_data_res': "'2m','2m',", 'e_we': '32, 112,', 'dx': '0.35,', 'dy': '0.35,', 'parent_grid_ratio': '1, 3,', 'pole_lon': '0.0', 'ref_lon': '20,', 'j_parent_start': '1, 17,', 'i_parent_start': '1, 31,'}
-    """
-    fname = 'get_namelist_vars'
-# List of characters which split pairs name/value
-    valuessep = ['=']
-# List of characters which indicate a comment
-    commentchars = ['#']
-    if values == 'h':
-        print fname + '_____________________________________________________________'
-        print get_namelist_vars.__doc__
-        quit()
-    expectargs = '[sectionname],[kout]'
-    check_arguments(fname,values,expectargs,',')
-    secname = values.split(',')[0]
-    kout = values.split(',')[1]
-    if not os.path.isfile(namelist):
-        print errormsg
-        print '  ' + fname + ": namelist file '" + namelist + "' does not exist !!"
-        quit(-1)
-    ncml = open(namelist, 'r')
-    sections = {}
-    namelistvals = {}
-    namelistsecs = {}
-    sectionnames = []
-    namessec = []
-    allnames = []
-    namelistvalssec = {}
-    namelistsecssec = {}
-    nmlname = ''
-    sectionname = ''
-    for line in ncml:
-        linevals = reduce_spaces(line)
-        Nvals = len(linevals)
-        if Nvals >= 1 and linevals[0][0:1] == '&':
-            if len(sectionnames) > 1:
-                sections[sectionname] = namessec
-            sectionname = linevals[0][1:len(linevals[0])+1]
-#            print '    ' + fname + ": new section '" + sectionname + "' !!!"
-            sectionnames.append(sectionname)
-            namessec = []
-            nmlname = ''
-        elif Nvals >= 1 and not searchInlist(commentchars,linevals[0][0:1]):
-            if Nvals >= 3 and searchInlist(valuessep,linevals[1]):
-                nmlname = linevals[0]
-                nmlval = numVector_String(linevals[2:Nvals],' ')
-            elif Nvals == 1:
-                for valsep in valuessep:
-                    if linevals[0].find(valsep) != -1:
-                        nmlname = linevals[0].split(valsep)[0]
-                        nmlval = linevals[0].split(valsep)[1]
-                        break
-            elif Nvals == 2:
-                for valsep in valuessep:
-                    if linevals[0].find(valsep) != -1:
-                        nmlname = linevals[0].split(valsep)[0]
-                        nmlval = linevals[1]
-                        break
-                    elif linevals[1].find(valsep) != -1:
-                        nmlname = linevals[0]
-                        nmlval = linevals[1].split(valsep)[0]
-                        break
-            else:
-                print warnmsg
-                print '  ' + fname + ': wrong number of values', Nvals,              \
-                  'in the namelist line!'
-                print '    line: ',line
-                print '    line values:',linevals
-#                quit(-1)
-            namelistvals[nmlname] = nmlval
-            namelistsecs[nmlname] = sectionname
-            namessec.append(nmlname)
-            allnames.append(nmlname)
-    if len(sectionname) > 1:
-        sections[sectionname] = namessec
-    if secname != 'all':
-        if not searchInlist(sections.keys(),secname):
-            print errormsg
-            print '  ' + fname + ": section '" + values + "' does not exist !!"
-            print '    only found:',sectionnames
-            quit(-1)
-        namestouse = []
-        for nml in allnames:
-            for nnml in sections[secname]:
-                namelistvalssec[nnml] = namelistvals[nnml]
-                namelistsecssec[nnml] = secname
-                if nml == nnml: namestouse.append(nml)
-    else:
-        namestouse = allnames
-        namelistvalssec = namelistvals
-        namelistsecssec = namelistsecs
-    if kout == 'tex3':
-        ofile='get_namelist_vars_3col.tex'
-        fobj = open(ofile, 'w')
-        vals = namestouse
-        Nvals = len(vals)
-        Nvals3 = int(Nvals/3)
-        latextab = '\\begin{center}\n\\begin{tabular}{lclclc}\n'
-        for il in range(2):
-            latextab = latextab + '{\\bfseries{name}} & {\\bfseries{value}} & '
-        latextab= latextab+ '{\\bfseries{name}} & {\\bfseries{value}} \\\\ \\hline\n'
-        if np.mod(Nvals,3) != 0:
-            Nvals0 = Nvals - np.mod(Nvals,3)
-        else:
-            Nvals0 = Nvals
-        for ival in range(0,Nvals0,3):
-            line = ''
-            print '  ival:',ival
-            for il in range(2):
-                line = line + '\\verb+' + vals[ival+il] + '+ & ' +                   \
-                   namelistvalssec[vals[ival+il]].replace('_','\\_') +' & '
-            line = line + '\\verb+' + vals[ival+2] + '+ & ' +                       \
-               namelistvalssec[vals[ival+2]].replace('_','\\_') + ' \\\\\n'
-            latextab = latextab + line
-        latextab = latextab + '%not multiple of three!!!!\n'
-        print 'mod:', np.mod(Nvals,3),Nvals0,Nvals
-        if np.mod(Nvals,3) != 0:
-            ival = Nvals0
-            line = ''
-            for il in range(np.mod(Nvals,3)):
-                print 'ival:',ival + il
-                line = line + '\\verb+' + vals[ival+il] + '+ & ' +                   \
-                      namelistvalssec[vals[ival+il]].replace('_','\\_') + ' & '
-            for il in range(2-np.mod(Nvals,3)):
-                line = line + ' & & '
-            latextab = latextab + line + ' & \\\\\n'
-        latextab = latextab + '\\end{tabular}\n\\end{center}\n'
-#        print latextab
-        fobj.write(latextab)
-        fobj.close()
-        print fname + ": successful writen '" + ofile + "' LaTeX tale file !!"
-        return
-    elif kout == 'column':
-        for dictv in namestouse:
-            print dictv + ' = ' + namelistvalssec[dictv]
-        return
-    elif kout == 'dict':
-        return namelistvalssec, namelistsecssec
-    else:
-        print errormsg
-        print '  ' + fname + ": kind of output '" + kout + "' not ready!!!"
-        quit(-1)
-    return
-def WRF_d0Nref(values, geofold):
-    """ Function for the generation of an extra WRF domain from a given one
-    field_stats(values, ncfile, varn)
-      [values]= [namelist],[minLON],[minLAT],[maxLON],[maxLAT],[Ndomref],[factor],[geogres]
-        [namelist]: namelist.wps to use
-        [minLON],[minLAT]: minimum longitude and latitude
-        [maxLON],[maxLAT]: maximum longitude and latitude
-        [Ndomref]: domain reference number
-        [factor]: resolution factor with respect d01
-        [geogres]: which topography resolution to use
-      [geofold]= folder where to found the geo_em.d0[N-1].nc file
-    """
-    fname='WRF_d0Nref'
-    if values == 'h':
-        print fname + '_____________________________________________________________'
-        print WRF_d0Nref.__doc__
-        quit()
-    expectargs = '[namelist],[minLON],[minLAT],[maxLON],[maxLAT],[Ndomref],[factor],'
-    expectargs = expectargs + '[geogres]'
-    check_arguments(fname,values,expectargs,',')
-    namelist = values.split(',')[0]
-    minLON = np.float(values.split(',')[1])
-    minLAT = np.float(values.split(',')[2])
-    maxLON = np.float(values.split(',')[3])
-    maxLAT = np.float(values.split(',')[4])
-    Ndomref = int(values.split(',')[5])
-    factor = int(values.split(',')[6])
-    geogres = values.split(',')[7]
-    onml = 'namelist.wps_d0' + str(Ndomref+1)
-# Namelist variables to which new value is a copy of the last one
-    repeatvalue = ['start_date' ,'end_date']
-# Namelist variables to which value is increased by 1 (single value)
-    add1 = ['max_dom']
-# Namelist variables to which new value is an increase by 1 of the last one
-    add1value = ['parent_id']
-#    ncfile = geofold + '/geo_em.d' + zeros(Ndomref,2) + '.nc'
-    ncfile = geofold + '/geo_em.d0' + str(Ndomref) + '.nc'
-    if not os.path.isfile(ncfile):
-        print errormsg
-        print '  ' + fname + ': geo_em file "' + ncfile + '" does not exist !!'
-        quit(-1)
-    ncobj = NetCDFFile(ncfile, 'r')
-    nlon='XLONG_M'
-    nlat='XLAT_M'
-    lonobj = ncobj.variables[nlon]
-    latobj = ncobj.variables[nlat]
-    lon = lonobj[0,:,:]
-    lat = latobj[0,:,:]
-    ncobj.close()
-    diffmin = np.sqrt((lon - minLON)**2. + (lat - minLAT)**2.)
-    diffmax = np.sqrt((lon - maxLON)**2. + (lat - maxLAT)**2.)
-    posdiffmin = index_mat(diffmin,np.min(diffmin))
-    posdiffmax = index_mat(diffmax,np.min(diffmax))
-#    print 'min vals. min:',np.min(diffmin),'max:',np.min(diffmax)
-#    print 'min:', minLON, ',', minLAT, ':', posdiffmin, '.', lon[tuple(posdiffmin)], \
-#      ',', lat[tuple(posdiffmin)]
-#    print 'max:', maxLON, ',', maxLAT, ':', posdiffmax, '.', lon[tuple(posdiffmax)], \
-#      ',', lat[tuple(posdiffmax)]
-    print '  ' + fname +': coordinates of the SW corner:', posdiffmin, ' NE corner:',\
-      posdiffmax
-    nmlobj = open(namelist, 'r')
-    onmlobj = open(onml, 'w')
-    for line in nmlobj:
-        listline = reduce_spaces(line)
-        nvals = len(listline)
-        if nvals > 0:
-            if searchInlist(add1,listline[0]):
-                nmlval = int(listline[nvals-1].replace(',','')) + 1
-                onmlobj.write(' ' + listline[0] + ' = ' + str(nmlval) + ',\n')
-            elif searchInlist(add1value,listline[0]):
-                if Ndomref == 1:
-                    nmlval = 1
-                else:
-                    nmlval = int(listline[nvals-1].replace(',','')) + 1
-                onmlobj.write(' ' + listline[0] + ' = ' +                            \
-                  numVector_String(listline[2:nvals+1], ' ') + ' ' +str(nmlval)+',\n')
-            elif searchInlist(repeatvalue,listline[0]):
-                nmlval = listline[nvals-1]
-                onmlobj.write(' ' + listline[0] + ' = ' +                            \
-                  numVector_String(listline[2:nvals+1],' ') + ' ' + str(nmlval)+'\n')
-            elif listline[0] == 'parent_grid_ratio':
-                nmlval = factor
-                onmlobj.write(' ' + listline[0] + ' = ' +                            \
-                  numVector_String(listline[2:nvals+1],' ') + ' ' + str(nmlval)+',\n')
-            elif listline[0] == 'i_parent_start':
-                nmlval = posdiffmin[1] - 10
-                onmlobj.write(' ' + listline[0] + ' = ' +                            \
-                  numVector_String(listline[2:nvals+1],' ') + ' ' + str(nmlval)+',\n')
-            elif listline[0] == 'j_parent_start':
-                nmlval = posdiffmin[0] - 10
-                onmlobj.write(' ' + listline[0] + ' = ' +                            \
-                  numVector_String(listline[2:nvals+1],' ') + ' ' + str(nmlval)+',\n')
-            elif listline[0] == 'e_we':
-                nmlval = (posdiffmax[1] - posdiffmin[1] + 20)*factor
-                nmlval = nmlval + np.mod(nmlval,factor) + 1
-                onmlobj.write(' ' + listline[0] + ' = ' +                            \
-                  numVector_String(listline[2:nvals+1],' ') + ' ' + str(nmlval)+',\n')
-            elif listline[0] == 'e_sn':
-                nmlval = (posdiffmax[0] - posdiffmin[0] + 20)*factor
-                nmlval = nmlval + np.mod(nmlval,factor) + 1
-                onmlobj.write(' ' + listline[0] + ' = ' +                            \
-                  numVector_String(listline[2:nvals+1],' ') + ' ' + str(nmlval)+',\n')
-            elif listline[0] == 'geog_data_res':
-                nmlval = "'" + geogres + "'"
-                onmlobj.write(' ' + listline[0] + ' = ' +                            \
-                  numVector_String(listline[2:nvals+1],' ') + ' ' + str(nmlval)+',\n')
-            else:
-                onmlobj.write(line)
-        else:
-            onmlobj.write('\n')
-    nmlobj.close()
-    onmlobj.close()
-    print fname + ": successfully written new namelist '" + onml + "' !!"
-    return
 def filter_2dim(values, ncfile, varn):
 …
 #selvar('x@lon,y@lat,time@time,shan@int', '/home/lluis/PY/test.nc', 'var')
-def coincident_CFtimes(tvalB, tunitA, tunitB):
-    """ Function to make coincident times for two different sets of CFtimes
-    tvalB= time values B
-    tunitA= time units times A to which we want to make coincidence
-    tunitB= time units times B
-    >>> coincident_CFtimes(np.arange(10),'seconds since 1949-12-01 00:00:00',
-      'hours since 1949-12-01 00:00:00')
-    [     0.   3600.   7200.  10800.  14400.  18000.  21600.  25200.  28800.  32400.]
-    >>> coincident_CFtimes(np.arange(10),'seconds since 1949-12-01 00:00:00',
-      'hours since 1979-12-01 00:00:00')
-    [  9.46684800e+08   9.46688400e+08   9.46692000e+08   9.46695600e+08
-.46699200e+08   9.46702800e+08   9.46706400e+08   9.46710000e+08
-.46713600e+08   9.46717200e+08]
-    """
-    import datetime as dt
-    fname = 'coincident_CFtimes'
-    trefA = tunitA.split(' ')[2] + ' ' + tunitA.split(' ')[3]
-    trefB = tunitB.split(' ')[2] + ' ' + tunitB.split(' ')[3]
-    tuA = tunitA.split(' ')[0]
-    tuB = tunitB.split(' ')[0]
-    if tuA != tuB:
-        if tuA == 'microseconds':
-            if tuB == 'microseconds':
-                tB = tvalB*1.
-            elif tuB == 'seconds':
-                tB = tvalB*10.e6
-            elif tuB == 'minutes':
-                tB = tvalB*60.*10.e6
-            elif tuB == 'hours':
-                tB = tvalB*3600.*10.e6
-            elif tuB == 'days':
-                tB = tvalB*3600.*24.*10.e6
-            else:
-                print errormsg
-                print '  ' + fname + ": combination of time untis: '" + tuA +        \
-                  "' & '" + tuB + "' not ready !!"
-                quit(-1)
-        elif tuA == 'seconds':
-            if tuB == 'microseconds':
-                tB = tvalB/10.e6
-            elif tuB == 'seconds':
-                tB = tvalB*1.
-            elif tuB == 'minutes':
-                tB = tvalB*60.
-            elif tuB == 'hours':
-                tB = tvalB*3600.
-            elif tuB == 'days':
-                tB = tvalB*3600.*24.
-            else:
-                print errormsg
-                print '  ' + fname + ": combination of time untis: '" + tuA +        \
-                  "' & '" + tuB + "' not ready !!"
-                quit(-1)
-        elif tuA == 'minutes':
-            if tuB == 'microseconds':
-                tB = tvalB/(60.*10.e6)
-            elif tuB == 'seconds':
-                tB = tvalB/60.
-            elif tuB == 'minutes':
-                tB = tvalB*1.
-            elif tuB == 'hours':
-                tB = tvalB*60.
-            elif tuB == 'days':
-                tB = tvalB*60.*24.
-            else:
-                print errormsg
-                print '  ' + fname + ": combination of time untis: '" + tuA +        \
-                  "' & '" + tuB + "' not ready !!"
-                quit(-1)
-        elif tuA == 'hours':
-            if tuB == 'microseconds':
-                tB = tvalB/(3600.*10.e6)
-            elif tuB == 'seconds':
-                tB = tvalB/3600.
-            elif tuB == 'minutes':
-                tB = tvalB/60.
-            elif tuB == 'hours':
-                tB = tvalB*1.
-            elif tuB == 'days':
-                tB = tvalB*24.
-            else:
-                print errormsg
-                print '  ' + fname + ": combination of time untis: '" + tuA +        \
-                  "' & '" + tuB + "' not ready !!"
-                quit(-1)
-        elif tuA == 'days':
-            if tuB == 'microseconds':
-                tB = tvalB/(24.*3600.*10.e6)
-            elif tuB == 'seconds':
-                tB = tvalB/(24.*3600.)
-            elif tuB == 'minutes':
-                tB = tvalB/(24.*60.)
-            elif tuB == 'hours':
-                tB = tvalB/24.
-            elif tuB == 'days':
-                tB = tvalB*1.
-            else:
-                print errormsg
-                print '  ' + fname + ": combination of time untis: '" + tuA +        \
-                  "' & '" + tuB + "' not ready !!"
-                quit(-1)
-        else:
-            print errormsg
-            print '  ' + fname + ": time untis: '" + tuA + "' not ready !!"
-            quit(-1)
-    else:
-        tB = tvalB*1.
-    if trefA != trefB:
-        trefTA = dt.datetime.strptime(trefA, '%Y-%m-%d %H:%M:%S')
-        trefTB = dt.datetime.strptime(trefB, '%Y-%m-%d %H:%M:%S')
-        difft = trefTB - trefTA
-        diffv = difft.days*24.*3600.*10.e6 + difft.seconds*10.e6 + difft.microseconds
-        print '  ' + fname + ': different reference refA:',trefTA,'refB',trefTB
-        print '    difference:',difft,':',diffv,'microseconds'
-        if tuA == 'microseconds':
-            tB = tB + diffv
-        elif tuA == 'seconds':
-            tB = tB + diffv/10.e6
-        elif tuA == 'minutes':
-            tB = tB + diffv/(60.*10.e6)
-        elif tuA == 'hours':
-            tB = tB + diffv/(3600.*10.e6)
-        elif tuA == 'dayss':
-            tB = tB + diffv/(24.*3600.*10.e6)
-        else:
-            print errormsg
-            print '  ' + fname + ": time untis: '" + tuA + "' not ready !!"
-            quit(-1)
-    return tB
-def 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)
-.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 read_ASCIIlist(filen):
-    """ Function to build a list from an ASCII lines file
-      filen= name of the ASCII file
-    """
-    import os
-    fname = 'read_ASCIIlist'
-    if not os.path.isfile(filen):
-        print errormsg
-        print '  ' + fname + ":  ASCII types file '" + filen + "' does not exist !!"
-        quit(-1)
-    of = open(filen, 'r')
-    filevals = []
-    for linef in of:
-        filevals.append(linef.replace('\n',''))
-    of.close()
-    return filevals
-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 read_SQradii(filename='SQradii.dat', Npt=-1, maxradii=None):
-    """ Function to read the outcome of the function `squared_radial'
-      filename= Name of the outcome file of the function
-      Npt= Number of points on each grid
-      maxradii= allowed maximum radii
-      'SQradii.dat': file with the radii and paris of integers after ran squared_radial(1000)
-    """
-    fname = 'read_SQradii'
-    if not os.path.isfile(filename):
-        print errormsg
-        print '  ' + fname + ": file '" + filename + "' not found!!"
-        quit(-1)
-    infile = open(filename, 'r')
-    gridradii = {}
-    Nradii = {}
-    radii = []
-    for line in infile:
-        values = values_line(line, ' ', ['[', ']'])
-        rad = np.float(values[0])
-# Limited by a given radii
-        if maxradii is not None and rad > maxradii: break
-        Nrad = int(values[2])
-        couples = []
-        Nradgood = 0
-        if Nrad > 1:
-           Ncouples = 0
-           for ic in range(Nrad):
-               ix = int(values[3 + ic*2])
-               iy = int(values[3 + ic*2 + 1])
-               if ix <= Npt and iy <= Npt:
-                   couples.append(np.array([ix,iy]))
-                   Nradgood = Nradgood + 1
-               elif ix > Npt and iy > Npt: break
-        else:
-             ix = int(values[3])
-             iy = int(values[4])
-             if ix <= Npt and iy <= Npt:
-                 couples.append(np.array([ix,iy]))
-                 Nradgood = Nradgood + 1
-             elif ix > Npt and iy > Npt: break
-        if Nradgood > 0:
-            Nradii[rad] = Nradgood
-            finalcouples = np.zeros((Nradgood,2), dtype=int)
-            for ic in range(Nradgood):
-                finalcouples[ic,:] = couples[ic]
-            gridradii[rad] = finalcouples
-#    for rad in sorted(gridradii.keys()):
-#        print rad, gridradii[rad]
-    return gridradii
-def grid_combinations(x,y):
-    """ Function to provide all the possible grid points combination for a given pair of values
-      x,y= pair of grid points
-    >>>grid_combinations(1,2)
-    [[1, 2], [-1, 2], [-1, -2], [1, -2], [2, 1], [-2, 1], [-2, -1], [2, -1]]
-    """
-    fname = 'grid_combinations'
-    gridcomb = []
-    gridcomb.append([x,y])
-    if x != 0: gridcomb.append([-x,y])
-    if x != 0 and y != 0: gridcomb.append([-x,-y])
-    if y != 0: gridcomb.append([x,-y])
-    if x != y:
-        gridcomb.append([y,x])
-        if y != 0: gridcomb.append([-y,x])
-        if x != 0 and y != 0: gridcomb.append([-y,-x])
-        if x != 0: gridcomb.append([y,-x])
-    return gridcomb
-def radii_points(xpos,ypos,Lrad,Nang,dx,dy):
-    """ Function to provide the grid points for radial cross sections, by angle and radii and `squared' radii
-      xpos= x position of the center
-      ypos= y position of the center
-      Lrad= length of the maximum radi in grid points
-      Nang= number of equivalent angles
-      dx= x size of the domain of values
-      dy= y size of the domain of values
-      >>> radpos, SQradpos = radii_points(12,12,1,8,25,25)
-      radpos:
-      [[[ 13.          12.        ]]
-       [[ 12.70710678  12.70710678]]
-       [[ 12.          13.        ]]
-       [[ 11.29289322  12.70710678]]
-       [[ 11.          12.        ]]
-       [[ 11.29289322  11.29289322]]
-       [[ 12.          11.        ]]
-       [[ 12.70710678  11.29289322]]]
-      SQradpos:
-      {1.0: [[13, 12], [11, 12], [12, 13], [12, 11]], 1.41421356237: [[13, 13], [11, 13], [11, 11], [13, 11]]}
-    """
-    fname = 'radii_points'
-    radiipos = np.zeros((Nang, Lrad, 2), dtype = np.float)
-    SQradiipos = {}
-# Getting the range of radii:
-    pairsSQradii = read_SQradii(Npt=Lrad)
-    SQradii = pairsSQradii.keys()
-    Nradii = len(SQradii)
-# Angle loop
-    for ia in range(Nang):
-        angle = 2.*np.pi*ia/Nang
-        posangle = np.zeros((Lrad,2), dtype=np.float)
-# Points at that given angle
-        pangle = radial_points(angle,Lrad)
-        posangle[:,0] = pangle[:,0] + xpos
-        posangle[:,1] = pangle[:,1] + ypos
-        for ip in range(Lrad):
-            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 < dx-1 and jjpos >= 0 and jjpos < dy-1:
-                radiipos[ia,ip,0] = posangle[ip,0]
-                radiipos[ia,ip,1] = posangle[ip,1]
-# SQ radii values
-# Radii loop (avoiding 0)
-    pairswithin = {}
-    for ir in range(1,Nradii):
-        pairs = pairsSQradii[SQradii[ir]]
-# pairs loop
-        Npairs = len(pairs.flatten())/2
-        pairsw = []
-        for ip in range(Npairs):
-            if Npairs == 0:
-                break
-            else:
-                allpairs = grid_combinations(pairs[ip,0],pairs[ip,1])
-            for iap in range(len(allpairs)):
-                ppos = allpairs[iap]
-                iipos = xpos + ppos[0]
-                jjpos = ypos + ppos[1]
-# Creation of a matrix with all the grid points at each time-step
-                if iipos >= 0 and iipos < dx and jjpos >= 0 and jjpos < dy:
-                    pairsw.append([iipos,jjpos])
-        SQradiipos[SQradii[ir]] = pairsw
-    return radiipos, SQradiipos
-def hor_weight_int(SWval, SEval, NWval, NEval, xpos, ypos):
-    """ Function to weighted by distance interpolate a horizontal value from it closest 4 neighbourgs
-      field= a 4 points representing the environment of a given point
-      xpos, ypos: position to which one want to interpolate (relative to 0,0 at the corner SW)
-    >>> hor_weight_int(1., 1., 1., 2., 0.75, 0.75)
-.44888128193
-    """
-    fname = 'hor_weight_int'
-    vals = np.array([SWval, SEval, NWval, NEval])
-    print vals
-    if xpos > 1.:
-        print errormsg
-        print '  ' + fname + ': Wrong position to interpolate! (',xpos,',',ypos,     \
-          ') must be within (1,1) !!'
-        quit(-1)
-    vertexs = np.zeros((4,2), dtype=np.float)
-    vertexs[0,0] = 0.
-    vertexs[0,1] = 0.
-    vertexs[1,0] = 1.
-    vertexs[1,1] = 0.
-    vertexs[2,0] = 0.
-    vertexs[2,1] = 1.
-    vertexs[3,0] = 1.
-    vertexs[3,1] = 1.
-    print vertexs
-    dist = np.sqrt((vertexs[:,0]-xpos)**2 + (vertexs[:,1]-ypos)**2)
-    print dist
-    done = False
-    for id in range(4):
-        if dist[id] == 0.:
-            intval = vals[id]
-            done = True
-    if not done: intval = np.sum(vals/dist)/np.sum(1./dist)
-    return intval
 def compute_tevolboxtraj_radialsec(values, ncfile, varn):
     """ Function to compute tevolboxtraj_radialsec: temporal evolution at a given point along
 …
 #  '/bdd/PCER/workspace/lfita/etudes/FF/medic051215_2km/run/wrfout/wrfout_d02_2005-12-13_00:00:00', 'QVAPOR')
 #quit()
 def DatesFiles(values, ncfile, varn):
 …
 #fileobj = NetCDFFile('/home/lluis/PY/ERAI_pl199501_130-133.nc', 'r')
 #TimeSplitmat(fileobj, 'time', 'time', 'all', 'd,H')
-def rmNOnum(string):
-    """ Removing from a string all that characters which are not numbers
-    # From: http://stackoverflow.com/questions/4289331/python-extract-numbers-from-a-string
-    """
-    fname = 'rmNOnum'
-    newstr = str(re.findall("[-+]?\d+[\.]?\d*", string)[0])
-    return newstr
-def values_fortran_fmt(lin,fFMT):
-    """ Function to give back the values of an ASCII line according to its fortran printing format
-      lin= ASCII line
-      fFMT= list with the fortran FORMAT formats
-    forline = 'Natchitoches (RGNL)        1 11 0011  ( 31.733, -93.100) (  28, 157) ( 31.761, -93.113)   41.2 meters'
-    formats = ['A26', 'I2', 'I3', 'A6', 'A2', 'F7.3', 'A1', 'F8.3', 'A3', 'I4', 'A1', 'I4',
-      'A3', 'F7.3', 'A1', 'F8.3', 'A2', 'F6.1', 'A7']
-    >>> values_fortran_fmt(forline, formats)
-    ['Natchitoches (RGNL)        ', 1, 11, ' 0011  ', ' ( ', 31.733, ', ', -93.1, ') ( ', 28, ', ', 157, ')
-      ( ', 31.761, ', ', -93.113, ')  ', 41.2, ' meters']
-    """
-    fname = 'values_fortran_fmt'
-    afmts = ['A', 'D', 'F', 'I']
-    if lin == 'h':
-        print fname + '_____________________________________________________________'
-        print values_fortran_fmt.__doc__
-        quit()
-    fvalues = []
-    ichar=0
-    ival = 0
-    for ft in fFMT:
-        Lft = len(ft)
-        if ft[0:1] == 'A' or ft[0:1] == 'a':
-            Lfmt = int(ft[1:Lft+1])
-            fvalues.append(lin[ichar:ichar+Lfmt+1])
-            ichar = ichar + Lfmt
-        elif ft[0:1] == 'F' or ft[0:1] == 'f':
-            if ft.find('.') != -1:
-                Lft = len(ft.split('.')[0])
-            Lfmt = int(ft[1:Lft])
-            fvalues.append(np.float(rmNOnum(lin[ichar:ichar+Lfmt+1])))
-            ichar = ichar + Lfmt
-        elif ft[0:1] == 'D' or ft[0:1] == 'd':
-            if ft.find('.') != -1:
-                Lft = len(ft.split('.')[0])
-            Lfmt = int(ft[1:Lft])
-            fvalues.append(np.float64(rmNOnum(lin[ichar:ichar+Lfmt+1])))
-            ichar = ichar + Lfmt
-        elif ft[0:1] == 'I' or ft[0:1] == 'i':
-            Lfmt = int(ft[1:Lft+1])
-            fvalues.append(int(rmNOnum(lin[ichar:ichar+Lfmt+1])))
-            ichar = ichar + Lfmt
-        elif ft.find('X') != -1 or ft.find('x') != -1:
-            print '    ' + fname + ': skipping space'
-            ichar = ichar + int(rmNOnum(ft))
-        else:
-            print errormsg
-            print '  ' + fname + ": format '" + ft[0:1] + "' not ready!!"
-            print '    Available formats:',afmts
-            quit(-1)
-        ival = ival + 1
-    return fvalues
-def reduce_last_spaces(string):
-    """ Function to reduce the last right spaces from a string
-      string= string to remove the spaces at the righ side of the string
-    >>> reduce_last_spaces('Hola Lluis      ')
-    'Hola Lluis'
-    """
-    fname = 'reduce_last_spaces'
-    Lstring = len(string)
-    firstspace = -1
-    for istr in range(Lstring-1,0,-1):
-        if string[istr:istr+1] == ' ':
-            firstspace = istr
-        else:
-            break
-    if firstspace == -1:
-        print errormsg
-        print '  ' + fname + ": string '" + string + "' is not ended by spaces!!"
-        print '    char. number of last character:',ord(string[Lstring:Lstring+1])
-        quit(-1)
-    newstr = string[0:firstspace]
-    return newstr
-def squared_radial(Npt):
-    """ Function to provide the series of radii as composite of pairs (x,y) of gid cells
-      Npt= largest amount of grid points on x and y directions
-    >>> squared_radial(5)
-    [[0.0, 0, 0], [1.0, 1, 0], [1.4142135623730951, 1, 1], [2.0, 2, 0], [2.2360679774997898, 2, 1], [2.8284271247461903, 2, 2], [3.0, 3, 0], [3.1622776601683795, 3, 1], [3.6055512754639891, 3, 2], [4.0, 4, 0], [4.1231056256176606, 4, 1], [4.2426406871192848, 3, 3], [4.4721359549995796, 4, 2], [5.0, array([4, 3]), array([5, 0])], [5.0990195135927845, 5, 1], [5.3851648071345037, 5, 2], [5.6568542494923806, 4, 4], [5.8309518948453007, 5, 3], [6.4031242374328485, 5, 4], [7.0710678118654755, 5, 5]]
-    """
-    fname = 'squared_radial'
-    radii = []
-    gridradii = {}
-    singradii = []
-    Nradii = {}
-# Computing all radii
-##
-    for i in range(Npt+1):
-        if np.mod(i,100) == 0: print 'done: ',i
-        for j in range(i+1):
-#            ir = np.float(i)
-#            jr = np.float(j)
-            ir = np.float64(i)
-            jr = np.float64(j)
-            rad = np.sqrt(ir*ir + jr*jr)
-            if not searchInlist(singradii, rad):
-                singradii.append(rad)
-                gridradii[rad] = np.array([i,j], dtype=int)
-                Nradii[rad] = 1
-            else:
-#                print warnmsg
-#                print '  ' + fname + ': repeated radii',rad,'!!'
-#                print '    Previous values:', gridradii[rad]
-                oldvalues = gridradii[rad]
-                Nradii[rad] = Nradii[rad] + 1
-                newvalues = np.zeros((Nradii[rad],2), dtype=int)
-                if Nradii[rad] == 2:
-                    newvalues[0,:] = oldvalues
-                    Nstart = 1
-                else:
-                    Nstart = 0
-                for Nr in range(Nstart,Nradii[rad]-1):
-                    newvalues[Nr,:] = oldvalues[Nr,:]
-                newvalues[Nradii[rad]-1,:] = np.array([i,j], dtype=int)
-                gridradii[rad] = newvalues
-# Sorting values
-##
-## Only required to downwrite the output file
-#    radiif = open('SQradii.dat', 'w')
-#    radii = []
-#    prevrad = 0.
-#    for rad in sorted(gridradii.keys()):
-#        dradii = rad - prevrad
-## Check fo repeated values
-#        radii.append([rad, gridradii[rad][0], gridradii[rad][1]])
-#        radiif.write(str(rad) + ' ' + str(dradii) + ' ' + str(Nradii[rad]) + ' [ '+  \
-#          numVector_String(gridradii[rad].flatten(), ' ') + ' ] \n')
-#        prevrad = rad.copy()
-#    radiif.close()
-    return gridradii
 def getvalues_lonlat(values, ncfile):
 …
 #SpatialWeightedMean('reglonlat,XLONG,XLAT,west_east,south_north','wrfout_d01_2001-11-11_00:00:00','HGT')
-def fillvalue_kind(vartype, fillval0):
-    """ Function to provide the fiven fill_Value for a kind of variable
-      [vartype]= type of the variable
-      [fillval0]= value to take as fill value, 'std' for the standard value
-           Float = 1.e20
-           Character = '-'
-           Integer = -99999
-           Integer16 = -99999
-           Float64 = 1.e20
-           Integer32 = -99999
-    """
-    fname = 'fillvalue_kind'
-    if vartype == type('c'):
-        if fillval0 == 'std':
-            fillval = fillvalueC
-        else:
-            fillval = fillval0
-    elif vartype == type(int(1)):
-        if fillval0 == 'std':
-            fillval = fillValueI
-        else:
-            fillval = int(fillval0)
-    elif vartype == type(np.int16(1)):
-        if fillval0 == 'std':
-            fillval = np.int(9999999)
-        else:
-            fillval = np.int16(fillval0)
-    elif vartype == type(np.int32(1)):
-        if fillval0 == 'std':
-            fillval = np.int32(fillValueI)
-        else:
-            fillval = np.int32(fillval0)
-    elif vartype == type(np.float(1.)):
-        if fillval0 == 'std':
-            fillval = np.float(fillValueF)
-        else:
-            fillval = np.float(fillval0)
-    elif vartype == type(np.float32(1.)):
-        if fillval0 == 'std':
-            fillval = np.float32(fillValueF)
-        else:
-            fillval = np.float32(fillval0)
-    elif vartype == type(np.float64(1.)):
-        if fillval0 == 'std':
-            fillval = np.float64(fillValueF)
-        else:
-            fillval = np.float64(fillval0)
-    else:
-        print errormsg
-        print '  ' + fname + ': variable type', vartype, 'not ready !!'
-        quit(-1)
-    return fillval
 def nctype(vartype):
     """ Function to provide the string for the variable creation in a netCDF file
 …
     return ncs
-def varval_ind(var, ind):
-    """ Function to provide a value from a variable providing a value for each index
-      var= variable
-      ind= list of indices for each dimension of the variable
-    """
-    fname = 'varval_ind'
-    if len(var.shape) != len(ind):
-        print errormsg
-        print '  ' + fname + ': different shape:', var.shape,'and indices:', ind, '!!'
-        quit(-1)
-    varslice=[]
-    for iv in ind:
-        varslice.append(iv)
-    val = varval_ind[slice(varslice)]
-    return val
 def VarVal_FillValue(values, filen, varn):
     """ Function to transform a given value from a given variable to _FillValue in a netCDF file
 …
 #lonlatProj(100000.,100000.,'lonlat',True)
 #quit()
-def CoarselonlatFind(ilon, ilat, lonv, latv, per):
-    """ Function to search a given value from a coarser version of the data
-      ilon, ilat= original 2D matrices with the longitudes and the latitudes
-      lonv, latv= longitude and latitude to find
-      per= period (as fraction over 1) of the fractions of the original grid to use to explore
-      >>> npt=100
-      >>> lonval0 = np.arange(0.,npt*1.,1.)
-      >>> latval0 = np.arange(0.,npt*1.,1.)
-      >>> lonval = np.zeros((npt,npt), dtype=np.float)
-      >>> latval = np.zeros((npt,npt), dtype=np.float)
-      >>> for i in range(npt):
-      >>>     lonval[:,i] = lonval0[i]
-      >>>     latval[i,:] = latval0[i]
-      >>> CoarselonlatFind(lonval, latval, 5.25, 5.25, .1)
-      (array([5, 5]), 0.35355339059327379)
-    """
-    fname = 'CoarselonlatFind'
-    dx = ilon.shape[1]
-    dy = ilon.shape[0]
-    nlon = np.min(ilon)
-    xlon = np.max(ilon)
-    nlat = np.min(ilat)
-    xlat = np.max(ilat)
-    if lonv < nlon or lonv > xlon:
-        print errormsg
-        print '  ' + fname + ': longitude outside data range!!'
-        print '    given value:', lonv,'outside [',nlon,',',xlon,']'
-        quit(-1)
-    if latv < nlat or latv > xlat:
-        print errormsg
-        print '  ' + fname + ': latitude outside data range!!'
-        print '    given value:', latv,'outside [',nlat,',',xlat,']'
-        quit(-1)
-    fracx = int(dx*per)
-    fracy = int(dy*per)
-    fraclon = ilon[::fracy,::fracx]
-    fraclat = ilat[::fracy,::fracx]
-    fracdx = fraclon.shape[1]
-    fracdy = fraclon.shape[0]
-#    print 'fraclon _______'
-#    print fraclon
-#    print 'fraclat _______'
-#    print fraclat
-# Fraction point
-    difffraclonlat = np.sqrt((fraclon-lonv)**2. + (fraclat-latv)**2.)
-    mindifffracLl = np.min(difffraclonlat)
-    ilatlonfrac = index_mat(difffraclonlat, mindifffracLl)
-#    print 'mindifffracLl:', mindifffracLl, 'ilatlonfrac:', ilatlonfrac
-#    print 'frac lon, lat:', fraclon[ilatlonfrac[0],ilatlonfrac[1]], ',',             \
-#      fraclat[ilatlonfrac[0],ilatlonfrac[1]]
-#    print 'values lon, lat:', lonv, latv
-# Providing fraction range
-    if fraclon[ilatlonfrac[0],ilatlonfrac[1]] >= lonv and                            \
-      fraclat[ilatlonfrac[0],ilatlonfrac[1]] >= latv:
-#        ifracbeg = [ilatlonfrac[0]-1, ilatlonfrac[1]-1]
-#        ifracend = [ilatlonfrac[0], ilatlonfrac[1]]
-        if ilatlonfrac[0] > 0:
-            iybeg = (ilatlonfrac[0]-1)*fracy
-            iyend = ilatlonfrac[0]*fracy+1
-        else:
-            iybeg = 0
-            iyend = fracy+1
-        if ilatlonfrac[1] > 0:
-            ixbeg = (ilatlonfrac[1]-1)*fracx
-            ixend = ilatlonfrac[1]*fracx+1
-        else:
-            ixbeg = 0
-            ixend = fracx+1
-    elif fraclon[ilatlonfrac[0],ilatlonfrac[1]] < lonv and                           \
-      fraclat[ilatlonfrac[0],ilatlonfrac[1]] >= latv:
-#        ifracbeg = [ilatlonfrac[0]-1, ilatlonfrac[1]]
-#        ifracend = [ilatlonfrac[0], ilatlonfrac[1]+1]
-        if ilatlonfrac[0] > 0:
-            iybeg = (ilatlonfrac[0]-1)*fracy
-            iyend = ilatlonfrac[0]*fracy+1
-        else:
-            iybeg = 0
-            iyend = fracy+1
-        if ilatlonfrac[1] < fracdx:
-            if ilatlonfrac[1] != 0:
-                ixbeg = (ilatlonfrac[1]-1)*fracx
-                ixend = ilatlonfrac[1]*fracx+1
-            else:
-                ixbeg = 0
-                ixend = fracx+1
-        else:
-            ixbeg = fracdx*fracx
-            ixend = dx+1
-    elif fraclon[ilatlonfrac[0],ilatlonfrac[1]] < lonv and                           \
-      fraclat[ilatlonfrac[0],ilatlonfrac[1]] < latv:
-#        ifracbeg = [ilatlonfrac[0], ilatlonfrac[1]]
-#        ifracend = [ilatlonfrac[0]+1, ilatlonfrac[1]+1]
-        if ilatlonfrac[0] < fracdy:
-            if ilatlonfrac[0] != 0:
-                iybeg = (ilatlonfrac[0]-1)*fracy
-                iyend = ilatlonfrac[0]*fracy+1
-            else:
-                iybeg = 0
-                iyend = fracy+1
-        else:
-            iybeg = fracdy*fracy
-            iyend = dy+1
-        if ilatlonfrac[1] < fracdx:
-            if ilatlonfrac[1] != 0:
-                ixbeg = (ilatlonfrac[1]-1)*fracx
-                ixend = ilatlonfrac[1]*fracx+1
-            else:
-                ixbeg = 0
-                ixend = fracx+1
-        else:
-            ixbeg = fracdx*fracx
-            ixend = dx+1
-    elif fraclon[ilatlonfrac[0],ilatlonfrac[1]] >= lonv and                          \
-      fraclat[ilatlonfrac[0],ilatlonfrac[1]] < latv:
-#        ifracbeg = [ilatlonfrac[0], ilatlonfrac[1]-1]
-#        ifracend = [ilatlonfrac[0]+1, ilatlonfrac[1]]
-        if ilatlonfrac[0] < fracdy:
-            if ilatlonfrac[0] != 0:
-                iybeg = (ilatlonfrac[0]-1)*fracy
-                iyend = ilatlonfrac[0]*fracy+1
-            else:
-                iybeg = 0
-                iyend = fracy+1
-        else:
-            iybeg = fracdy*fracy
-            iyend = dy+1
-        if ilatlonfrac[1] > 0:
-            ixbeg = (ilatlonfrac[1]-1)*fracx
-            ixend = ilatlonfrac[1]*fracx+1
-        else:
-            ixbeg = 0
-            ixend = fracx+1
-#    ibeg = [iybeg, ixbeg]
-#    iend = [iyend, ixend]
-#    print 'find window; beginning', ibeg, 'end:', iend
-    lon = ilon[iybeg:iyend,ixbeg:ixend]
-    lat = ilat[iybeg:iyend,ixbeg:ixend]
-#    print 'lon _______'
-#    print lon
-#    print 'lat _______'
-#    print lat
-# Find point
-    difflonlat = np.sqrt((lon-lonv)**2. + (lat-latv)**2.)
-    mindiffLl = np.min(difflonlat)
-    ilatlon = index_mat(difflonlat, mindiffLl)
-    ilatlon[0] = ilatlon[0] + iybeg
-    ilatlon[1] = ilatlon[1] + ixbeg
-#    print 'mindiffLl:', mindiffLl, 'ilatlon:', ilatlon
-#    print 'lon, lat:', lon[ilatlon[0],ilatlon[1]], ',', lat[ilatlon[0],ilatlon[1]]
-#    quit()
-    return ilatlon, mindiffLl
-def CoarselonlatFindAll(onc, ilon, ilat, longvs, latvs, per, frac, out):
-    """ Function to search all values from a coarser version of the data
-      onc= netCDF object from an existing file
-      ilon, ilat= original 2D matrices with the longitudes and the latitudes
-      lonv, latv= longitude and latitude to find
-      per= period (as fraction over 1) of the fractions of the original grid to use to explore
-      frac= frequency of points at which syncronize file with calculations
-      out= True/False if outcome is required
-      >>> npt=100
-      >>> lonval0 = np.arange(0.,npt*1.,1.)
-      >>> latval0 = np.arange(0.,npt*1.,1.)
-      >>> lonval = np.zeros((npt,npt), dtype=np.float)
-      >>> latval = np.zeros((npt,npt), dtype=np.float)
-      >>> lonvalues = np.arange(0.25,10.,0.25)
-      >>> latvalues = np.arange(0.25,10.,0.25)
-      >>> for i in range(npt):
-      >>>     lonval[:,i] = lonval0[i]
-      >>>     latval[i,:] = latval0[i]
-      >>> CoarselonlatFindAll(onewnc, lonval, latval, lovalues, latvalues, .1, 100)
-      (array([5, 5]), 0.35355339059327379)
-    """
-    fname = 'CoarselonlatFindAll'
-    dx = ilon.shape[1]
-    dy = ilon.shape[0]
-    nlon = np.min(ilon)
-    xlon = np.max(ilon)
-    nlat = np.min(ilat)
-    xlat = np.max(ilat)
-    Nallpts = len(longvs)
-    fracx = int(dx*per)
-    fracy = int(dy*per)
-    fraclon = ilon[::fracy,::fracx]
-    fraclat = ilat[::fracy,::fracx]
-    fracdx = fraclon.shape[1]
-    fracdy = fraclon.shape[0]
-#    print 'fraclon _______'
-#    print fraclon
-#    print 'fraclat _______'
-#    print fraclat
-    if out:
-        ilatlon = np.zeros((2,Nallpts), dtype=int)
-        mindiffLl = np.zeros((Nallpts), dtype=np.float)
-    iptpos=0
-    for iv in range(Nallpts):
-        lonv = longvs[iv]
-        latv = latvs[iv]
-        if lonv < nlon or lonv > xlon:
-            print errormsg
-            print '  ' + fname + ': longitude outside data range!!'
-            print '    given value:', lonv,'outside [',nlon,',',xlon,']'
-            quit(-1)
-        if latv < nlat or latv > xlat:
-            print errormsg
-            print '  ' + fname + ': latitude outside data range!!'
-            print '    given value:', latv,'outside [',nlat,',',xlat,']'
-            quit(-1)
-# Fraction point
-        difffraclonlat = np.sqrt((fraclon-lonv)**2. + (fraclat-latv)**2.)
-        mindifffracLl = np.min(difffraclonlat)
-        ilatlonfrac = index_mat(difffraclonlat, mindifffracLl)
-#        print 'values lon, lat:', lonv, latv
-#        print 'mindifffracLl:', mindifffracLl, 'ilatlonfrac:', ilatlonfrac
-#        print 'frac lon, lat:', fraclon[ilatlonfrac[0],ilatlonfrac[1]], ',',         \
-#          fraclat[ilatlonfrac[0],ilatlonfrac[1]]
-# Providing fraction range
-        if fraclon[ilatlonfrac[0],ilatlonfrac[1]] >= lonv and                        \
-          fraclat[ilatlonfrac[0],ilatlonfrac[1]] >= latv:
-#            ifracbeg = [ilatlonfrac[0]-1, ilatlonfrac[1]-1]
-#            ifracend = [ilatlonfrac[0], ilatlonfrac[1]]
-            if ilatlonfrac[0] > 0:
-                iybeg = (ilatlonfrac[0]-1)*fracy
-                iyend = ilatlonfrac[0]*fracy+1
-            else:
-                iybeg = 0
-                iyend = fracy+1
-            if ilatlonfrac[1] > 0:
-                ixbeg = (ilatlonfrac[1]-1)*fracx
-                ixend = ilatlonfrac[1]*fracx+1
-            else:
-                ixbeg = 0
-                ixend = fracx+1
-        elif fraclon[ilatlonfrac[0],ilatlonfrac[1]] < lonv and                       \
-          fraclat[ilatlonfrac[0],ilatlonfrac[1]] >= latv:
-#            ifracbeg = [ilatlonfrac[0]-1, ilatlonfrac[1]]
-#            ifracend = [ilatlonfrac[0], ilatlonfrac[1]+1]
-            if ilatlonfrac[0] > 0:
-                iybeg = (ilatlonfrac[0]-1)*fracy
-                iyend = ilatlonfrac[0]*fracy+1
-            else:
-                iybeg = 0
-                iyend = fracy+1
-            if ilatlonfrac[1] < fracdx:
-                if ilatlonfrac[1] != 0:
-                    ixbeg = (ilatlonfrac[1]-1)*fracx
-                    ixend = ilatlonfrac[1]*fracx+1
-                else:
-                    ixbeg = 0
-                    ixend = fracx+1
-            else:
-                ixbeg = fracdx*fracx
-                ixend = dx+1
-        elif fraclon[ilatlonfrac[0],ilatlonfrac[1]] < lonv and                       \
-          fraclat[ilatlonfrac[0],ilatlonfrac[1]] < latv:
-#            ifracbeg = [ilatlonfrac[0], ilatlonfrac[1]]
-#            ifracend = [ilatlonfrac[0]+1, ilatlonfrac[1]+1]
-            if ilatlonfrac[0] < fracdy:
-                if ilatlonfrac[0] != 0:
-                    iybeg = (ilatlonfrac[0]-1)*fracy
-                    iyend = ilatlonfrac[0]*fracy+1
-                else:
-                    iybeg = 0
-                    iyend = fracy+1
-            else:
-                iybeg = fracdy*fracy
-                iyend = dy+1
-            if ilatlonfrac[1] < fracdx:
-                if ilatlonfrac[1] != 0:
-                    ixbeg = (ilatlonfrac[1]-1)*fracx
-                    ixend = ilatlonfrac[1]*fracx+1
-                else:
-                    ixbeg = 0
-                    ixend = fracx+1
-            else:
-                ixbeg = fracdx*fracx
-                ixend = dx+1
-        elif fraclon[ilatlonfrac[0],ilatlonfrac[1]] >= lonv and                      \
-          fraclat[ilatlonfrac[0],ilatlonfrac[1]] < latv:
-#            ifracbeg = [ilatlonfrac[0], ilatlonfrac[1]-1]
-#            ifracend = [ilatlonfrac[0]+1, ilatlonfrac[1]]
-            if ilatlonfrac[0] < fracdy:
-                if ilatlonfrac[0] != 0:
-                    iybeg = (ilatlonfrac[0]-1)*fracy
-                    iyend = ilatlonfrac[0]*fracy+1
-                else:
-                    iybeg = 0
-                    iyend = fracy+1
-            else:
-                iybeg = fracdy*fracy
-                iyend = dy+1
-            if ilatlonfrac[1] > 0:
-                ixbeg = (ilatlonfrac[1]-1)*fracx
-                ixend = ilatlonfrac[1]*fracx+1
-            else:
-                ixbeg = 0
-                ixend = fracx+1
-#        ibeg = [iybeg, ixbeg]
-#        iend = [iyend, ixend]
-#        print 'find window; beginning', ibeg, 'end:', iend
-        lon = ilon[iybeg:iyend,ixbeg:ixend]
-        lat = ilat[iybeg:iyend,ixbeg:ixend]
-#        print 'Looking within _______'
-#        print 'lon:',lon[0,0],',',lon[-1,-1]
-#        print 'lat:',lat[0,0],',',lat[-1,-1]
-#        print 'lon _______'
-#        print lon
-#        print 'lat _______'
-#        print lat
-# Find point
-        difflonlat = np.sqrt((lon-lonv)**2. + (lat-latv)**2.)
-        if out:
-            mindiffLl[iv] = np.min(difflonlat)
-            ilatlon[:,iv] = index_mat(difflonlat, mindiffLl)
-            ilatlon[0,iv] = ilatlon[0,iv] + iybeg
-            ilatlon[1,iv] = ilatlon[1,iv] + ixbeg
-        else:
-            mindiffLl = np.min(difflonlat)
-            ilatlon = index_mat(difflonlat, mindiffLl)
-            ilatlon[0] = ilatlon[0] + iybeg
-            ilatlon[1] = ilatlon[1] + ixbeg
-#        print 'mindiffLl:', mindiffLl, 'ilatlon:', ilatlon
-#        print 'lon, lat:', lon[ilatlon[0],ilatlon[1]], ',', lat[ilatlon[0],ilatlon[1]]
-#        quit()
-#        if mindiffLl == 0. and type(newvar[ilatlon[0],ilatlon[1]]) == type(amsk):
-#            percendone(iv,Ninpts,0.5,'done:')
-        if mindiffLl == 0.:
-            iptpos = iptpos + 1
-# Speeding it up!
-#            if mindiffLl > mindiff:
-#                print errormsg
-#                print '  ' + fname + ': for point #', iv,'lon,lat in ' +             \
-#                  'incomplet map:', lonvs[iv], ',', latvs[iv], 'there is ' +         \
-#                  'not a set of lon,lat in the completed map closer than: ',         \
-#                  mindiff, '!!'
-#                print '    minimum difference:', mindiffLl
-#                quit(-1)
-# Speeding it up!
-#            if ilatlon[0] >= 0 and ilatlon[1] >= 0:
-#                newvar[ilatlon[0],ilatlon[1]] = ovar[iv]
-#                newvarin[ilatlon[0],ilatlon[1]] = iv
-#                newvarinpt[iv] = 1
-#                newvarindiff[iv] = mindiffLl
-            newvar[ilatlon[0],ilatlon[1]] = ovar[iv]
-            newvarin[ilatlon[0],ilatlon[1]] = iv
-            newvarinpt[iv] = 1
-            newvarindiff[iv] = mindiffLl
-#                print 'Lluis iv:', newvarin[ilatlon[0],ilatlon[1]],                  \
-#                  'localized:', newvarinpt[iv], 'values:',                           \
-#                  newvar[ilatlon[0],ilatlon[1]], 'invalues:', ovar[iv],              \
-#                  'mindist:', newvarindiff[iv], 'point:',ilatlon
-#            else:
-#                print errormsg
-#                print '  ' + fname + ': point iv:', iv, 'at', lonvs[iv], ',',        \
-#                  latvs[iv],' not relocated !!'
-#                print '    mindiffl:', mindiffLl, 'ilon:', ilatlon[1],               \
-#                  'ilatlon:', ilatlon[1]
-#                quit(-1)
-            if np.mod(iptpos,fracs) == 0:
-                print 'Lluis syncronizing!'
-                onc.sync()
-    if out:
-        return ilatlon, mindiffLl
-    else:
-        return
-def lonlatFind(ilon, ilat, lonv, latv):
-    """ Function to search a given value from a lon, lat 2D data
-      ilon, ilat= original 2D matrices with the longitudes and the latitudes
-      lonv, latv= longitude and latitude to find
-      >>> npt=100
-      >>> lonval0 = np.arange(0.,npt*1.,1.)
-      >>> latval0 = np.arange(0.,npt*1.,1.)
-      >>> lonval = np.zeros((npt,npt), dtype=np.float)
-      >>> latval = np.zeros((npt,npt), dtype=np.float)
-      >>> for i in range(npt):
-      >>>     lonval[:,i] = lonval0[i]
-      >>>     latval[i,:] = latval0[i]
-      >>> lonlatFind(lonval, latval, 5.25, 5.25)
-      (array([5, 5]), 0.35355339059327379)
-    """
-    fname = 'lonlatFind'
-    dx = ilon.shape[1]
-    dy = ilon.shape[0]
-    nlon = np.min(ilon)
-    xlon = np.max(ilon)
-    nlat = np.min(ilat)
-    xlat = np.max(ilat)
-    if lonv < nlon or lonv > xlon:
-        print errormsg
-        print '  ' + fname + ': longitude outside data range!!'
-        print '    given value:', lonv,'outside [',nlon,',',xlon,']'
-        quit(-1)
-    if latv < nlat or latv > xlat:
-        print errormsg
-        print '  ' + fname + ': latitude outside data range!!'
-        print '    given value:', latv,'outside [',nlat,',',xlat,']'
-        quit(-1)
-# Find point
-    difflonlat = np.sqrt((ilon-lonv)**2. + (ilat-latv)**2.)
-    mindiffLl = np.min(difflonlat)
-    ilatlon = index_mat(difflonlat, mindiffLl)
-#    print 'mindiffLl:', mindiffLl, 'ilatlon:', ilatlon
-#    print 'lon, lat:', lon[ilatlon[0],ilatlon[1]], ',', lat[ilatlon[0],ilatlon[1]]
-    return ilatlon, mindiffLl
 def Partialmap_Entiremap(values, filen, varn):

Note: See TracChangeset for help on using the changeset viewer.

Context Navigation

Changeset 747 in lmdz_wrf

Legend:

trunk/tools/nc_var_tools.py

Download in other formats: