diff --git a/python/ibmppp/ibmppp.py b/python/ibmppp/ibmppp.py
index b634fc6..6fd5e65 100644
--- a/python/ibmppp/ibmppp.py
+++ b/python/ibmppp/ibmppp.py
@@ -621,6 +621,108 @@ class ibmppp:
         fid.write('</Xdmf> \n')
         fid.close()
 
+    def saveFieldSingleFile(self,filename,key,append=False,downsample=False,xdmf=False,verbose=False):
+      '''Writes field data into a single h5 file'''
+      # Determine filename of output file
+      file_out = filename+'.h5'
+      # Compute global nx,ny,nz
+      nxg = self.__procGrid[key][1][-1]-self.__procGrid[key][0][0]+1
+      nyg = self.__procGrid[key][3][-1]-self.__procGrid[key][2][0]+1
+      nzg = self.__procGrid[key][5][-1]-self.__procGrid[key][4][0]+1
+      # Adjust these values if downsampling is requested
+      if downsample:
+        nxg = np.ceil(nxg/2).astype('int')
+        nyg = np.ceil(nyg/2).astype('int')
+        nzg = np.ceil(nzg/2).astype('int')
+      # Setup the correct slicing
+      if not downsample:
+        ishift = 0
+        jshift = 0
+        kshift = 0
+        ifb = self.__localChunkBounds[key][0]
+        jfb = self.__localChunkBounds[key][2]
+        kfb = self.__localChunkBounds[key][4]
+        nxl = self.__localChunkSize[key][0]
+        nyl = self.__localChunkSize[key][1]
+        nzl = self.__localChunkSize[key][2]
+        step = 1
+      else:
+        ishift = 1 if self.__localChunkBounds[key][0]%2==0 else 0
+        jshift = 1 if self.__localChunkBounds[key][2]%2==0 else 0
+        kshift = 1 if self.__localChunkBounds[key][4]%2==0 else 0
+        ifb = np.ceil((self.__localChunkBounds[key][0]+ishift)/2).astype('int')
+        jfb = np.ceil((self.__localChunkBounds[key][2]+jshift)/2).astype('int')
+        kfb = np.ceil((self.__localChunkBounds[key][4]+kshift)/2).astype('int')
+        nxl = np.ceil(self.__localChunkSize[key][0]/2).astype('int')
+        nyl = np.ceil(self.__localChunkSize[key][1]/2).astype('int')
+        nzl = np.ceil(self.__localChunkSize[key][2]/2).astype('int')
+        step = 2
+      isl_ch = slice(self.__nghx+ishift,self.__nghx+self.__localChunkSize[key][0],step)
+      jsl_ch = slice(self.__nghy+jshift,self.__nghy+self.__localChunkSize[key][1],step)
+      ksl_ch = slice(self.__nghz+kshift,self.__nghz+self.__localChunkSize[key][2],step)
+      isl_h5 = slice(ifb-1,ifb+nxl-1)
+      jsl_h5 = slice(jfb-1,jfb+nyl-1)
+      ksl_h5 = slice(kfb-1,kfb+nzl-1)
+      # If append flag is set, we add a dataset to the existing file
+      if append:
+        ioflag = 'a'
+      else:
+        ioflag = 'w'
+      # I avoid parallel IO because h5py is not available everywhere with
+      # MPI support. So rank 0 does all the writing and we communicate manually.
+      if self.__rank==0:
+        if verbose:
+          print('[saveFieldSingleFile] rank {} / {} writing to file'.format(self.__rank,self.__nproc-1),end='\r')
+        fid = h5py.File(file_out,ioflag)
+        gid = fid.create_group('/'+key)
+        gid.create_dataset('nx',data=nxg)
+        gid.create_dataset('ny',data=nyg)
+        gid.create_dataset('nz',data=nzg)
+        did = gid.create_dataset('data',(nzg,nyg,nxg),dtype=self.__precision)
+        did[ksl_h5,jsl_h5,isl_h5] = np.transpose(self.field[key][isl_ch,jsl_ch,ksl_ch])
+        fid.close()
+        self.__comm.send(True, dest=self.__rank+1, tag=1)
+      else:
+        msg = self.__comm.recv(source=self.__rank-1,tag=1)
+        if verbose:
+          print('[saveFieldSingleFile] rank {} / {} writing to file'.format(self.__rank,self.__nproc-1),end='\r')
+        fid = h5py.File(file_out,'a')
+        fid['/'+key+'/data'][ksl_h5,jsl_h5,isl_h5] = np.transpose(self.field[key][isl_ch,jsl_ch,ksl_ch])
+        fid.close()
+        if self.__rank!=self.__nproc-1:
+          self.__comm.send(True, dest=self.__rank+1, tag=1)
+      # Create an XDMF file for the field
+      if xdmf and self.__rank==0:
+        if verbose:
+          print('[saveFieldSingleFile] generating XDMF',end='\r')
+        # Construct XDMF filename
+        file_xdmf = filename+'_'+key+'.xdmf'
+        filename_h5 = os.path.basename(file_out)
+        # Open file and write header
+        fid = open(file_xdmf,'w') 
+        fid.write('<?xml version="1.0"?>\n')
+        fid.write('<Xdmf Version="2.0">\n')
+        fid.write('   <Domain>\n')
+        fid.write('      <Grid Name="{}">\n'.format(key))
+        #fid.write('         <Time Value="%.2f" />\n',params.general.simtime)
+        dx = self.__dx[key]*step
+        x0 = self.grid[key][0][0]
+        y0 = self.grid[key][1][0]
+        z0 = self.grid[key][2][0]
+        fid.write('          <Topology TopologyType="3DCORECTMesh" NumberOfElements="{:d} {:d} {:d}"/>\n'.format(nzg,nyg,nxg))
+        fid.write('          <Geometry Origin="" Type="ORIGIN_DXDYDZ">\n')
+        fid.write('            <DataItem DataType="Float" Dimensions="3" Format="XML" Precision="8"> {:12f} {:12f} {:12f}</DataItem>\n'.format(z0,y0,x0))
+        fid.write('            <DataItem DataType="Float" Dimensions="3" Format="XML" Precision="8"> {:12f} {:12f} {:12f}</DataItem>\n'.format(dx,dx,dx))
+        fid.write('          </Geometry>\n')
+        fid.write('          <Attribute Name="{}" Dimensions="{:d} {:d} {:d}">\n'.format(key,nzg,nyg,nxg))
+        fid.write('          <DataItem Format="HDF5" DataType="Float" Dimensions="{:d} {:d} {:d}">\n'.format(nzg,nyg,nxg))
+        fid.write('          {}:/{}/data </DataItem>\n'.format(filename_h5,key))
+        fid.write('          </Attribute>\n')
+        fid.write('      </Grid>\n')
+        fid.write('   </Domain>\n')
+        fid.write('</Xdmf> \n')
+        fid.close()
+
     def saveStatistics(self,filename):
       '''Writes all gathered statistics to a h5 file.'''
       if self.__rank==0:
@@ -1201,17 +1303,52 @@ class ibmppp:
       '''Subtract mean from field.'''
       # Check if this operation can be executed
       if self.__rank==0:
-        if not key in spatialMean or spatialMean[key] is None:
+        if not key in self.spatialMean or self.spatialMean[key] is None:
           raise ValueError('Statistics not available')
       # Broadcast statistics
       self.spatialMean[key] = self.__comm.bcast(self.spatialMean[key],root=0)
       # Get dimensions
       dims = self.spatialMean[key].shape
-      # Slice chunk
-      if dims[0]==1:
-        slx = slice(0,1)
-      else:
-        slx = slice(self.ib)
+      # Slice statistics for each processor
+      ib,ie = (0,1) if dims[0]==1 else (self.__localChunkBounds[key][0]-1,self.__localChunkBounds[key][1])
+      jb,je = (0,1) if dims[1]==1 else (self.__localChunkBounds[key][2]-1,self.__localChunkBounds[key][3])
+      kb,ke = (0,1) if dims[2]==1 else (self.__localChunkBounds[key][4]-1,self.__localChunkBounds[key][5])
+      # Subtract mean from field
+      self.field[key][
+        self.__nghx:-self.__nghx,
+        self.__nghy:-self.__nghy,
+        self.__nghz:-self.__nghz] -= self.spatialMean[key][ib:ie,jb:je,kb:ke]
+      # Exchange ghost cells and set boundary conditions
+      self.exchangeGhostCells(key)
+      self.imposeBoundaryConditions(key)
+    
+    def normalizeByStandardDeviation(self,key):
+      '''Normalize field by its standard deviation.'''
+      # Check if this operation can be executed
+      if self.__rank==0:
+        if (not key in self.spatialMean or self.spatialMean[key] is None or
+            not key in self.spatialMsqr or self.spatialMsqr[key] is None):
+          raise ValueError('Statistics not available')
+      # Broadcast statistics
+      self.spatialMean[key] = self.__comm.bcast(self.spatialMean[key],root=0)
+      self.spatialMsqr[key] = self.__comm.bcast(self.spatialMsqr[key],root=0)
+      # Compute standard deviation
+      assert(self.spatialMean[key].shape==self.spatialMsqr[key].shape)
+      spatialStdDev = np.sqrt(self.spatialMsqr[key]-np.square(self.spatialMean[key]))
+      # Get dimensions
+      dims = spatialStdDev.shape
+      # Slice statistics for each processor
+      ib,ie = (0,1) if dims[0]==1 else (self.__localChunkBounds[key][0]-1,self.__localChunkBounds[key][1])
+      jb,je = (0,1) if dims[1]==1 else (self.__localChunkBounds[key][2]-1,self.__localChunkBounds[key][3])
+      kb,ke = (0,1) if dims[2]==1 else (self.__localChunkBounds[key][4]-1,self.__localChunkBounds[key][5])
+      # Subtract mean from field
+      self.field[key][
+        self.__nghx:-self.__nghx,
+        self.__nghy:-self.__nghy,
+        self.__nghz:-self.__nghz] /= spatialStdDev[ib:ie,jb:je,kb:ke]
+      # Exchange ghost cells and set boundary conditions
+      self.exchangeGhostCells(key)
+      self.imposeBoundaryConditions(key)
 
     def collocateVelocity(self):
       '''Moves velocity components onto the pressure grid.'''