Index: /trunk/tools/documentation/ncmanage/compute_opersvarsfiles.html
===================================================================
--- /trunk/tools/documentation/ncmanage/compute_opersvarsfiles.html	(revision 1560)
+++ /trunk/tools/documentation/ncmanage/compute_opersvarsfiles.html	(revision 1561)
@@ -38,4 +38,23 @@
       &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;'divc',[modval1]: [prevalues] divide by [modval1]<BR>
       &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;'forwrdderiv',[N],[ord],[dim]: un-scaled forward [N]-derivative of order [ord] along dimension [dim] of [var]<BR>
+      &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;'ifreq_anom',[stepdimn],[stepvardimn],[istep]: computing anomalies by substracting sub-means at each<BR>
+      &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;step by averaging from there all values taken every [istep] along dimension [stepdim]. <BR>
+      &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;mean(j) = sum(matA[j+k*istep]_k=0,Nstep)/Nstep; Nstep = len(stepdimn)/istep; j=[0,istep]<BR>
+      &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;[stepdimn]= name of the dimension along which to sample<BR>
+      &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;[stepvardimn]= name of the variable-dimension with the values for [stepdimn]<BR>
+      &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;[istep]= frequency to sample<BR>
+      &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;'ifreq_mean',[stepdimn],[stepvardimn],[istep]: computing sub-means at each step by averaging from <BR>
+      &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;there all values taken every [istep] along dimension [stepdim]. <BR>
+      &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;mean(j) = sum(matA[j+k*istep]_k=0,Nstep)/Nstep; Nstep = len(stepdimn)/istep; j=[0,istep]<BR>
+      &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;[stepdimn]= name of the dimension along which to sample<BR>
+      &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;[stepvardimn]= name of the variable-dimension with the values for [stepdimn]<BR>
+      &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;[istep]= frequency to sample<BR>
+      &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;'ifreq_normmeanstd',[stepdimn],[stepvardimn],[istep]: nbormalizing anomalies by substracting mean(j)/dtsv(j) <BR>
+      &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;sub-stats at each step by averaging from there all values taken every [istep] along dimension[stepdim]. <BR>
+      &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;mean(j) = sum(matA[j+k*istep]_k=0,Nstep)/Nstep; std(j) = std(matA[j+k*istep]_k=0,Nstep); <BR>
+      &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Nstep = len(stepdimn)/istep; j=[0,istep]<BR>
+      &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;[stepdimn]= name of the dimension along which to sample<BR>
+      &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;[stepvardimn]= name of the variable-dimension with the values for [stepdimn]<BR>
+      &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;[istep]= frequency to sample<BR>
       &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;'inv': inverting [prevalues] (1/[prevalues])<BR>
       &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;'lowthres',[modval1],[modval2]: if [prevalues] < [modval1]; prevalues = [modval2]<BR>
@@ -49,5 +68,5 @@
       &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;'norm_meanstd',[NOTnormdims]: normalization of data as: (val-<val>)/stdev(val) except along <BR>
       &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;dimensions [NOTnormdims] (':' list of dimension names, or 'any' for using all dimensions)<BR>
-       &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;dimensions [NOTnormdims] (':' list of dimension names)<BR>
+       &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;dimensions [NOTnormdims] (':' list of dimension names)<BR>
       &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;'pot': powering with [var] ([prevalues] ** [var])<BR>
       &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;'potc',[modval1]: [prevalues] ** [modval1]<BR>
@@ -72,13 +91,16 @@
     <DIV CLASS="valins">
     * Transforming temperature from Kelvin to &deg;C<BR>
-    &nbsp;&nbsp;  $ python ${pyHOME}/nc_var.py -o compute_opersvarsfiles -S 'west_east|XLONG|-1;south_north|XLAT|-1;Time|Times|3@subc,273.15|wrfout_d01_2001-11-11_00:00:00|T2' -v 'tempC,air!temperature,C'<BR>
+    &nbsp;&nbsp;$ python ${pyHOME}/nc_var.py -o compute_opersvarsfiles -S 'west_east|XLONG|-1;south_north|XLAT|-1;Time|Times|3@subc,273.15|wrfout_d01_2001-11-11_00:00:00|T2' -v 'tempC,air!temperature,C'<BR>
     * Computing the x-derivative of first order<BR>
-    &nbsp;&nbsp;  $ python ${pyHOME}/nc_var.py -o compute_opersvarsfiles -S 'lon|lon|-1;lat|lat|-1;time_counter|time_counter|-1@forwrdderiv,1,1,2|histday.nc|t2m' -v 'tasderiv,x-derivative|of|air|temperature,K' <BR>
+    &nbsp;&nbsp;$ python ${pyHOME}/nc_var.py -o compute_opersvarsfiles -S 'lon|lon|-1;lat|lat|-1;time_counter|time_counter|-1@forwrdderiv,1,1,2|histday.nc|t2m' -v 'tasderiv,x-derivative|of|air|temperature,K' <BR>
     * Normalizing a variable by substracting its mean and weighting by its standard-deviation<BR>
-    &nbsp;&nbsp; $ python ${pyHOME}/nc_var.py -o compute_opersvarsfiles -S 'west_east|XLONG|-1;south_north|XLAT|-1;Time|WRFtime|-1@norm_meanstd,Time|wrfout_d01_1995-01-01_00:00:00|T2' -v 'tasnorm,normalized!2m!temperature!substracting!mean!and!weighting!by!standard!deviation,K' <BR>
+    &nbsp;&nbsp;$ python ${pyHOME}/nc_var.py -o compute_opersvarsfiles -S 'west_east|XLONG|-1;south_north|XLAT|-1;Time|WRFtime|-1@norm_meanstd,Time|wrfout_d01_1995-01-01_00:00:00|T2' -v 'tasnorm,normalized!2m!temperature!substracting!mean!and!weighting!by!standard!deviation,K' <BR>
     * Getting height of the frist level from surface from a WRF file<BR>
-    &nbsp;&nbsp; $ python $pyHOME/nc_var.py -o compute_opersvarsfiles -S 'west_east|XLONG|-1;south_north|XLAT|-1;bottom_top_stag|ZNW|1@addc,0|wrfout_d01_1995-01-01_00:00:00|PH%west_east|XLONG|-1;south_north|XLAT|-1;bottom_top_stag|ZNW|1@add|wrfout_d01_1995-01-01_00:00:00|PHB%contoperation@divc,9.81%west_east|XLONG|-1;south_north|XLAT|-1;Time|WRFtime|-1@sub|wrfout_d01_1995-01-01_00:00:00|HGT' -v 'height1lev,height!above!surface!of!first!level,m'<BR>
+    &nbsp;&nbsp;$ python $pyHOME/nc_var.py -o compute_opersvarsfiles -S 'west_east|XLONG|-1;south_north|XLAT|-1;bottom_top_stag|ZNW|1@addc,0|wrfout_d01_1995-01-01_00:00:00|PH%west_east|XLONG|-1;south_north|XLAT|-1;bottom_top_stag|ZNW|1@add|wrfout_d01_1995-01-01_00:00:00|PHB%contoperation@divc,9.81%west_east|XLONG|-1;south_north|XLAT|-1;Time|WRFtime|-1@sub|wrfout_d01_1995-01-01_00:00:00|HGT' -v 'height1lev,height!above!surface!of!first!level,m'<BR>
     * Computing wind-direction from a WRF file<BR>
-    &nbsp;&nbsp; $ python $pyHOME/nc_var.py -o compute_opersvarsfiles -S 'west_east|XLONG|-1;south_north|XLAT|-1;Time|WRFtime|-1@addc,0|wrfout_d01_1995-01-01_00:00:00|V10%west_east|XLONG|-1;south_north|XLAT|-1;Time|WRFtime|-1@arctan|wrfout_d01_1995-01-01_00:00:00|U10%contoperation@mulc,57.2957795131' -v 'wsdir,2m!wind!direction,Degrees'
+    &nbsp;&nbsp;$ python $pyHOME/nc_var.py -o compute_opersvarsfiles -S 'west_east|XLONG|-1;south_north|XLAT|-1;Time|WRFtime|-1@addc,0|wrfout_d01_1995-01-01_00:00:00|V10%west_east|XLONG|-1;south_north|XLAT|-1;Time|WRFtime|-1@arctan|wrfout_d01_1995-01-01_00:00:00|U10%contoperation@mulc,57.2957795131' -v 'wsdir,2m!wind!direction,Degrees'
+    * Computing normalized anomalyes of WRF outputs by the mean and standard deviations of the outputs at the same hour<BR>
+    &nbsp;&nbsp;$ python nc_var.py -o compute_opersvarsfiles -S 'Time|WRFtime|-1;south_north|XLAT|-1;west_east|XLONG|-1@ifreq_normmeanstd,Time,WRFtime,8|/home/lluis/PY/wrfout_d01_1995-01-01_00:00:00|T2' -v 'tas_ifreqnormmeanstd,tas!anomaly!by!substracting!frequency!mean!at!every!8!time-steps,K'
+
     </DIV>
   </BODY>
Index: /trunk/tools/generic_tools.py
===================================================================
--- /trunk/tools/generic_tools.py	(revision 1560)
+++ /trunk/tools/generic_tools.py	(revision 1561)
@@ -10206,4 +10206,16 @@
         'divc',[modval1]: divide by [modval1]
         'forwrdderiv',[N],[ord],[dim]: un-scaled forward [N]-derivative of order [ord] along dimension [dim] of matA
+        'ifreq_anom',[dimns],[istep],[stepdimn]: computing anomalues by substracting sub-means at each step by
+             averaging from there all values taken every [istep] along dimension [stepdim]. 
+             mean(j) = sum(matA[j+k*istep]_k=0,Nstep)/Nstep; Nstep = len(stepdimn)/istep; j=[0,istep]
+          [dimns]= ':' list of dimensions of the matrice
+          [istep]= frequency to sample
+          [stepdimn]= name of the dimension along which to sample
+        'ifreq_mean',[dimns],[istep],[stepdimn]: computing sub-means at each step by averaging from there 
+           all values taken every [istep] along dimension [stepdim]. mean(j) = sum(matA[j+k*istep]_k=0,Nstep)/Nstep;
+           Nstep = len(stepdimn)/istep; j=[0,istep]
+          [dimns]= ':' list of dimensions of the matrice
+          [istep]= frequency to sample
+          [stepdimn]= name of the dimension along which to sample
         'inv': inverting matA (1/matA)
         'lowthres',[modval1],[modval2]: if [val] < [modval1]; val = [modval2]
@@ -10350,6 +10362,7 @@
     # Dictionary with the description of the operations
     Lopn = {'add': '+', 'addc': '+', 'arctan': 'arctangent', 'div': '/', 'divc': '/',\
-      'inv': '^(-1)', 'lowthres': '<', 'mul': '*', 'mulc': '*',                      \
-      'norm_meanminmax': 'Norm_mean_minmax',                                         \
+      'ifreq_mean': 'ifreq_mean', 'ifreq_mean': 'ifreq_mean', 'inv': '^(-1)',        \
+      'lowthres': '<', 'mul': '*',                                                   \
+      'mulc': '*', 'norm_meanminmax': 'Norm_mean_minmax',                            \
       'norm_minmax': 'Norm_minmax', 'norm_meanstd': 'Norm_mean_devstd', 'pot': '^',  \
       'potc': '^', 'repl': 'repl_excpt', 'replbig': 'repl_excpt_bigger',             \
@@ -10362,5 +10375,5 @@
     # Available operations
     opavail = ['add', 'addc', 'arctan', 'centerderiv', 'div', 'divc', 'forwrdderiv', \
-      'inv',    \
+      'ifreq_anom', 'ifreq_mean', 'inv',                                             \
       'lowthres', 'lowthres@oper', 'mul', 'mulc', 'norm_meanminmax', 'norm_minmax',  \
       'norm_meanstd', 'pot', 'potc', 'repl', 'replbig', 'repless', 'sub', 'subc',    \
@@ -10431,4 +10444,31 @@
             quit(-1)
         derivcoeff = forwrdderiv[Nderiv-1,order-1,:]
+    elif opN[0:10] == 'ifreq_anom':
+        opn = 'ifreq_anom'
+        dimns = opN.split(',')[1].split(':')
+        istep = int(opN.split(',')[2])
+        stepdimn = opN.split(',')[3]
+        Lopn['ifreq_anom'] = 'anomaly by substracting ifreq_mean = sum(mat[j+k*' +   \
+          str(istep)+']_k=0,Nstep)/Nstep; ' + 'Nstep=len(' + stepdimn + ')/' +       \
+          str(istep) + '; j=[0,' + str(istep) + ']'
+    elif opN[0:10] == 'ifreq_mean':
+        opn = 'ifreq_mean'
+        dimns = opN.split(',')[1].split(':')
+        istep = int(opN.split(',')[2])
+        stepdimn = opN.split(',')[3]
+        Lopn['ifreq_mean'] = 'sum(mat[j+k*' +str(istep)+']_k=0,Nstep)/Nstep; ' +     \
+          'Nstep=len(' + stepdimn + ')/' + str(istep) + '; j=[0,' + str(istep) + ']'
+    elif opN[0:17] == 'ifreq_normmeanstd':
+        opn = 'ifreq_normmeanstd'
+        dimns = opN.split(',')[1].split(':')
+        istep = int(opN.split(',')[2])
+        stepdimn = opN.split(',')[3]
+        Lopn['ifreq_normmeanstd'] = 'normalization by substracting ' +               \
+          'ifreq_mean/ifreq_stdv ifreq_mean= sum(mat[j+k*' + str(istep) +            \
+          ']_k=0,Nstep)/Nstep; ifreq_stdv=stdv(mat[j+k*' + str(istep)+']_k=0,Nstep)'+\
+          'Nstep=len(' + stepdimn + ')/' + str(istep) + '; j=[0,' + str(istep) + ']'
+    elif opN[0:4] == 'inv':
+        opn = 'inv'
+        Lopn['inv'] = 'inverse'
     elif opN[0:8] == 'lowthres' and opN.find('@') == -1:
         opn = 'lowthres'
@@ -10662,4 +10702,109 @@
               matA[tuple(mat5)]*derivcoeff[5] + matA[tuple(mat6)]*derivcoeff[6] +    \
               matA[tuple(mat7)]*derivcoeff[7] + matA[tuple(mat8)]*derivcoeff[8]
+    elif opn == 'ifreq_anom':
+        istepdim = dimns.index(stepdimn)
+        lenstepdim = matA.shape[istepdim]
+        Nsteps = lenstepdim/istep
+
+        meannewmatsize = list(matA.shape)
+        meannewmatsize[istepdim] = istep
+        meannewmat = np.zeros(tuple(meannewmatsize), dtype=valtype)
+        matAshape = matA.shape
+        NmatAdims = len(matAshape)
+
+        newmat = np.zeros(tuple(matAshape), dtype=valtype)
+        for j in range(istep):
+            slicemat = []
+            slicenewmat = []
+            meandims = []
+            for idim in range(NmatAdims):
+                if idim == istepdim:
+                    slicemat.append(slice(j,j+Nsteps*istep,istep))
+                    slicenewmat.append(j)
+                    meandims.append(idim)
+                else:
+                    slicemat.append(slice(0,matAshape[idim]))
+                    slicenewmat.append(slice(0,matAshape[idim]))
+            mdims=tuple(meandims)
+            meannewmat[tuple(slicenewmat)]= np.mean(matA[tuple(slicemat)], axis=mdims)
+
+            # Computing anomaly
+            for it in range(Nsteps):
+                slicemat = []
+                for idim in range(NmatAdims):
+                    if idim == istepdim:
+                        slicemat.append(j+it*istep)
+                    else:
+                        slicemat.append(slice(0,matAshape[idim]))
+                newmat[tuple(slicemat)] = matA[tuple(slicemat)] -                    \
+                  meannewmat[tuple(slicenewmat)]
+        additional = [meannewmat]
+
+    elif opn == 'ifreq_mean':
+        istepdim = dimns.index(stepdimn)
+        lenstepdim = matA.shape[istepdim]
+        Nsteps = lenstepdim/istep
+
+        newmatsize = list(matA.shape)
+        newmatsize[istepdim] = istep
+        newmat = np.zeros(tuple(newmatsize), dtype=valtype)
+        matAshape = matA.shape
+        NmatAdims = len(matAshape)
+
+        for j in range(istep):
+            slicemat = []
+            slicenewmat = []
+            meandims = []
+            for idim in range(NmatAdims):
+                if idim == istepdim:
+                    slicemat.append(slice(j,j+Nsteps*istep,istep))
+                    slicenewmat.append(j)
+                    meandims.append(idim)
+                else:
+                    slicemat.append(slice(0,matAshape[idim]))
+                    slicenewmat.append(slice(0,matAshape[idim]))
+            mdims=tuple(meandims)
+            newmat[tuple(slicenewmat)] = np.mean(matA[tuple(slicemat)], axis=mdims)
+    elif opn == 'ifreq_normmeanstd':
+        istepdim = dimns.index(stepdimn)
+        lenstepdim = matA.shape[istepdim]
+        Nsteps = lenstepdim/istep
+
+        meannewmatsize = list(matA.shape)
+        meannewmatsize[istepdim] = istep
+        meannewmat = np.zeros(tuple(meannewmatsize), dtype=valtype)
+        stdnewmat = np.zeros(tuple(meannewmatsize), dtype=valtype)
+        matAshape = matA.shape
+        NmatAdims = len(matAshape)
+
+        newmat = np.zeros(tuple(matAshape), dtype=valtype)
+        for j in range(istep):
+            slicemat = []
+            slicenewmat = []
+            meandims = []
+            for idim in range(NmatAdims):
+                if idim == istepdim:
+                    slicemat.append(slice(j,j+Nsteps*istep,istep))
+                    slicenewmat.append(j)
+                    meandims.append(idim)
+                else:
+                    slicemat.append(slice(0,matAshape[idim]))
+                    slicenewmat.append(slice(0,matAshape[idim]))
+            mdims=tuple(meandims)
+            meannewmat[tuple(slicenewmat)]= np.mean(matA[tuple(slicemat)], axis=mdims)
+            stdnewmat[tuple(slicenewmat)]= np.std(matA[tuple(slicemat)], axis=mdims)
+
+            # Computing anomaly
+            for it in range(Nsteps):
+                slicemat = []
+                for idim in range(NmatAdims):
+                    if idim == istepdim:
+                        slicemat.append(j+it*istep)
+                    else:
+                        slicemat.append(slice(0,matAshape[idim]))
+                newmat[tuple(slicemat)] = (matA[tuple(slicemat)] -                   \
+                  meannewmat[tuple(slicenewmat)])/stdnewmat[tuple(slicenewmat)]
+        additional = [meannewmat, stdnewmat]
+
     elif opn == 'inv':
         newmat = retype(1., valtype) / matA
Index: /trunk/tools/nc_var_tools.py
===================================================================
--- /trunk/tools/nc_var_tools.py	(revision 1560)
+++ /trunk/tools/nc_var_tools.py	(revision 1561)
@@ -813,5 +813,5 @@
     Nnewdims = len(newdims)
     
-    if Nnewdims == 0:
+    if Nnewdims == 0 or len(newvar.shape) == 0:
         print '%2s %f' % ( 'NC',newvar )
     elif Nnewdims == 1:
@@ -7206,4 +7206,23 @@
              'divc',[modval1]: [prevalues] divide by [modval1]
              'forwrdderiv',[N],[ord],[dim]: un-scaled forward [N]-derivative of order [ord] along dimension [dim] of [var]
+             'ifreq_anom',[stepdimn],[stepvardimn],[istep]: computing anomalies by substracting sub-means at each
+                 step by averaging from there all values taken every [istep] along dimension [stepdim]. 
+                 mean(j) = sum(matA[j+k*istep]_k=0,Nstep)/Nstep; Nstep = len(stepdimn)/istep; j=[0,istep]
+               [stepdimn]= name of the dimension along which to sample
+               [stepvardimn]= name of the variable-dimension with the values for [stepdimn]
+               [istep]= frequency to sample
+             'ifreq_mean',[stepdimn],[stepvardimn],[istep]: computing sub-means at each step by averaging from 
+                 there all values taken every [istep] along dimension [stepdim]. 
+                 mean(j) = sum(matA[j+k*istep]_k=0,Nstep)/Nstep; Nstep = len(stepdimn)/istep; j=[0,istep]
+               [stepdimn]= name of the dimension along which to sample
+               [stepvardimn]= name of the variable-dimension with the values for [stepdimn]
+               [istep]= frequency to sample
+             'ifreq_normmeanstd',[stepdimn],[stepvardimn],[istep]: nbormalizing anomalies by substracting mean(j)/dtsv(j) 
+                 sub-stats at each step by averaging from there all values taken every [istep] along dimension[stepdim]. 
+                 mean(j) = sum(matA[j+k*istep]_k=0,Nstep)/Nstep; std(j) = std(matA[j+k*istep]_k=0,Nstep); 
+                 Nstep = len(stepdimn)/istep; j=[0,istep]
+               [stepdimn]= name of the dimension along which to sample
+               [stepvardimn]= name of the variable-dimension with the values for [stepdimn]
+               [istep]= frequency to sample
              'inv': inverting [prevalues] (1/[prevalues])
              'lowthres',[modval1],[modval2]: if [prevalues] < [modval1]; prevalues = [modval2]
@@ -7250,5 +7269,10 @@
     opinione = ['mul', 'div', 'pot']
     # Operations which need to deal with dimension's names and indices
-    opdim=['norm_meanminmax', 'norm_meanstd', 'norm_minmax']
+    opdim = ['ifreq_anom', 'ifreq_mean', 'ifreq_normmeanstd', 'norm_meanminmax',     \
+      'norm_meanstd', 'norm_minmax']
+    # Operations which require a new dimension
+    opnewdim = ['ifreq_mean']
+    # istep Operations
+    opistep = ['ifreq_anom', 'ifreq_mean', 'ifreq_normmeanstd']
 
     varname = varinfo.split(',')[0]
@@ -7310,8 +7334,11 @@
     Dadditionalunits = {}
     # Dictionary with the Headers for the names of the additional variables
-    DHadditionalvars = {'norm_meanminmax': ['mean','min','max'],                     \
+    DHadditionalvars = {'ifreq_anom': ['ifreq_mean'], 'ifreq_normmeanstd':           \
+      ['ifreq_mean', 'ifreq_std'], 'norm_meanminmax': ['mean','min','max'],          \
       'norm_meanstd': ['mean','std'], 'norm_minmax': ['min','max']}
 
     prevalues = None
+    slicewrft = None
+
     for ifv in range(Nfileopers):
         operN = opers[ifv]
@@ -7429,5 +7456,5 @@
                     Ddv[dimn] = cftimes
                     Dtdv[dimn] = np.float
-                    Dadv[dimn] = {'stadard_name': 'time', 'long_name': 'time',       \
+                    Dadv[dimn] = {'standard_name': 'time', 'long_name': 'time',      \
                       'units': tunitsval + ' since ' + refdateS}
 
@@ -7484,7 +7511,91 @@
                 usedims = 'any'
             for iv in range(len(normi)): normi[iv] = str(normi[iv])
-            opern = opernS + ',' + ':'.join(dimvarvals) + ',' +                  \
+            opern = opernS + ',' + ':'.join(dimvarvals) + ',' +                      \
               operN.split(',')[1] + ',' + ':'.join(normi)
+
+            if gen.searchInlist(opistep, opernS):
+                opern = opernS + ',' + ':'.join(dimvarvals) + ',' +                  \
+                  operN.split(',')[3] + ',' + operN.split(',')[1]
+                usedims = list(dimvarvals)
+
+            if gen.searchInlist(opnewdim, opernS):
+                # Creation of a new dimension for the operation
+                dimv = list(dimvarvals)
+                idim = dimv.index(operN.split(',')[1])
+                newdimn = opernS + dimv[idim]
+                print infmsg
+                print '    ' + fname + ": operation '" + opernS + "' will create " + \
+                  " a new dimension '" + newdimn + "' !!"
+                dimv[idim] = newdimn
+                dimvarvals = tuple(dimv)
+                vardimnstep = operN.split(',')[2]
+                istep = int(operN.split(',')[3])
+                if vardimnstep != 'WRFtime':
+                    if not objnc.variables.has_key[vardimnstep]:
+                        print errormsg
+                        print '  ' + fname + ": file '" + filen + "' does not " +    \
+                          "have variable '" + vardimnstep + "' to be used for '" +   \
+                          "ifreq_mean' !!"
+                        print '    available ones:', objnc.variables.keys()
+                        quit(-1)
+                    odimv = objnc.variables[vardimnstep]
+                    varslice, Dddv[newdimn] = SliceVarDict(odimv, dictslice)
+                    allvalues = odimv[tuple(varslice)]
+                    Dtdv[newdimn] = odimv.dtype
+                    dicattrs = {}
+                    for attrn in odimv.ncattrs():
+                        attrval = odimv.getncattr(attrn)
+                        dicattrs[attrn] = attrval
+                    Dadv[newdimn] = dicattrs
+                else:
+                    if slicewrft is None:
+                        print infmsg
+                        print '  '+fname+ ": computing time-variable from 'WRF'-Time"
+                        datesv = []
+                        wrft = objnc.variables['Times']
+                        slicewrft, dwrfd = SliceVarDict(wrft, dictslice)
+                        fdates = wrft[tuple(slicewrft)]
+                        refdate = '19491201000000'
+                        tunitsval = 'minutes'
+                        yrref=refdate[0:4]
+                        monref=refdate[4:6]
+                        dayref=refdate[6:8]
+                        horref=refdate[8:10]
+                        minref=refdate[10:12]
+                        secref=refdate[12:14]
+        
+                        refdateS= yrref +'-'+ monref + '-' + dayref + ' ' + horref + \
+                          ':' + minref + ':' + secref
+        
+                        dt = fdates.shape[0]
+                        allvalues = np.zeros((dt), dtype=np.float)
+                        for it in range(dt):
+                            wrfdates = gen.datetimeStr_conversion(wrft[it,:],        \
+                             'WRFdatetime', 'matYmdHMS')
+                            allvalues[it] = gen.realdatetime1_CFcompilant(wrfdates,  \
+                              refdate, tunitsval)
+                        tunits = tunitsval + ' since ' + refdateS
+                    else:
+                        # Using previous computed WRFtime
+                        allvalues = cftimes
+
+                    Dtdv[newdimn] = np.float
+                    Dadv[newdimn] = {'standard_name': 'time', 'long_name': 'time',   \
+                      'units': tunitsval + ' since ' + refdateS, 'calendar':         \
+                      'standard'}
+
+                # Filling-up new dimension characteristics
+                Dds[newdimn] = istep
+
+                Dddv[newdimn] = [newdimn]
+                newdimvals = np.zeros((istep), dtype=Dtdv[newdimn])
+                dsize = allvalues.shape[0]
+                Nsteps = dsize/istep
+                for it in range(istep):
+                    newdimvals[it] = np.mean(allvalues[slice(it,istep*Nsteps,istep)])
+                Ddv[newdimn] = newdimvals[:]
+
             Dadditionaldim[opernS+'@'+varn] = usedims
+
             if varn != 'continuation':
                 Dadditionalunits[opernS+'@'+varn] = get_varunits(varobj)
@@ -7593,4 +7704,5 @@
                     newvn = Hvarn + '_' + operS + '_' + vnS
                     newvar = objofile.createVariable(newvn, 'f4', tuple(newvdims))
+                    print fname + ': LLuis shapes:', newvar.shape, 'newvdims:', newvdims, 'addvs:', addvs[iHv].shape
                     newvar[:] = addvs[iHv]
                     Lattr = Hvarn + ' value for operation ' + operS + ' with variable ' + vnS
