diff --git a/Common/CommandLine.cpp b/Common/CommandLine.cpp
index 62e1b059f32350cb3e3702de2e67f580bc3ae06d..9993d8315864965f6881f47de46034026178d157 100644
--- a/Common/CommandLine.cpp
+++ b/Common/CommandLine.cpp
@@ -472,6 +472,8 @@ void Get_Options(int argc, char *argv[])
             CTX.mesh.format = FORMAT_P3D;
           else if(!strcmp(argv[i], "cgns"))
             CTX.mesh.format = FORMAT_CGNS;
+          else if(!strcmp(argv[i], "diff"))
+            CTX.mesh.format = FORMAT_DIFF;
           else if(!strcmp(argv[i], "med"))
             CTX.mesh.format = FORMAT_MED;
           else
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index 3010a9e2c291593d3608a8e7926bf943a3e690b9..5d9855290620b2017039124aef7627b457b81241 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -2157,11 +2157,10 @@ int GModel::readVTK(const std::string &name, bool bigEndian)
 int GModel::writeDIFF(const std::string &name, bool binary, bool saveAll,
                       double scalingFactor)
 {
-  Msg::Warning("DIFF export is experimental:");
-  Msg::Warning("*   only tetrahedra P1 or P2 are taken into account");
-  Msg::Warning("*   only ASCII output is implemented");
-
-  if(binary) return 0;
+  if(binary){
+    Msg::Error("Binary DIFF output is not implemented");
+    return 0;
+  }
 
   FILE *fp = fopen(name.c_str(), binary ? "wb" : "w");
   if(!fp){
@@ -2170,8 +2169,6 @@ int GModel::writeDIFF(const std::string &name, bool binary, bool saveAll,
   }
   
   if(noPhysicalGroups()) saveAll = true;
-  int nbP1 = 4; // number of P1 nodes
-  int nbP2 = 6; // number of P2 nodes which are not P1 (edges' midpoints)
 
   // get the number of vertices and index the vertices in a continuous
   // sequence
@@ -2182,12 +2179,15 @@ int GModel::writeDIFF(const std::string &name, bool binary, bool saveAll,
   getEntities(entities);
 
   // loop over all elements we need to save
-  int numElements = 0;
+  int numElements = 0, maxNumNodesPerElement = 0;
   for(unsigned int i = 0; i < entities.size(); i++){
     if(entities[i]->physicals.size() || saveAll){
       for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){
-        if(entities[i]->getMeshElement(j)->getTypeForDIFF())
+        MElement *e = entities[i]->getMeshElement(j);
+        if(e->getTypeForDIFF()){
           numElements++;
+          maxNumNodesPerElement = std::max(maxNumNodesPerElement, e->getNumVertices());
+        }
       }
     }
   }
@@ -2198,12 +2198,14 @@ int GModel::writeDIFF(const std::string &name, bool binary, bool saveAll,
   fprintf(fp, " Number of elements   =  %d\n", numElements);
   fprintf(fp, " Number of nodes      =  %d\n\n", numVertices);
   fprintf(fp, " All elements are of the same type : dpTRUE\n");
-  fprintf(fp, " Max number of nodes in an element: %d \n", nbP1 + nbP2);
+  fprintf(fp, " Max number of nodes in an element: %d \n", maxNumNodesPerElement);
   fprintf(fp, " Only one subdomain el              : dpFALSE\n");
   fprintf(fp, " Lattice data                     ? 0\n\n\n\n");
   int nbi = getNumFaces();
   fprintf(fp, " %d Boundary indicators:  ", nbi);
   for(fiter it = firstFace(); it != lastFace(); ++it){
+    // Jacques: are you sure about this? "nbi" might not reflect the
+    // number of tags written if saveAll==false...
     if(saveAll || (*it)->physicals.size()){
       fprintf(fp, " %d", (*it)->tag());
     }
@@ -2221,11 +2223,11 @@ int GModel::writeDIFF(const std::string &name, bool binary, bool saveAll,
   for(unsigned int i = 0; i < entities.size(); i++)
     for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++){
       entities[i]->mesh_vertices[j]->writeDIFF(fp, binary, scalingFactor);
-      fprintf(fp," [%d] ",entities[i]->faces().size());
+      fprintf(fp, " [%d] ", entities[i]->faces().size());
       std::list<GFace*> lf = entities[i]->faces();
       std::list<GFace*>::iterator it = lf.begin();
       for(unsigned int k = 0; k < lf.size(); k++){
-        fprintf(fp," %d ",(*it)->tag());
+        fprintf(fp," %d ", (*it)->tag());
         it++;
       }
       fprintf(fp,"\n");
@@ -2242,12 +2244,13 @@ int GModel::writeDIFF(const std::string &name, bool binary, bool saveAll,
   fprintf(fp,     "#\n");
 
   // write mesh elements
+  int num = 0;
   for(unsigned int i = 0; i < entities.size(); i++){
     if(entities[i]->physicals.size() || saveAll){
       for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){
-        int type = entities[i]->getMeshElement(j)->getTypeForDIFF();
-	if(type)
-          entities[i]->getMeshElement(j)->writeDIFF(fp, binary, entities[i]->tag());
+        MElement *e = entities[i]->getMeshElement(j);
+        if(e->getTypeForDIFF())
+          e->writeDIFF(fp, ++num, binary, entities[i]->tag());
       }
     }
   }
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index a2775ae2939b11e73a8bed0012ff3d859ca6acdd..f3b4df472a8cda77fe1c8e5faa401d4f4b4a3e93 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -617,16 +617,10 @@ void MElement::writeBDF(FILE *fp, int format, int elementary)
   }
 }
 
-void MElement::writeDIFF(FILE *fp, bool binary, int physical_property)
+void MElement::writeDIFF(FILE *fp, int num, bool binary, int physical_property)
 {
   int type = getTypeForDIFF();
   if(!type) return;
-  static int first = 1;
-  static int start = 0;
-  if(first){
-    start = getNum() - 1;
-    first = 0;
-  }
 
   setVolumePositive();
   int n = getNumVertices();
@@ -639,9 +633,9 @@ void MElement::writeDIFF(FILE *fp, bool binary, int physical_property)
   }
   else{
     if(type == MSH_TET_10)
-      fprintf(fp, "%d %s", getNum() - start, "ElmT10n3D ");
+      fprintf(fp, "%d %s", num, "ElmT10n3D ");
     else
-      fprintf(fp, "%d %s", getNum() - start, "ElmT4n3D ");
+      fprintf(fp, "%d %s", num, "ElmT4n3D ");
     fprintf(fp, " %d ", physical_property);
     for(int i = 0; i < n; i++)
       fprintf(fp, " %d", getVertexVTK(i)->getIndex() - 1);
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 9b07bc1d47aa89f5fde290e5301708f4c5c6c04c..b5b6cde1556ff6f4cbcd732f3885f907b8578b3c 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -224,7 +224,7 @@ class MElement
   virtual void writeVTK(FILE *fp, bool binary=false, bool bigEndian=false);
   virtual void writeMESH(FILE *fp, int elementary=1);
   virtual void writeBDF(FILE *fp, int format=0, int elementary=1);
-  virtual void writeDIFF(FILE *fp, bool binary=false, int physical_property=1);
+  virtual void writeDIFF(FILE *fp, int num, bool binary=false, int physical_property=1);
  
   // info for specific IO formats (returning 0 means that the element
   // is not implemented in that format)