Index: /trunk/tools/nc_var_tools.py
===================================================================
--- /trunk/tools/nc_var_tools.py	(revision 270)
+++ /trunk/tools/nc_var_tools.py	(revision 271)
@@ -9763,5 +9763,5 @@
           [time] [xpoint] [ypoint]
         [Tbeg]: equivalent first time-step of the trajectory within the netCDF file
-        [lonname],[latname],[timename]: longitude, latitude and time variables names
+        [lonname],[latname],[zname],[timename]: longitude, latitude, z and time variables names
         [timekind]: kind of time
           'cf': cf-compilant
@@ -9770,5 +9770,5 @@
         [circler]: radius in grid points of a centerd circle
       ncfile= netCDF file to use
-      varn= variable name
+      varn= ',' list of variables' name ('all', for all variables)
         EMERGENCY version, assuming 3D [time,lat,lon] variable !
     """
@@ -9788,5 +9788,6 @@
     lonn = values.split(',')[1]
     latn = values.split(',')[2]
-    timn = values.split(',')[3]
+    zn = values.split(',')[3]
+    timn = values.split(',')[4]
     timekind = values.split(',')[4]
     boxs = int(values.split(',')[5])
@@ -9827,4 +9828,12 @@
     timobj = objfile.variables[timn]
 
+    if varn.find(',') != -1:
+        varns = varn.split(',')
+    else:
+        if varn == 'all':
+            varns = objfile.variables
+        else:
+            varns = varn
+
     if len(lonobj.shape) == 3:
         lonv = lonobj[0,:,:]
@@ -9844,6 +9853,4 @@
           ncfile + '" !!!!!'
         quit(-1)
-
-    varobj = objfile.variables[varn]
 
 # Selecting accordingly a trajectory
@@ -9860,22 +9867,27 @@
     trajobj = open(trajfile,'r')
 
+# Trajectory values/slices
+##
+    iline = 0
+    it = 0
+
+    xrangeslice = []
+    yrangeslice = []
+    xrangeslice2D = []
+    yrangeslice2D = []
+
+    cxrangeslice = []
+    cyrangeslice = []
+    cxrangeslice2D = []
+    cyrangeslice2D = []
+
+    cslicev = []
+    cslice2D = []
+    cslicevnoT = []
+
     gtrajvals = np.zeros((Ttraj,3), dtype=int)
     trajvals = np.zeros((Ttraj,3), dtype=np.float)
-    varvals = np.ones(tuple([Ttraj,boxs,boxs]), dtype=np.float)
-    lonvals = np.ones(tuple([Ttraj,boxs,boxs]), dtype=np.float)
-    latvals = np.ones(tuple([Ttraj,boxs,boxs]), dtype=np.float)
-    rvarvals = np.ones(tuple([Ttraj,Nrad*2+1,Nrad*2+1]), dtype=np.float)
-    rlonvals = np.ones(tuple([Ttraj,Nrad*2+1,Nrad*2+1]), dtype=np.float)
-    rlatvals = np.ones(tuple([Ttraj,Nrad*2+1,Nrad*2+1]), dtype=np.float)
-
-    statvarvals = np.ones(tuple([Ttraj,6]), dtype=np.float)
-    rstatvarvals = np.ones(tuple([Ttraj,6]), dtype=np.float)
-
-    iline = 0
-    it = 0
+
     for line in trajobj:
-        slicev = []
-        slice2D = []
-        slicevnoT = []
 
 ## Skipping first line
@@ -9950,74 +9962,15 @@
                     xend2D = boxs
 
-                slicev.append(slice(yinit, yend))
-                slicev.append(slice(xinit, xend))
-
-                slicevnoT.append(slice(yinit, yend))
-                slicevnoT.append(slice(xinit, xend))
-
-                slice2D.append(slice(yinit2D, yend2D))
-                slice2D.append(slice(xinit2D, xend2D))
-
-                varvalst[tuple(slice2D)] = varobj[tuple(slicev)]
-                varvals[it,:,:] = varvalst
-
-# box stats values
-                maskedvals = ma.masked_values (varvalst, fillValue)
-                statvarvals[it,0] = varvalst[box2,box2]
-                statvarvals[it,1] = maskedvals.min()
-                statvarvals[it,2] = maskedvals.max()
-                statvarvals[it,3] = maskedvals.mean()
-                maskedvals2 = maskedvals*maskedvals
-                statvarvals[it,4] = maskedvals2.mean()
-                statvarvals[it,5] = np.sqrt(statvarvals[it,4] -                      \
-                  statvarvals[it,3]*statvarvals[it,3])
-
-                varvalst[tuple(slice2D)] = lonv[tuple(slicevnoT)]
-                lonvals[it,:,:] = varvalst
-                varvalst[tuple(slice2D)] = latv[tuple(slicevnoT)]           
-                latvals[it,:,:] = varvalst
-
-            else:
-                slicev.append(slice(gtrajvals[it,2]-box2, gtrajvals[it,2]+box2 + 1))
-                slicev.append(slice(gtrajvals[it,1]-box2, gtrajvals[it,1]+box2 + 1))
-                slicevnoT.append(slice(gtrajvals[it,2]-box2, gtrajvals[it,2]+box2+1))
-                slicevnoT.append(slice(gtrajvals[it,1]-box2, gtrajvals[it,1]+box2+1))
-
-                slice2D.append(slice(gtrajvals[it,2]-box2, gtrajvals[it,2]+box2 + 1))
-                slice2D.append(slice(gtrajvals[it,1]-box2, gtrajvals[it,1]+box2 + 1))
-
-                varvalst = varobj[tuple(slicev)]
-# box values
-
-                varvals[it,:,:] = varvalst
-#            print 'values at time t______'
-#            print varvalst
-
-# box stats values
-                statvarvals[it,0] = varvalst[box2,box2]
-                statvarvals[it,1] = np.min(varvalst)
-                statvarvals[it,2] = np.max(varvalst)
-                statvarvals[it,3] = np.mean(varvalst)
-                statvarvals[it,4] = np.mean(varvalst*varvalst)
-                statvarvals[it,5] = np.sqrt(statvarvals[it,4] -                      \
-                  statvarvals[it,3]*statvarvals[it,3])
-
-                varvalst = lonv[tuple(slicevnoT)]
-                lonvals[it,:,:] = varvalst
-                varvalst = latv[tuple(slicevnoT)]           
-                latvals[it,:,:] = varvalst
+                yrangeslice.append([yinit, yend])
+                xrangeslice.append([xinit, xend])
+                yrangeslice2D.append([yinit2D, yend2D])
+                xrangeslice2D.append([xinit2D, xend2D])
 
 # circle values
-            slicev = []
-            slice2D = []
-            slicevnoT = []
-
-            slicev.append(gtrajvals[it,0])
+            cslicev.append(gtrajvals[it,0])
             circdist = radius_dist(dimy, dimx, gtrajvals[it,2], gtrajvals[it,1])
 
             if gtrajvals[it,2]-Nrad < 0 or gtrajvals[it,2]+Nrad + 1 > dimy + 1 \
               or gtrajvals[it,1]-Nrad < 0 or gtrajvals[it,1]+Nrad + 1 > dimx + 1:
-
-                rvarvalst = np.ones((Nrad*2+1, Nrad*2+1), dtype=np.float)*fillValue
 
                 if gtrajvals[it,2]-Nrad < 0: 
@@ -10049,69 +10002,8 @@
                     xend2D = 2*Nrad+1
 
-                slicev.append(slice(yinit, yend))
-                slicev.append(slice(xinit, xend))
-
-                slicevnoT.append(slice(yinit, yend))
-                slicevnoT.append(slice(xinit, xend))
-
-                slice2D.append(slice(yinit2D, yend2D))
-                slice2D.append(slice(xinit2D, xend2D))
-
-                rvarvalst[tuple(slice2D)] = varobj[tuple(slicev)]
-                rvarvalst[tuple(slice2D)] = np.where(circdist[tuple(slice2D)] >      \
-                  np.float(Nrad), fillValue, rvarvalst[tuple(slice2D)])
-
-                rvarvals[it,:,:] = rvarvalst
-
-# circle stats values
-                maskedvals = ma.masked_values (rvarvalst, fillValue)
-                rstatvarvals[it,0] = rvarvalst[Nrad,Nrad]
-                rstatvarvals[it,1] = maskedvals.min()
-                rstatvarvals[it,2] = maskedvals.max()
-                rstatvarvals[it,3] = maskedvals.mean()
-                maskedvals2 = maskedvals*maskedvals
-                rstatvarvals[it,4] = maskedvals2.mean()
-                rstatvarvals[it,5] = np.sqrt(rstatvarvals[it,4] -                    \
-                  rstatvarvals[it,3]*rstatvarvals[it,3])
-
-                rvarvalst[tuple(slice2D)] = lonv[tuple(slicevnoT)]
-                rlonvals[it,:,:] = rvarvalst
-                rvarvalst[tuple(slice2D)] = latv[tuple(slicevnoT)]           
-                rlatvals[it,:,:] = rvarvalst
-
-            else:
-                slicev.append(slice(gtrajvals[it,2]-Nrad, gtrajvals[it,2]+Nrad + 1))
-                slicev.append(slice(gtrajvals[it,1]-Nrad, gtrajvals[it,1]+Nrad + 1))
-                slicevnoT.append(slice(gtrajvals[it,2]-Nrad, gtrajvals[it,2]+Nrad+1))
-                slicevnoT.append(slice(gtrajvals[it,1]-Nrad, gtrajvals[it,1]+Nrad+1))
-
-                slice2D.append(slice(gtrajvals[it,2]-Nrad, gtrajvals[it,2]+Nrad + 1))
-                slice2D.append(slice(gtrajvals[it,1]-Nrad, gtrajvals[it,1]+Nrad + 1))
-
-                rvarvalst = varobj[tuple(slicev)]
-                cdist = circdist[tuple(slicevnoT)]
-# circle values
-                rvarvalst = np.where(cdist > np.float(Nrad), fillValue, rvarvalst)
-
-                rvarvals[it,:,:] = rvarvalst
-
-# circle stats values
-                maskedvals = ma.masked_values (rvarvalst, fillValue)
-                rstatvarvals[it,0] = rvarvalst[Nrad,Nrad]
-                rstatvarvals[it,1] = maskedvals.min()
-                rstatvarvals[it,2] = maskedvals.max()
-                rstatvarvals[it,3] = maskedvals.mean()
-                maskedvals2 = maskedvals*maskedvals
-                rstatvarvals[it,4] = maskedvals2.mean()
-                rstatvarvals[it,5] = np.sqrt(rstatvarvals[it,4] -                    \
-                  rstatvarvals[it,3]*rstatvarvals[it,3])
-
-                rvarvalst = lonv[tuple(slicevnoT)]
-                rlonvals[it,:,:] = rvarvalst
-                rvarvalst = latv[tuple(slicevnoT)]           
-                rlatvals[it,:,:] = rvarvalst
-
-            it = it + 1
-#            print 'statistics:',rstatvarvals[it,:] 
+                cyrangeslice.append([yinit, yend])
+                cxrangeslice.append([xinit, xend])
+                cyrangeslice2D.append([yinit2D, yend2D])
+                cxrangeslice2D.append([xinit2D, xend2D])
 
         iline = iline + 1
@@ -10153,81 +10045,504 @@
     newvar[:] = latvals
 
+    statnames = ['minbox', 'maxbox', 'meanbox', 'mean2box', 'stdevbox']
+    vstlname = ['minimum value within', 'maximum value within', 'mean value within', \
+      'squared mean value within', 'standard deviation value within']
+    cstatnames = ['mincircle','maxcircle','meancircle','mean2circle','stdevcircle']
+
+# Getting values
+##
+    ivar = 0
+    for vn in varns:
+        varobj = objfile.variables[vn]
+        Nvardims = len(varobj.shape)
+
+        slicev = []
+        slice2D = []
+        slicevnoT = []
+        cslicev = []
+        cslice2D = []
+        cslicevnoT = []
+
+        if Nvardims == 4:
+# Too optimistic, but at this stage...
+            dimt = varobj.shape[0]
+            dimz = varobj.shape[1]
+ 
+            if not objofile.dimensions.has_key('z'):
+                newdim = objofile.createDimension('z', dimt)
+
+                varzobj = objfile.variables[zn]
+                newvar = objofile.createVariable(zn, 'f8', ('z'))
+
+                if searchInlist(varzobj.ncattrs(),'standard_name'):
+                    vsname = varzobj.getncattr('standard_name')
+                else:
+                    vsname = variables_values(zn)[1]
+                if searchInlist(varzobj.ncattrs(),'long_name'):
+                    vlname = varzobj.getncattr('long_name')
+                else:
+                    vlname = variables_values(zn)[4].replace('|',' ')
+                if searchInlist(varzobj.ncattrs(),'units'):
+                    vunits = varzobj.getncattr('units')
+                else:
+                    vunits = variables_values(zn)[5].replace('|',' ')
+
+                newattr = basicvardef(newvar, vsname, vlname, vunits)
+                if len(varzobj.shape) == 1:
+                    newvar[:] = varzobj[:]
+                elif len(varzobj.shape) == 2:
+                    newvar[:] = varzobj[0,:]
+                elif len(varzobj.shape) == 3:
+                    newvar[:] = varzobj[0,0,:]
+
+            varvals = np.ones(tuple([Ttraj,dimz,boxs,boxs]), dtype=np.float)
+            lonvals = np.ones(tuple([Ttraj,dimz,boxs,boxs]), dtype=np.float)
+            latvals = np.ones(tuple([Ttraj,dimz,boxs,boxs]), dtype=np.float)
+            rvarvals = np.ones(tuple([Ttraj,dimz,Nrad*2+1,Nrad*2+1]), dtype=np.float)
+            rlonvals = np.ones(tuple([Ttraj,dimz,Nrad*2+1,Nrad*2+1]), dtype=np.float)
+            rlatvals = np.ones(tuple([Ttraj,dimz,Nrad*2+1,Nrad*2+1]), dtype=np.float)
+
+            statvarvals = np.ones(tuple([Ttraj,6]), dtype=np.float)
+            rstatvarvals = np.ones(tuple([Ttraj,6]), dtype=np.float)
+
+            for it in range(dimt):
+# box values
+                slicev.appen(slice(0:dimz))
+                slicev.appen(slice(yrangeslice[it]))
+                slicev.appen(slice(xrangeslice[it]))
+
+                slicevnoT.appen(slice(0:dimz))
+                slicevnoT.appen(slice(yrangeslice[it]))
+                slicevnoT.appen(slice(xrangeslice[it]))
+
+                slice2D.appen(slice(0:dimz))
+                slice2D.appen(slice(yrangeslice2D[it]))
+                slice2D.appen(slice(xrangeslice2D[it]))
+
+                if gtrajvals[it,2]-box2 < 0 or gtrajvals[it,2]+box2 + 1 > dimy + 1   \
+                  or gtrajvals[it,1]-box2 < 0 or gtrajvals[it,1]+box2 + 1 > dimx + 1:
+
+                    rvarvalst = np.ones((dimz,Nrad*2+1,Nrad*2+1),dtype=np.float)*    \
+                      fillValue
+
+                    varvalst[tuple(slice2D)] = varobj[tuple(slicev)]
+                    varvals[it,:,:,:] = varvalst
+    
+# box stats values
+                    maskedvals = ma.masked_values (varvalst, fillValue)
+                    statvarvals[it,:,0] = varvalst[box2,box2]
+                    statvarvals[it,:,1] = maskedvals.min()
+                    statvarvals[it,:,2] = maskedvals.max()
+                    statvarvals[it,:,3] = maskedvals.mean()
+                    maskedvals2 = maskedvals*maskedvals
+                    statvarvals[it,:,4] = maskedvals2.mean()
+                    statvarvals[it,:,5] = np.sqrt(statvarvals[it,:,4] -              \
+                      statvarvals[it,:,3]*statvarvals[it,:,3])
+
+                    varvalst[tuple(slice2D)] = lonv[tuple(slicevnoT)]
+                    lonvals[it,:,:,:] = varvalst
+                    varvalst[tuple(slice2D)] = latv[tuple(slicevnoT)]           
+                    latvals[it,:,:,:] = varvalst
+
+                else:
+                    slicev.append(slice(gtrajvals[it,2]-box2, gtrajvals[it,2]+box2+1))
+                    slicev.append(slice(gtrajvals[it,1]-box2, gtrajvals[it,1]+box2+1))
+                    slicevnoT.append(slice(gtrajvals[it,2]-box2, gtrajvals[it,2]+    \
+                      box2+1))
+                    slicevnoT.append(slice(gtrajvals[it,1]-box2, gtrajvals[it,1]+    \
+                      box2+1))
+
+                    slice2D.append(slice(gtrajvals[it,2]-box2, gtrajvals[it,2]+box2 +\
+                      1))
+                    slice2D.append(slice(gtrajvals[it,1]-box2, gtrajvals[it,1]+box2 +\
+                      1))
+
+                    varvalst = varobj[tuple(slicev)]
+# box values
+
+                    varvals[it,:,:,:] = varvalst
+#                    print 'values at time t______'
+#                    print varvalst
+
+# box stats values
+                    statvarvals[it,:,0] = varvalst[box2,box2]
+                    statvarvals[it,:,1] = np.min(varvalst)
+                    statvarvals[it,:,2] = np.max(varvalst)
+                    statvarvals[it,:,3] = np.mean(varvalst)
+                    statvarvals[it,:,4] = np.mean(varvalst*varvalst)
+                    statvarvals[it,:,5] = np.sqrt(statvarvals[it,:,4] -              \
+                      statvarvals[it,:,3]*statvarvals[it,:,3])
+
+                    varvalst = lonv[tuple(slicevnoT)]
+                    lonvals[it,:,:] = varvalst
+                    varvalst = latv[tuple(slicevnoT)]           
+                    latvals[it,:,:] = varvalst
+
+# Circle values
+                cslicev.appen(slice(0:dimz))
+                cslicev.appen(slice(cyrangeslice[it]))
+                cslicev.appen(slice(cxrangeslice[it]))
+
+                cslicevnoT.appen(slice(0:dimz))
+                cslicevnoT.appen(slice(cyrangeslice[it]))
+                cslicevnoT.appen(slice(cxrangeslice[it]))
+
+                cslice2D.appen(slice(0:dimz))
+                cslice2D.appen(slice(cyrangeslice2D[it]))
+                cslice2D.appen(slice(cxrangeslice2D[it]))
+
+                if gtrajvals[it,2]-Nrad < 0 or gtrajvals[it,2]+Nrad + 1 > dimy + 1   \
+                  or gtrajvals[it,1]-Nrad < 0 or gtrajvals[it,1]+Nrad + 1 > dimx + 1:
+
+                    rvarvalst = np.ones((Nrad*2+1,Nrad*2+1),dtype=np.float)*fillValue
+                    rvarvalst[tuple(cslice2D)] = varobj[tuple(cslicev)]
+                    rvarvalst[tuple(cslice2D)] = np.where(circdist[tuple(cslice2D)] >\
+                      np.float(Nrad), fillValue, rvarvalst[tuple(cslice2D)])
+
+                    rvarvals[it,:,:,:] = rvarvalst
+
+# circle stats values
+                    maskedvals = ma.masked_values (rvarvalst, fillValue)
+                    rstatvarvals[it,:,0] = rvarvalst[Nrad,Nrad]
+                    rstatvarvals[it,:,1] = maskedvals.min()
+                    rstatvarvals[it,:,2] = maskedvals.max()
+                    rstatvarvals[it,:,3] = maskedvals.mean()
+                    maskedvals2 = maskedvals*maskedvals
+                    rstatvarvals[it,:,4] = maskedvals2.mean()
+                    rstatvarvals[it,:,5] = np.sqrt(rstatvarvals[it,:,4] -            \
+                      rstatvarvals[it,:,3]*rstatvarvals[it,:,3])
+
+                    rvarvalst[tuple(cslice2D)] = lonv[tuple(cslicevnoT)]
+                    rlonvals[it,:,:,:] = rvarvalst
+                    rvarvalst[tuple(cslice2D)] = latv[tuple(cslicevnoT)]           
+                    rlatvals[it,:,:,:] = rvarvalst
+
+                else:
+                    slicev.append(slice(gtrajvals[it,2]-Nrad,gtrajvals[it,2]+Nrad+1))
+                    slicev.append(slice(gtrajvals[it,1]-Nrad,gtrajvals[it,1]+Nrad+1))
+                    slicevnoT.append(slice(gtrajvals[it,2]-Nrad, gtrajvals[it,2]+    \
+                      Nrad+1))
+                    slicevnoT.append(slice(gtrajvals[it,1]-Nrad, gtrajvals[it,1]+    \
+                      Nrad+1))
+
+                    slice2D.append(slice(gtrajvals[it,2]-Nrad,gtrajvals[it,2]+Nrad+1))
+                    slice2D.append(slice(gtrajvals[it,1]-Nrad,gtrajvals[it,1]+Nrad+1))
+
+                    rvarvalst = varobj[tuple(cslicev)]
+                    cdist = circdist[tuple(cslicevnoT)]
+# circle values
+                    rvarvalst = np.where(cdist > np.float(Nrad),fillValue,rvarvalst)
+
+                    rvarvals[it,:,:,:] = rvarvalst
+
+# circle stats values
+                    maskedvals = ma.masked_values (rvarvalst, fillValue)
+                    rstatvarvals[it,:,0] = rvarvalst[Nrad,Nrad]
+                    rstatvarvals[it,:,1] = maskedvals.min()
+                    rstatvarvals[it,:,2] = maskedvals.max()
+                    rstatvarvals[it,:,3] = maskedvals.mean()
+                    maskedvals2 = maskedvals*maskedvals
+                    rstatvarvals[it,:,4] = maskedvals2.mean()
+                    rstatvarvals[it,:,5] = np.sqrt(rstatvarvals[it,4] -              \
+                      rstatvarvals[it,:,3]*rstatvarvals[it,3])
+
+                    rvarvalst = lonv[tuple(cslicevnoT)]
+                    rlonvals[it,:,:,:] = rvarvalst
+                    rvarvalst = latv[tuple(cslicevnoT)]           
+                    rlatvals[it,:,:,:] = rvarvalst
+
+#            print 'statistics:',rstatvarvals[it,:] 
+
 # variable box values
-    newvar = objofile.createVariable(varn + 'box', 'f4', ('time', 'y', 'x'),         \
-      fill_value=fillValue)
-    if searchInlist(varobj.ncattrs(),'standard_name'):
-        vsname = varobj.getncattr('standard_name')
-    else:
-        vsname = variables_values(varn)[1]
-    if searchInlist(varobj.ncattrs(),'long_name'):
-        vlname = varobj.getncattr('long_name')
-    else:
-        vlname = variables_values(varn)[4].replace('|',' ')
-    if searchInlist(varobj.ncattrs(),'units'):
-        vunits = varobj.getncattr('units')
-    else:
-        vunits = variables_values(varn)[5].replace('|',' ')
-
-    newattr = basicvardef(newvar, vsname, vlname, vunits)
-    newattr = newvar.setncattr('projection','lon lat')
-    newvar[:] = varvals
+            newvar = objofile.createVariable(varn + 'box', 'f4', ('time','z','y','x'),\
+              fill_value=fillValue)
+            if searchInlist(varobj.ncattrs(),'standard_name'):
+                vsname = varobj.getncattr('standard_name')
+            else:
+                vsname = variables_values(varn)[1]
+            if searchInlist(varobj.ncattrs(),'long_name'):
+                vlname = varobj.getncattr('long_name')
+            else:
+                vlname = variables_values(varn)[4].replace('|',' ')
+            if searchInlist(varobj.ncattrs(),'units'):
+                vunits = varobj.getncattr('units')
+            else:
+                vunits = variables_values(varn)[5].replace('|',' ')
+
+            newattr = basicvardef(newvar, vsname, vlname, vunits)
+            newattr = newvar.setncattr('projection','lon lat')
+            newvar[:] = varvals
 
 # center of the trajectory
-    newvar = objofile.createVariable('trj_' + varn, 'f', ('time'),                  \
-      fill_value=fillValue)
-    newattr = basicvardef(newvar, 'trj_' + vsname, 'value along the trajectory of '+\
-      varn, vunits)
-    newvar[:] = statvarvals[:,0]
+            newvar = objofile.createVariable('trj_' + varn, 'f', ('time','z'),       \
+              fill_value=fillValue)
+            newattr = basicvardef(newvar, 'trj_' + vsname, 'value along the ' +      \
+              'trajectory of '+ varn, vunits)
+            newvar[:] = statvarvals[:,:,0]
 
 # variable box statistics
-    ist = 0
-
-    statnames = ['minbox', 'maxbox', 'meanbox', 'mean2box', 'stdevbox']
-    vstlname = ['minimum value within', 'maximum value within',                      \
-     'mean value within', 'squared mean value within',                               \
-     'standard deviation value within']
-
-    for statn in statnames:
-        newvar = objofile.createVariable(statn + '_' + varn, 'f', ('time'),          \
-          fill_value=fillValue)
-        newattr = basicvardef(newvar, statn + '_' + vsname, vstlname[ist] +          \
-          ' the box ('+str(boxs)+'x'+str(boxs)+') of ' + varn, vunits)
-        newvar[:] = statvarvals[:,ist+1]
-#        newattr = set_attributek(newvar,'_FillValue',fillValue,'npfloat')
-        ist = ist + 1
+            ist = 0
+            for statn in statnames:
+                newvar = objofile.createVariable(statn + '_' + varn,'f',('time','z'),\
+                  fill_value=fillValue)
+                newattr = basicvardef(newvar, statn + '_' + vsname, vstlname[ist] +  \
+                  ' the box ('+str(boxs)+'x'+str(boxs)+') of ' + varn, vunits)
+                newvar[:] = statvarvals[:,:,ist+1]
+#                newattr = set_attributek(newvar,'_FillValue',fillValue,'npfloat')
+                ist = ist + 1
         
-
 # variable circle values
-    newvar = objofile.createVariable(varn + 'circle', 'f4', ('time', 'yr', 'xr'),      \
-      fill_value=fillValue)
-    if searchInlist(varobj.ncattrs(),'standard_name'):
-        vsname = varobj.getncattr('standard_name')
-    else:
-        vsname = variables_values(varn)[1]
-    if searchInlist(varobj.ncattrs(),'long_name'):
-        vlname = varobj.getncattr('long_name')
-    else:
-        vlname = variables_values(varn)[4].replace('|',' ')
-    if searchInlist(varobj.ncattrs(),'units'):
-        vunits = varobj.getncattr('units')
-    else:
-        vunits = variables_values(varn)[5].replace('|',' ')
-
-    newattr = basicvardef(newvar, vsname, vlname, vunits)
-    newattr = newvar.setncattr('projection','lon lat')
-    newvar[:] = rvarvals
+            newvar = objofile.createVariable(varn + 'circle','f4',('time','z','yr',  \
+              'xr'), fill_value=fillValue)
+            if searchInlist(varobj.ncattrs(),'standard_name'):
+                vsname = varobj.getncattr('standard_name')
+            else:
+                vsname = variables_values(varn)[1]
+            if searchInlist(varobj.ncattrs(),'long_name'):
+                vlname = varobj.getncattr('long_name')
+            else:
+                vlname = variables_values(varn)[4].replace('|',' ')
+            if searchInlist(varobj.ncattrs(),'units'):
+                vunits = varobj.getncattr('units')
+            else:
+                vunits = variables_values(varn)[5].replace('|',' ')
+
+            newattr = basicvardef(newvar, vsname, vlname, vunits)
+            newattr = newvar.setncattr('projection','lon lat')
+            newvar[:] = rvarvals
 
 # variable circle statistics
-    ist = 0
-    statnames = ['mincircle', 'maxcircle', 'meancircle', 'mean2circle', 'stdevcircle']
-
-    for statn in statnames:
-        newvar = objofile.createVariable(statn + '_' + varn, 'f', ('time'),          \
-          fill_value=fillValue)
-        newattr = basicvardef(newvar, statn + '_' + vsname, vstlname[ist] +          \
-          ' the circle of radius ('+ str(circler)+') of ' + varn, vunits)
-        newvar[:] = rstatvarvals[:,ist+1]
-#        newattr = set_attributek(newvar,'_FillValue',fillValue,'npfloat')
-        ist = ist + 1
+            ist = 0
+            for statn in cstatnames:
+                newvar = objofile.createVariable(statn + '_' + varn,'f',('time','z'),\
+                  fill_value=fillValue)
+                newattr = basicvardef(newvar, statn + '_' + vsname, vstlname[ist] +  \
+                  ' the circle of radius ('+ str(circler)+') of ' + varn, vunits)
+                newvar[:] = rstatvarvals[:,:,ist+1]
+#                newattr = set_attributek(newvar,'_FillValue',fillValue,'npfloat')
+                ist = ist + 1
+
+        elif Nvardims == 3:
+# Too optimistic, but at this stage...
+            dimt = varobj.shape[0]
+ 
+            varvals = np.ones(tuple([Ttraj,boxs,boxs]), dtype=np.float)
+            lonvals = np.ones(tuple([Ttraj,boxs,boxs]), dtype=np.float)
+            latvals = np.ones(tuple([Ttraj,boxs,boxs]), dtype=np.float)
+            rvarvals = np.ones(tuple([Ttraj,Nrad*2+1,Nrad*2+1]), dtype=np.float)
+            rlonvals = np.ones(tuple([Ttraj,Nrad*2+1,Nrad*2+1]), dtype=np.float)
+            rlatvals = np.ones(tuple([Ttraj,Nrad*2+1,Nrad*2+1]), dtype=np.float)
+
+            statvarvals = np.ones(tuple([Ttraj,6]), dtype=np.float)
+            rstatvarvals = np.ones(tuple([Ttraj,6]), dtype=np.float)
+
+            for it in range(dimt):
+# box values
+                slicev.appen(slice(yrangeslice[it]))
+                slicev.appen(slice(xrangeslice[it]))
+
+                slicevnoT.appen(slice(yrangeslice[it]))
+                slicevnoT.appen(slice(xrangeslice[it]))
+
+                slice2D.appen(slice(yrangeslice2D[it]))
+                slice2D.appen(slice(xrangeslice2D[it]))
+
+                if gtrajvals[it,2]-box2 < 0 or gtrajvals[it,2]+box2 + 1 > dimy + 1   \
+                  or gtrajvals[it,1]-box2 < 0 or gtrajvals[it,1]+box2 + 1 > dimx + 1:
+
+                    rvarvalst = np.ones((Nrad*2+1, Nrad*2+1),dtype=np.float)*fillValue
+
+                    varvalst[tuple(slice2D)] = varobj[tuple(slicev)]
+                    varvals[it,:,:] = varvalst
+    
+# box stats values
+                    maskedvals = ma.masked_values (varvalst, fillValue)
+                    statvarvals[it,0] = varvalst[box2,box2]
+                    statvarvals[it,1] = maskedvals.min()
+                    statvarvals[it,2] = maskedvals.max()
+                    statvarvals[it,3] = maskedvals.mean()
+                    maskedvals2 = maskedvals*maskedvals
+                    statvarvals[it,4] = maskedvals2.mean()
+                    statvarvals[it,5] = np.sqrt(statvarvals[it,4] -                   \
+                      statvarvals[it,3]*statvarvals[it,3])
+
+                    varvalst[tuple(slice2D)] = lonv[tuple(slicevnoT)]
+                    lonvals[it,:,:] = varvalst
+                    varvalst[tuple(slice2D)] = latv[tuple(slicevnoT)]           
+                    latvals[it,:,:] = varvalst
+
+                else:
+                    slicev.append(slice(gtrajvals[it,2]-box2, gtrajvals[it,2]+box2+1))
+                    slicev.append(slice(gtrajvals[it,1]-box2, gtrajvals[it,1]+box2+1))
+                    slicevnoT.append(slice(gtrajvals[it,2]-box2, gtrajvals[it,2]+    \
+                      box2+1))
+                    slicevnoT.append(slice(gtrajvals[it,1]-box2, gtrajvals[it,1]+    \
+                      box2+1))
+
+                    slice2D.append(slice(gtrajvals[it,2]-box2, gtrajvals[it,2]+box2 +\
+                      1))
+                    slice2D.append(slice(gtrajvals[it,1]-box2, gtrajvals[it,1]+box2 +\
+                      1))
+
+                    varvalst = varobj[tuple(slicev)]
+# box values
+
+                    varvals[it,:,:] = varvalst
+#                    print 'values at time t______'
+#                    print varvalst
+
+# box stats values
+                    statvarvals[it,0] = varvalst[box2,box2]
+                    statvarvals[it,1] = np.min(varvalst)
+                    statvarvals[it,2] = np.max(varvalst)
+                    statvarvals[it,3] = np.mean(varvalst)
+                    statvarvals[it,4] = np.mean(varvalst*varvalst)
+                    statvarvals[it,5] = np.sqrt(statvarvals[it,4] -                  \
+                      statvarvals[it,3]*statvarvals[it,3])
+
+                    varvalst = lonv[tuple(slicevnoT)]
+                    lonvals[it,:,:] = varvalst
+                    varvalst = latv[tuple(slicevnoT)]           
+                    latvals[it,:,:] = varvalst
+
+# Circle values
+                cslicev.appen(slice(cyrangeslice[it]))
+                cslicev.appen(slice(cxrangeslice[it]))
+
+                cslicevnoT.appen(slice(cyrangeslice[it]))
+                cslicevnoT.appen(slice(cxrangeslice[it]))
+
+                cslice2D.appen(slice(cyrangeslice2D[it]))
+                cslice2D.appen(slice(cxrangeslice2D[it]))
+                if gtrajvals[it,2]-Nrad < 0 or gtrajvals[it,2]+Nrad + 1 > dimy + 1   \
+                  or gtrajvals[it,1]-Nrad < 0 or gtrajvals[it,1]+Nrad + 1 > dimx + 1:
+
+                    rvarvalst = np.ones((Nrad*2+1,Nrad*2+1),dtype=np.float)*fillValue
+                    rvarvalst[tuple(cslice2D)] = varobj[tuple(cslicev)]
+                    rvarvalst[tuple(cslice2D)] = np.where(circdist[tuple(cslice2D)] >\
+                      np.float(Nrad), fillValue, rvarvalst[tuple(cslice2D)])
+
+                    rvarvals[it,:,:] = rvarvalst
+
+# circle stats values
+                    maskedvals = ma.masked_values (rvarvalst, fillValue)
+                    rstatvarvals[it,0] = rvarvalst[Nrad,Nrad]
+                    rstatvarvals[it,1] = maskedvals.min()
+                    rstatvarvals[it,2] = maskedvals.max()
+                    rstatvarvals[it,3] = maskedvals.mean()
+                    maskedvals2 = maskedvals*maskedvals
+                    rstatvarvals[it,4] = maskedvals2.mean()
+                    rstatvarvals[it,5] = np.sqrt(rstatvarvals[it,4] -                \
+                      rstatvarvals[it,3]*rstatvarvals[it,3])
+
+                    rvarvalst[tuple(cslice2D)] = lonv[tuple(cslicevnoT)]
+                    rlonvals[it,:,:] = rvarvalst
+                    rvarvalst[tuple(cslice2D)] = latv[tuple(cslicevnoT)]           
+                    rlatvals[it,:,:] = rvarvalst
+
+                else:
+                    slicev.append(slice(gtrajvals[it,2]-Nrad,gtrajvals[it,2]+Nrad+1))
+                    slicev.append(slice(gtrajvals[it,1]-Nrad,gtrajvals[it,1]+Nrad+1))
+                    slicevnoT.append(slice(gtrajvals[it,2]-Nrad, gtrajvals[it,2]+    \
+                      Nrad+1))
+                    slicevnoT.append(slice(gtrajvals[it,1]-Nrad, gtrajvals[it,1]+    \
+                      Nrad+1))
+
+                    slice2D.append(slice(gtrajvals[it,2]-Nrad,gtrajvals[it,2]+Nrad+1))
+                    slice2D.append(slice(gtrajvals[it,1]-Nrad,gtrajvals[it,1]+Nrad+1))
+
+                    rvarvalst = varobj[tuple(cslicev)]
+                    cdist = circdist[tuple(cslicevnoT)]
+# circle values
+                    rvarvalst = np.where(cdist > np.float(Nrad),fillValue,rvarvalst)
+
+                    rvarvals[it,:,:] = rvarvalst
+
+# circle stats values
+                    maskedvals = ma.masked_values (rvarvalst, fillValue)
+                    rstatvarvals[it,0] = rvarvalst[Nrad,Nrad]
+                    rstatvarvals[it,1] = maskedvals.min()
+                    rstatvarvals[it,2] = maskedvals.max()
+                    rstatvarvals[it,3] = maskedvals.mean()
+                    maskedvals2 = maskedvals*maskedvals
+                    rstatvarvals[it,4] = maskedvals2.mean()
+                    rstatvarvals[it,5] = np.sqrt(rstatvarvals[it,4] -                \
+                      rstatvarvals[it,3]*rstatvarvals[it,3])
+
+                    rvarvalst = lonv[tuple(cslicevnoT)]
+                    rlonvals[it,:,:] = rvarvalst
+                    rvarvalst = latv[tuple(cslicevnoT)]           
+                    rlatvals[it,:,:] = rvarvalst
+
+#            print 'statistics:',rstatvarvals[it,:] 
+
+# variable box values
+            newvar = objofile.createVariable(varn + 'box', 'f4', ('time', 'y', 'x'), \
+              fill_value=fillValue)
+            if searchInlist(varobj.ncattrs(),'standard_name'):
+                vsname = varobj.getncattr('standard_name')
+            else:
+                vsname = variables_values(varn)[1]
+            if searchInlist(varobj.ncattrs(),'long_name'):
+                vlname = varobj.getncattr('long_name')
+            else:
+                vlname = variables_values(varn)[4].replace('|',' ')
+            if searchInlist(varobj.ncattrs(),'units'):
+                vunits = varobj.getncattr('units')
+            else:
+                vunits = variables_values(varn)[5].replace('|',' ')
+
+            newattr = basicvardef(newvar, vsname, vlname, vunits)
+            newattr = newvar.setncattr('projection','lon lat')
+            newvar[:] = varvals
+
+# center of the trajectory
+            newvar = objofile.createVariable('trj_' + varn, 'f', ('time'),           \
+              fill_value=fillValue)
+            newattr = basicvardef(newvar, 'trj_' + vsname, 'value along the ' +      \
+              'trajectory of '+ varn, vunits)
+            newvar[:] = statvarvals[:,0]
+
+# variable box statistics
+            ist = 0
+            for statn in statnames:
+                newvar = objofile.createVariable(statn + '_' + varn, 'f', ('time'),  \
+                  fill_value=fillValue)
+                newattr = basicvardef(newvar, statn + '_' + vsname, vstlname[ist] +  \
+                  ' the box ('+str(boxs)+'x'+str(boxs)+') of ' + varn, vunits)
+                newvar[:] = statvarvals[:,ist+1]
+#                newattr = set_attributek(newvar,'_FillValue',fillValue,'npfloat')
+                ist = ist + 1
+        
+# variable circle values
+            newvar = objofile.createVariable(varn + 'circle','f4',('time','yr','xr'),\
+              fill_value=fillValue)
+            if searchInlist(varobj.ncattrs(),'standard_name'):
+                vsname = varobj.getncattr('standard_name')
+            else:
+                vsname = variables_values(varn)[1]
+            if searchInlist(varobj.ncattrs(),'long_name'):
+                vlname = varobj.getncattr('long_name')
+            else:
+                vlname = variables_values(varn)[4].replace('|',' ')
+            if searchInlist(varobj.ncattrs(),'units'):
+                vunits = varobj.getncattr('units')
+            else:
+                vunits = variables_values(varn)[5].replace('|',' ')
+
+            newattr = basicvardef(newvar, vsname, vlname, vunits)
+            newattr = newvar.setncattr('projection','lon lat')
+            newvar[:] = rvarvals
+
+# variable circle statistics
+            ist = 0
+            for statn in cstatnames:
+                newvar = objofile.createVariable(statn + '_' + varn, 'f', ('time'),  \
+                  fill_value=fillValue)
+                newattr = basicvardef(newvar, statn + '_' + vsname, vstlname[ist] +  \
+                  ' the circle of radius ('+ str(circler)+') of ' + varn, vunits)
+                newvar[:] = rstatvarvals[:,ist+1]
+#                newattr = set_attributek(newvar,'_FillValue',fillValue,'npfloat')
+                ist = ist + 1
 
 # global attributes
@@ -10253,8 +10568,11 @@
 
 #compute_tevolboxtraj('h', 'wrfout*', 'PSFC')
-#compute_tevolboxtraj('control/trajectory.dat@0,XLONG,XLAT,Times,wrf,3,3',            \
+#compute_tevolboxtraj('control/trajectory.dat@0,XLONG,XLAT,ZNU,Times,wrf,3,3',       \
 #  '../../superstorm/control/wrfout/wrfout_d01_2001-11-09_00:00:00', 'PSFC')
 #compute_tevolboxtraj('control/trajectory.dat@0,lon,lat,time,cf,5,5',                 \
 #  'control/wss.nc', 'wss')
+# in /home/lluis/etudes/WL_HyMeX/superstorm
+compute_tevolboxtraj('control/trajectory.dat@0,XLONG,XLAT,ZNU,Times,wrf,3,3',       \
+  'control/wrfout/wrfout_d01_2001-11-09_00:00:00', 'PSFC')
 
 def numVector_String(vec,char):
