diff --git a/Common/CreateFile.cpp b/Common/CreateFile.cpp
index 367c6cc1cc1bc71c0f4712836c82bba0da982f30..b3f97f4ed35bb8e52269969ed31bf07d3db3d66e 100644
--- a/Common/CreateFile.cpp
+++ b/Common/CreateFile.cpp
@@ -48,6 +48,7 @@ int GuessFileFormatFromFileName(const char *name)
   else if(!strcmp(ext, ".med"))  return FORMAT_MED;
   else if(!strcmp(ext, ".mesh")) return FORMAT_MESH;
   else if(!strcmp(ext, ".bdf"))  return FORMAT_BDF;
+  else if(!strcmp(ext, ".diff")) return FORMAT_DIFF;
   else if(!strcmp(ext, ".nas"))  return FORMAT_BDF;
   else if(!strcmp(ext, ".p3d"))  return FORMAT_P3D;
   else if(!strcmp(ext, ".wrl"))  return FORMAT_VRML;
@@ -82,6 +83,7 @@ void GetDefaultFileName(int format, char *name)
   case FORMAT_MED:  strcpy(ext, ".med"); break;
   case FORMAT_MESH: strcpy(ext, ".mesh"); break;
   case FORMAT_BDF:  strcpy(ext, ".bdf"); break;
+  case FORMAT_DIFF: strcpy(ext, ".diff"); break;
   case FORMAT_P3D:  strcpy(ext, ".p3d"); break;
   case FORMAT_VRML: strcpy(ext, ".wrl"); break;
   case FORMAT_GIF:  strcpy(ext, ".gif"); break;
@@ -168,6 +170,11 @@ void CreateOutputFile(const char *filename, int format)
 				CTX.mesh.save_all, CTX.mesh.scaling_factor);
     break;
 
+  case FORMAT_DIFF:
+    GModel::current()->writeDIFF(name, CTX.mesh.binary, CTX.mesh.save_all,
+                                 CTX.mesh.scaling_factor);
+    break;
+
   case FORMAT_P3D:
     GModel::current()->writeP3D(name, CTX.mesh.save_all, CTX.mesh.scaling_factor);
     break;
diff --git a/Common/GmshDefines.h b/Common/GmshDefines.h
index 1f8031462f1f5d60aa534cfa537a15e4b05eda42..f8def9cb1662362781d6e2e519b72bda76a9a500 100644
--- a/Common/GmshDefines.h
+++ b/Common/GmshDefines.h
@@ -36,6 +36,7 @@
 #define FORMAT_BDF           31
 #define FORMAT_CGNS          32
 #define FORMAT_MED           33
+#define FORMAT_DIFF          34
 
 // Element types in .msh file format
 #define MSH_LIN_2  1
diff --git a/Fltk/Callbacks.cpp b/Fltk/Callbacks.cpp
index dedc32b4465b15d594a1e2fcc227e4b2f25d4a6e..f459f2de3d148d4b0b009084427449a663431332 100644
--- a/Fltk/Callbacks.cpp
+++ b/Fltk/Callbacks.cpp
@@ -651,6 +651,7 @@ int _save_geo(const char *name){ return geo_dialog(name); }
 int _save_cgns(const char *name){ return cgns_write_dialog(name); }
 int _save_unv(const char *name){ return unv_dialog(name); }
 int _save_vtk(const char *name){ return generic_mesh_dialog(name, "VTK Options", FORMAT_VTK, true); }
+int _save_diff(const char *name){ return generic_mesh_dialog(name, "DIFFPACK Options", FORMAT_DIFF, true); }
 int _save_med(const char *name){ return generic_mesh_dialog(name, "MED Options", FORMAT_MED, false); }
 int _save_mesh(const char *name){ return generic_mesh_dialog(name, "MESH Options", FORMAT_MESH, false); }
 int _save_bdf(const char *name){ return bdf_dialog(name); }
@@ -681,6 +682,7 @@ int _save_auto(const char *name)
   case FORMAT_MED  : return _save_med(name);
   case FORMAT_MESH : return _save_mesh(name);
   case FORMAT_BDF  : return _save_bdf(name);
+  case FORMAT_DIFF : return _save_diff(name);
   case FORMAT_P3D  : return _save_p3d(name);
   case FORMAT_STL  : return _save_stl(name);
   case FORMAT_VRML : return _save_vrml(name);
@@ -725,6 +727,7 @@ void file_save_as_cb(CALLBACK_ARGS)
     {"CGNS" TT "*.cgns", _save_cgns},
 #endif
     {"I-deas universal mesh" TT "*.unv", _save_unv},
+    {"DIFFPACK 3D mesh" TT "*.diff", _save_diff},
     {"VTK mesh" TT "*.vtk", _save_vtk},
 #if defined(HAVE_MED)
     {"MED file" TT "*.med", _save_med},
diff --git a/Geo/GModel.h b/Geo/GModel.h
index df1c3860eab179140aa989ec45f872dbe5bce158..97ea56612a889b4b318006053771dda62e486ebb 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -327,6 +327,10 @@ class GModel
   int writeVTK(const std::string &name, bool binary=false,
                bool saveAll=false, double scalingFactor=1.0,
 	       bool bigEndian=false);
+
+  // DIFFPACK format
+  int writeDIFF(const std::string &name, bool binary=false,
+               bool saveAll=false, double scalingFactor=1.0);
 };
 
 #endif
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index 264254cdf39fa7baf3c9f3ffccf0338e6c49399c..04480a73410a128e625836ea25104622a2a81654 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -2153,3 +2153,113 @@ int GModel::readVTK(const std::string &name, bool bigEndian)
   fclose(fp);
   return 1;
 }
+
+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;
+
+  FILE *fp = fopen(name.c_str(), binary ? "wb" : "w");
+  if(!fp){
+    Msg::Error("Unable to open file '%s'", name.c_str());
+    return 0;
+  }
+  
+  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
+  int numVertices = indexMeshVertices(saveAll);
+
+  // get all the entities in the model
+  std::vector<GEntity*> entities;
+  getEntities(entities);
+
+  // loop over all elements we need to save
+  int numElements = 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())
+          numElements++;
+      }
+    }
+  }
+  fprintf(fp, "\n\n");
+  fprintf(fp, " Finite element mesh (GridFE):\n\n");
+  fprintf(fp, " Number of space dim. =   3\n");
+  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, " Only one subdomain el              : dpFALSE\n");
+  fprintf(fp, " Lattice data                     ? 0\n\n\n\n");
+  int nbi = getNumFaces();
+  fprintf(fp, " %d Boundary indicators:  ", getNumFaces());
+  for(fiter it = firstFace(); it != lastFace(); ++it){
+    if(saveAll || (*it)->physicals.size()){
+      fprintf(fp, " %d", (*it)->tag());
+    }
+  } 
+  fprintf(fp, "\n\n\n");
+  fprintf(fp,"  Nodal coordinates and nodal boundary indicators,\n");
+  fprintf(fp,"  the columns contain:\n");
+  fprintf(fp,"   - node number\n");
+  fprintf(fp,"   - coordinates\n");
+  fprintf(fp,"   - no of boundary indicators that are set (ON)\n");
+  fprintf(fp,"   - the boundary indicators that are set (ON) if any.\n");
+  fprintf(fp,"#\n");
+
+  // write mesh vertices
+  fprintf(fp, "%d\n", numVertices);
+  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);
+      for(unsigned int k = 0; k < entities[i]->physicals.size(); k++){
+        fprintf(fp," %d %d\n", entities[i]->tag(), entities[i]->physicals[k]);
+      }
+    }
+  /* 
+  for(viter it = firstVertex(); it != lastVertex(); ++it){
+    for(unsigned int j = 0; j < (*it)->mesh_vertices.size(); j++) { 
+      (*it)->mesh_vertices[j]->writeDIFF(fp, binary, scalingFactor);
+      fprintf(fp," %d\n",(*it)->physicals.size());
+    }
+  }
+  */
+  fprintf(fp, "\n");
+  
+  fprintf(fp, "\n");
+  fprintf(fp,     "  Element types and connectivity\n");
+  fprintf(fp,     "  the columns contain:\n");
+  fprintf(fp,     "   - element number\n");
+  fprintf(fp,     "   - element type\n");
+  fprintf(fp,     "   - subdomain number \n");
+  fprintf(fp,     "   - the global node numbers of the nodes in the element.\n");
+  fprintf(fp,     "#\n");
+
+  static int element_number = 0;
+  fprintf(fp, "%d\n", numElements);
+  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){
+          element_number++;
+          int physical_property = 1;
+          entities[i]->getMeshElement(j)->writeDIFF(fp, binary, physical_property);
+        }
+      }
+    }
+  }
+  fprintf(fp, "\n");
+  
+  fclose(fp);
+  return 1;
+}
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index faff1b1e966a71213ab8b27848d106babac69675..ffb98b2a5d2d0851842391ca20eac3769101f8c0 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -618,6 +618,76 @@ void MElement::writeBDF(FILE *fp, int format, int elementary)
   }
 }
 
+void MElement::writeDIFF(FILE *fp, bool binary, int physical_property)
+{
+  /*  nbn = 4;
+    if(s->VSUP) {
+      type = TETRAHEDRON_2;
+      nbs = 6;
+      if(s->Volume_Simplexe() < 0) {
+        Vertex *temp;
+        temp = s->V[0]; s->V[0] = s->V[1]; s->V[1] = temp;
+        temp = s->VSUP[1]; s->VSUP[1] = s->VSUP[2]; s->VSUP[2] = temp;
+        temp = s->VSUP[5]; s->VSUP[5] = s->VSUP[3]; s->VSUP[3] = temp;
+      }
+    }
+    else{
+      type = TETRAHEDRON;
+      if(s->Volume_Simplexe() < 0) {
+        Vertex *temp;
+        temp = s->V[0];
+        s->V[0] = s->V[1];
+        s->V[1] = temp;
+      }
+    }
+  fprintf(MSHFILE, "%d %d 2 %d %d",
+            MSH_ELEMENT_NUM++, type, MSH_PHYSICAL_NUM ? MSH_PHYSICAL_NUM : s->iEnt,
+            s->iEnt);
+
+  if(DIFF_PHYSICAL_ORI > 0) {
+    for(i = 0; i < nbn; i++)
+      fprintf(DIFFFILE, " %d", s->V[i]->Num);
+    for(i = 0; i < nbs; i++)
+      fprintf(DIFFFILE, " %d", s->VSUP[i]->Num);
+  }
+  else {
+    for(i = 0; i < nbn; i++)
+      fprintf(DIFFFILE, " %d", s->V[nbn - i - 1]->Num);
+    for(i = 0; i < nbs; i++)
+      fprintf(DIFFFILE, " %d", s->VSUP[nbs - i - 1]->Num);
+  }
+  fprintf(DIFFFILE, "\n");*/
+
+  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();
+  if(binary){
+    int verts[60];
+    verts[0] = n;
+    for(int i = 0; i < n; i++)
+      verts[i + 1] = getVertexVTK(i)->getIndex() - 1;
+    fwrite(verts, sizeof(int), n + 1, fp);
+  }
+  else{
+    if(type == MSH_TET_10)
+      fprintf(fp, "%d %s", getNum() - start, "ElmT10n3D ");
+    else
+      fprintf(fp, "%d %s", getNum() - start, "ElmT4n3D ");
+    fprintf(fp, " %d ", physical_property);
+    for(int i = 0; i < n; i++)
+      fprintf(fp, " %d", getVertexVTK(i)->getIndex() - 1);
+    fprintf(fp, "\n");
+  }
+}
+
 int MElement::getInfoMSH(const int typeMSH, const char **const name)
 {
   switch(typeMSH){
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 3e0db7e9cc19f56be5fd028b4157ac8a1f4381da..9979f8eb997c307e37773d466b6eb3c3a0bede7a 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -224,12 +224,14 @@ 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);
+ 
   // info for specific IO formats (returning 0 means that the element
   // is not implemented in that format)
   virtual int getTypeForMSH() const { return 0; }
   virtual int getTypeForUNV() const { return 0; }
   virtual int getTypeForVTK() const { return 0; }
+  virtual int getTypeForDIFF() const { return 0; }
   virtual const char *getStringForPOS() const { return 0; }
   virtual const char *getStringForBDF() const { return 0; }
 
@@ -1172,6 +1174,7 @@ class MTetrahedron : public MElement {
   virtual int getTypeForMSH() const { return MSH_TET_4; }
   virtual int getTypeForUNV() const { return 111; } // solid linear tetrahedron
   virtual int getTypeForVTK() const { return 10; }
+  virtual int getTypeForDIFF() const { return MSH_TET_4; }
   virtual const char *getStringForPOS() const { return "SS"; }
   virtual const char *getStringForBDF() const { return "CTETRA"; }
   virtual void revert()
@@ -1334,6 +1337,7 @@ class MTetrahedron10 : public MTetrahedron {
   virtual int getTypeForMSH() const { return MSH_TET_10; }
   virtual int getTypeForUNV() const { return 118; } // solid parabolic tetrahedron
   //virtual int getTypeForVTK() const { return 24; }
+  virtual int getTypeForDIFF() const { return MSH_TET_10; }
   virtual const char *getStringForPOS() const { return "SS2"; }
   virtual const char *getStringForBDF() const { return "CTETRA"; }
   virtual void revert()
diff --git a/Geo/MVertex.cpp b/Geo/MVertex.cpp
index d39a7eb15432933d69465806fb06622439f0f46f..2aaaa0fac7d4cae32a0322b915e95984b86fc413 100644
--- a/Geo/MVertex.cpp
+++ b/Geo/MVertex.cpp
@@ -141,6 +141,14 @@ void MVertex::writeBDF(FILE *fp, int format, double scalingFactor)
   }
 }
 
+void MVertex::writeDIFF(FILE *fp, bool binary, double scalingFactor)
+{
+  if(_index < 0) return; // negative index vertices are never saved
+
+  fprintf(fp, " %d ( %25.16E , %25.16E , %25.16E )\n",
+          getIndex(), x() * scalingFactor, y() * scalingFactor, z() * scalingFactor);
+}
+
 std::set<MVertex*, MVertexLessThanLexicographic>::iterator 
 MVertex::linearSearch(std::set<MVertex*, MVertexLessThanLexicographic> &pos)
 {
diff --git a/Geo/MVertex.h b/Geo/MVertex.h
index e2ef9cdc47755a6a000450bb552da923e32f4f7d..9c9f5be8acbb114b8d4d54cb009de17de63fc0c8 100644
--- a/Geo/MVertex.h
+++ b/Geo/MVertex.h
@@ -113,6 +113,7 @@ class MVertex{
 		bool bigEndian=false);
   void writeMESH(FILE *fp, double scalingFactor=1.0);
   void writeBDF(FILE *fp, int format=0, double scalingFactor=1.0);
+  void writeDIFF(FILE *fp, bool binary, double scalingFactor=1.0);
 };
 
 class MEdgeVertex : public MVertex{
diff --git a/doc/CREDITS.txt b/doc/CREDITS.txt
index 8829a70894f181813cfa1949aeec303e3f6a7f25..656b0c0efeba483117dffe4300803ea055ff8806 100644
--- a/doc/CREDITS.txt
+++ b/doc/CREDITS.txt
@@ -20,7 +20,7 @@ with the MacOS port and the tensor display code; Pierre Badel for help
 with the GSL integration; Marc Ume for the original list code; Matt
 Gundry for the Plot3d mesh format; Jozef Vesely for help with the
 Tetgen integration; Koen Hillewaert for high order element mappings
-and other improvements.
+and other improvements; Jacques Lechelle for the DIFFPACK mesh format.
 
 The AVL tree code (Common/avl.*) and the YUV image code
 (Graphics/gl2yuv.*) are copyright (C) 1988-1993, 1995 The Regents of