diff --git a/tools/vlsvdiff.cpp b/tools/vlsvdiff.cpp
index ef1eb7a5e..c0433230c 100644
--- a/tools/vlsvdiff.cpp
+++ b/tools/vlsvdiff.cpp
@@ -421,14 +421,13 @@ bool convertMesh(vlsvinterface::Reader& vlsvReader,
       cerr << "ERROR, failed to get array info for '" << _varToExtract << "' at " << __FILE__ << " " << __LINE__ << endl;
       return false;
    }
+   
+   if (gridName==gridType::SpatialGrid) {
       char *variableBuffer = new char[variableVectorSize * variableDataSize];
       float *variablePtrFloat = reinterpret_cast<float *>(variableBuffer);
       double *variablePtrDouble = reinterpret_cast<double *>(variableBuffer);
       uint *variablePtrUint = reinterpret_cast<uint *>(variableBuffer);
       int *variablePtrInt = reinterpret_cast<int *>(variableBuffer);
-
-
-   if (gridName==gridType::SpatialGrid){
   
       // Read the mesh array one node (of a spatial cell) at a time
       // and create a map which contains each cell's CellID and variable to be extracted
@@ -488,32 +487,49 @@ bool convertMesh(vlsvinterface::Reader& vlsvReader,
          if (storeCellOrder == true) {
             cellOrder[CellID] = i;
          }
-       }
-  
-   }else if (gridName==gridType::fsgrid){
- 
- 
+      }
+      delete variableBuffer;
+   } else if (gridName==gridType::fsgrid) {
+      
+
+      // Get Spatial Grid's  max refinement Level
+      int maxRefLevel=0;
+      list<pair<string, string>> meshAttributesIn;
+      meshAttributesIn.push_back(make_pair("name", "SpatialGrid"));
+      map<string,string> meshAttributesOut;
+      if (vlsvReader.getArrayAttributes("MESH", meshAttributesIn,meshAttributesOut) == false)
+      {
+         cerr << "ERROR, failed to get array info for '" << _varToExtract << "' at " << __FILE__ << " " << __LINE__ << endl;
+         return false;
+      }
+
+      std::map<string, string>::iterator attributesOutIt;
+      attributesOutIt = meshAttributesOut.find("max_refinement_level");
+      if (attributesOutIt != meshAttributesOut.end())
+      {
+         maxRefLevel = stoi(attributesOutIt->second);
+      }
       int numtasks;
-      int xcells,ycells,zcells; 
+      int xcells,ycells,zcells;
       vlsvReader.readParameter("numWritingRanks",numtasks);
       vlsvReader.readParameter("xcells_ini",xcells);
       vlsvReader.readParameter("ycells_ini",ycells);
       vlsvReader.readParameter("zcells_ini",zcells);
+      xcells*=pow(2,maxRefLevel);
+      ycells*=pow(2,maxRefLevel);
+      zcells*=pow(2,maxRefLevel);
       std::array<int,3> GlobalBox={xcells,ycells,zcells};
       std::array<int,3> thisDomainDecomp;
       
       //Compute Domain Decomposition Scheme for this vlsv file
       computeDomainDecomposition(GlobalBox,numtasks,thisDomainDecomp);
 
-
       std::array<int32_t,3> taskSize,taskStart;
       std::array<int32_t,3> taskEnd;
       int readOffset=0;
-      size_t readSize;
       int index,my_x,my_y,my_z;
       orderedData->clear();
 
-      //Read into buffer
       for (int task=0; task<numtasks; task++){
 
          my_x=task/thisDomainDecomp[2]/thisDomainDecomp[1];
@@ -533,54 +549,109 @@ bool convertMesh(vlsvinterface::Reader& vlsvReader,
          taskEnd[1]= taskStart[1]+taskSize[1];
          taskEnd[2]= taskStart[2]+taskSize[2];
          
-         readSize= taskSize[0] * taskSize[1] * taskSize[2];
-         std::vector<Real> readIn(variableVectorSize * variableDataSize*readSize);
-
-         
-         int counter2=0;
-         uint64_t globalindex;
-         int64_t counter=0;
-         for(int z=taskStart[2]; z<taskEnd[2]; z++) {
-            for(int y=taskStart[1]; y<taskEnd[1]; y++) {
-               for(int x=taskStart[0]; x<taskEnd[0]; x++) {
-
-                  //Get global index
-                  globalindex= x + y*xcells + z*xcells*ycells;
-
-                  if (vlsvReader.readArray("VARIABLE", variableAttributes, readOffset+counter,1, variableBuffer) == false) {
-                     cerr << "ERROR, failed to read variable '" << _varToExtract << "' at " << __FILE__ << " " << __LINE__ << endl;
-                     variableSuccess = false; 
-                     abort();
-                     break;
-                  }
-
-                  // Get the variable value
-                  Real extract = NAN;
+         int64_t readSize=  taskSize[0] * taskSize[1] * taskSize[2] ;
+         //Allocate vector for reading
+         std::vector<Real> buffer(readSize*variableVectorSize);
+
+         if ( variableDataSize==sizeof(Real)){
+            if (vlsvReader.readArray("VARIABLE", variableAttributes, readOffset, readSize,  (char*)buffer.data()) == false) {
+               cerr << "ERROR, failed to read variable '" << _varToExtract << "' at " << __FILE__ << " " << __LINE__ << endl;
+               variableSuccess = false; 
+               break;
+            }
+         }else{
+            std::vector<float> tmpbuffer(readSize * variableVectorSize);
+            if (vlsvReader.readArray("VARIABLE", variableAttributes, readOffset, readSize, (char *)tmpbuffer.data()) == false){
+               cerr << "ERROR, failed to read variable '" << _varToExtract << "' at " << __FILE__ << " " << __LINE__ << endl;
+               variableSuccess = false;
+               break;
+            }
+            for (int i = 0; i < readSize * variableVectorSize; i++){
+               buffer[i] = tmpbuffer[i];
+            }
+         }
 
-                  switch (variableDataType) {
+         uint64_t globalindex,counter=0;;
+         for (int z=taskStart[2]; z<taskEnd[2]; z++){
+               for (int y=taskStart[1]; y< taskEnd[1]; y++){
+                  for (int x=taskStart[0]; x<taskEnd[0]; x++){
+                     globalindex= x + y*xcells + z*xcells*ycells;
+                     Real data;
+                     switch (variableDataType){
                      case datatype::type::FLOAT:
-                        if(variableDataSize == sizeof(float)) extract = (Real)(variablePtrFloat[compToExtract]);
-                        if(variableDataSize == sizeof(double)) extract = (Real)(variablePtrDouble[compToExtract]);
+                        if (variableDataSize == sizeof(float))
+                           memcpy(&data, &buffer[counter + compToExtract], sizeof(float));
+                        if (variableDataSize == sizeof(double))
+                           memcpy(&data, &buffer[counter + compToExtract], sizeof(double));
                         break;
                      case datatype::type::UINT:
-                        extract = (Real)(variablePtrUint[compToExtract]);
+                        memcpy(&data, &buffer[counter + compToExtract], sizeof(uint));
                         break;
                      case datatype::type::INT:
-                        extract = (Real)(variablePtrInt[compToExtract]);
+                        memcpy(&data, &buffer[counter + compToExtract], sizeof(int));
                         break;
                      case datatype::type::UNKNOWN:
                         cerr << "ERROR, BAD DATATYPE AT " << __FILE__ << " " << __LINE__ << endl;
                         break;
+                     }
+                     //Add to map
+                     orderedData->insert(pair<uint64_t, Real>(globalindex, data));
+                     counter+=variableVectorSize;
                   }
-                  orderedData->insert(pair<uint64_t, Real>(globalindex, extract));
-                  counter++;
-               
                }
             }
-         }
          readOffset+=readSize;
       }
-   }else{
+   } else if (gridName==gridType::ionosphere) {
+      char *variableBuffer = new char[variableArraySize * variableVectorSize * variableDataSize];
+      float *variablePtrFloat = reinterpret_cast<float *>(variableBuffer);
+      double *variablePtrDouble = reinterpret_cast<double *>(variableBuffer);
+      uint *variablePtrUint = reinterpret_cast<uint *>(variableBuffer);
+      int *variablePtrInt = reinterpret_cast<int *>(variableBuffer);
+      
+      if (compToExtract + 1 > variableVectorSize) {
+         cerr << "ERROR invalid component, this variable has size " << variableVectorSize << endl;
+         abort();
+      }
+      
+      if (storeCellOrder == true) {
+         cellOrder.clear();
+      }
+      
+      orderedData->clear();
+      
+      if (vlsvReader.readArray("VARIABLE", variableAttributes, 0, variableArraySize, variableBuffer) == false) {
+         cerr << "ERROR, failed to read variable '" << _varToExtract << "' at " << __FILE__ << " " << __LINE__ << endl;
+         variableSuccess = false;
+      }
+      
+      for (uint64_t i=0; i<variableArraySize; ++i) {
+         // Get the variable value
+         Real extract = NAN;
+         
+         switch (variableDataType) {
+            case datatype::type::FLOAT:
+               if(variableDataSize == sizeof(float)) extract = (Real)(variablePtrFloat[i*variableVectorSize + compToExtract]);
+               if(variableDataSize == sizeof(double)) extract = (Real)(variablePtrDouble[i*variableVectorSize + compToExtract]);
+               break;
+            case datatype::type::UINT:
+               extract = (Real)(variablePtrUint[i*variableVectorSize + compToExtract]);
+               break;
+            case datatype::type::INT:
+               extract = (Real)(variablePtrInt[i*variableVectorSize + compToExtract]);
+               break;
+            case datatype::type::UNKNOWN:
+               cerr << "ERROR, BAD DATATYPE AT " << __FILE__ << " " << __LINE__ << endl;
+               break;
+         }
+         // Put those into the map
+         orderedData->insert(pair<uint64_t, Real>(i, extract));
+         if (storeCellOrder == true) {
+            cellOrder[i] = i;
+         }
+      }
+      delete variableBuffer;
+   } else {
     cerr<<"meshName not recognized\t" << __FILE__ << " " << __LINE__ <<endl;
     abort();
    }
@@ -591,7 +662,6 @@ bool convertMesh(vlsvinterface::Reader& vlsvReader,
    if (variableSuccess == false) {
       cerr << "ERROR reading array VARIABLE " << varToExtract << endl;
    }
-   delete variableBuffer;
    return meshSuccess && variableSuccess;
 }
 
@@ -737,14 +807,12 @@ bool pDistance(const map<uint, Real>& orderedData1,
             value = abs(it1->second - it2->second);
             *absolute = max(*absolute, value);
             length    = max(length, abs(it1->second));
-         
-            }
-         if (gridName==gridType::SpatialGrid){  
+         }
+         if (gridName==gridType::SpatialGrid) {
             array[cellOrder.at(it1->first)] = value;
-         }else if (gridName==gridType::fsgrid) {   
-
+         } else if (gridName==gridType::fsgrid || gridName==gridType::ionosphere) {
             array.at(it1->first)=value;
-            }  
+         }
       }
    } else if (p == 1) {
       for (map<uint,Real>::const_iterator it1=orderedData1.begin(); it1!=orderedData1.end(); ++it1) {
@@ -754,13 +822,12 @@ bool pDistance(const map<uint, Real>& orderedData1,
             value = abs(it1->second - it2->second);
             *absolute += value;
             length    += abs(it1->second);
-         
-            }
-         if (gridName==gridType::SpatialGrid){  
+         }
+         if (gridName==gridType::SpatialGrid) {
             array[cellOrder.at(it1->first)] = value;
-         }else if (gridName==gridType::fsgrid){   
+         } else if (gridName==gridType::fsgrid || gridName==gridType::ionosphere) {
             array[it1->first]=value;
-            }  
+         }
       }
    } else {
       for (map<uint,Real>::const_iterator it1=orderedData1.begin(); it1!=orderedData1.end(); ++it1) {
@@ -770,13 +837,12 @@ bool pDistance(const map<uint, Real>& orderedData1,
             value = pow(abs(it1->second - it2->second), p);
             *absolute += value;
             length    += pow(abs(it1->second), p);
-         
-            }
-         if (gridName==gridType::SpatialGrid){  
+         }
+         if (gridName==gridType::SpatialGrid) {
             array[cellOrder.at(it1->first)] = pow(value,1.0/p);
-         }else if (gridName==gridType::fsgrid){   
+         }else if (gridName==gridType::fsgrid || gridName==gridType::ionosphere) {
             array[it1->first]=pow(value,1.0/p);
-            }  
+         }
       }
       *absolute = pow(*absolute, 1.0 / p);
       length = pow(length, 1.0 / p);
@@ -1788,13 +1854,13 @@ int main(int argn,char* args[]) {
    
 
    //Figure out Meshname
-   if (attributes["--meshname"] == "SpatialGrid") { 
-      gridName=gridType::SpatialGrid ;
-   }else if (attributes["--meshname"]=="fsgrid"){
-      gridName=gridType::fsgrid ;
-   }else if (attributes["--meshname"]=="ionosphere"){
-      gridName=gridType::ionosphere ;
-   }else{
+   if (attributes["--meshname"] == "SpatialGrid") {
+      gridName=gridType::SpatialGrid;
+   } else if (attributes["--meshname"]=="fsgrid") {
+      gridName=gridType::fsgrid;
+   } else if (attributes["--meshname"]=="ionosphere") {
+      gridName=gridType::ionosphere;
+   } else {
       std::cout<<attributes["--meshname"]<<std::endl;
       std::cerr<<"Wrong grid type"<<std::endl;
       abort();