diff --git a/Geo/GModelIO_Geo.cpp b/Geo/GModelIO_Geo.cpp
index 206bd826d389cc5036afbd9fd22ebcc5384c5fd6..f37c46e45c88f9e5058491fa42d84b1f23c2d17b 100644
--- a/Geo/GModelIO_Geo.cpp
+++ b/Geo/GModelIO_Geo.cpp
@@ -3,6 +3,7 @@
 // See the LICENSE.txt file for license information. Please report all
 // bugs and problems to <gmsh@geuz.org>.
 
+#include <sstream>
 #include <stdlib.h>
 #include "GmshConfig.h"
 #include "GmshMessage.h"
@@ -31,7 +32,7 @@ void GModel::_createGEOInternals()
 
 void GModel::_deleteGEOInternals()
 {
-  delete _geo_internals;
+  if(_geo_internals) delete _geo_internals;
   _geo_internals = 0;
 }
 
@@ -43,6 +44,7 @@ int GModel::readGEO(const std::string &name)
 
 int GModel::exportDiscreteGEOInternals()
 {
+  if(_geo_internals) delete _geo_internals;
   _geo_internals = new GEO_Internals;
 
   for(viter it = firstVertex(); it != lastVertex(); it++){
@@ -329,8 +331,21 @@ int GModel::writeGEO(const std::string &name, bool printLabels)
 {
   FILE *fp = fopen(name.c_str(), "w");
 
-  for(viter it = firstVertex(); it != lastVertex(); it++)
-    (*it)->writeGEO(fp);
+  std::map<double, std::string> meshSizeParameters;
+  int cpt = 0;
+  for(viter it = firstVertex(); it != lastVertex(); it++){
+    double val = (*it)->prescribedMeshSizeAtVertex();
+    if(meshSizeParameters.find(val) == meshSizeParameters.end()){
+      std::ostringstream paramName;
+      paramName << "cl" << ++cpt;
+      fprintf(fp, "%s = %.16g;\n", paramName.str().c_str(),val);
+      meshSizeParameters.insert(std::make_pair(val, paramName.str()));
+    }
+  }
+  for(viter it = firstVertex(); it != lastVertex(); it++){
+    double val = (*it)->prescribedMeshSizeAtVertex();
+    (*it)->writeGEO(fp, meshSizeParameters[val]);
+  }
   for(eiter it = firstEdge(); it != lastEdge(); it++)
     (*it)->writeGEO(fp);
   for(fiter it = firstFace(); it != lastFace(); it++)
diff --git a/Geo/GVertex.cpp b/Geo/GVertex.cpp
index decf099956f5ec796433e53ff2b66289328b1d70..b5320c0c3251ee69a47a307634b1357c246cfc2f 100644
--- a/Geo/GVertex.cpp
+++ b/Geo/GVertex.cpp
@@ -60,10 +60,17 @@ std::string GVertex::getAdditionalInfoString()
   return sstream.str();
 }
 
-void GVertex::writeGEO(FILE *fp)
+void GVertex::writeGEO(FILE *fp, const std::string &meshSizeParameter)
 {
-  fprintf(fp, "Point(%d) = {%.16g, %.16g, %.16g, %.16g};\n",
-          tag(), x(), y(), z(), prescribedMeshSizeAtVertex());
+  if(meshSizeParameter.size())
+    fprintf(fp, "Point(%d) = {%.16g, %.16g, %.16g, %s};\n",
+            tag(), x(), y(), z(), meshSizeParameter.c_str());
+  else if(prescribedMeshSizeAtVertex() != MAX_LC)
+    fprintf(fp, "Point(%d) = {%.16g, %.16g, %.16g, %.16g};\n",
+            tag(), x(), y(), z(), prescribedMeshSizeAtVertex());
+  else
+    fprintf(fp, "Point(%d) = {%.16g, %.16g, %.16g};\n",
+            tag(), x(), y(), z());
 }
 
 unsigned int GVertex::getNumMeshElements()
diff --git a/Geo/GVertex.h b/Geo/GVertex.h
index 73a8c8911a5f32445d697d57bc115fb862b69874..3b6800dd207830888f5aa741435a79b76b300efa 100644
--- a/Geo/GVertex.h
+++ b/Geo/GVertex.h
@@ -66,7 +66,7 @@ class GVertex : public GEntity
   virtual std::string getAdditionalInfoString();
 
   // export in GEO format
-  virtual void writeGEO(FILE *fp);
+  virtual void writeGEO(FILE *fp, const std::string &meshSizeParameter="");
 
   // get number of elements in the mesh
   unsigned int getNumMeshElements();
diff --git a/Geo/Geo.cpp b/Geo/Geo.cpp
index 9267a08a1612b10cb10f007d8ea5d5a86ffce6c5..b1b41e874fd23792f71e00ae39987276c603a702 100644
--- a/Geo/Geo.cpp
+++ b/Geo/Geo.cpp
@@ -3389,7 +3389,6 @@ void setVolumeSurfaces(Volume *v, List_T *loops)
             List_Add(v->SurfacesByTag, &is);
           }
           else{
-	    printf("surface %d is not in GModel !\n", abs(is));
             Msg::Error("Unknown surface %d", is);
             return;
           }