diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 8a177d58030fd74da7ab5f2ec23c9fb1d7ceda10..3164e16666383224d233dbfadefe0a71ad02da52 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -1,4 +1,4 @@
-// $Id: GModel.cpp,v 1.80 2008-03-29 21:36:29 geuzaine Exp $
+// $Id: GModel.cpp,v 1.81 2008-03-30 17:45:11 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -52,9 +52,6 @@ GModel::GModel(std::string name)
 #if !defined(HAVE_GMSH_EMBEDDED)
   _fields = new FieldManager();
 #endif
-
-  //printf("sizeof(MVertex) = %d\n", sizeof(MVertex));
-  //printf("sizeof(MElement) = %d\n", sizeof(MElement));
 }
 
 GModel::~GModel()
@@ -77,6 +74,13 @@ GModel *GModel::current()
   return list.back();
 }
 
+GModel *GModel::findByName(std::string name)
+{
+  for(unsigned int i = 0; i < list.size(); i++)
+    if(list[i]->getName() == name) return list[i];
+  return 0;
+}
+
 void GModel::destroy()
 {
   for(riter it = firstRegion(); it != lastRegion(); ++it)
diff --git a/Geo/GModel.h b/Geo/GModel.h
index 1add954174032b201cdf25d32ca54360024b2870..ba187d520bb692855b5ef145ef2df243a0bee3f2 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -80,6 +80,9 @@ class GModel
   // returns the current model
   static GModel *current();
 
+  // returns the model of name
+  static GModel *findByName(std::string name);
+
   // Deletes everything in a GModel
   void destroy();
 
@@ -254,8 +257,10 @@ class GModel
   int readCGNS(const std::string &name);
   int writeCGNS(const std::string &name, double scalingFactor=1.0);
 
-  // Med interface ("Modele d'Echange de Donnees")
-  int readMED(const std::string &name);
+  // Med "Modele d'Echange de Donnees" file format (the static routine
+  // is allowed to load multiple models/meshes)
+  static int readMED(const std::string &name);
+  int readMED(const std::string &name, int meshIndex);
   int writeMED(const std::string &name, 
 	       bool saveAll=false, double scalingFactor=1.0);
 };
diff --git a/Geo/GModelIO_MED.cpp b/Geo/GModelIO_MED.cpp
index 9d25a2058b4e82e28d7e347f2291978130deec63..2a067ad69c0dd67eec0f3ea6f7eca0d47bf92f83 100644
--- a/Geo/GModelIO_MED.cpp
+++ b/Geo/GModelIO_MED.cpp
@@ -1,4 +1,4 @@
-// $Id: GModelIO_MED.cpp,v 1.20 2008-03-29 21:36:29 geuzaine Exp $
+// $Id: GModelIO_MED.cpp,v 1.21 2008-03-30 17:45:12 geuzaine Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -60,6 +60,30 @@ static int getElementTypeForMED(int msh, med_geometrie_element &med)
 }
 
 int GModel::readMED(const std::string &name)
+{
+  med_idt fid = MEDouvrir((char*)name.c_str(), MED_LECTURE);
+  if(fid < 0) {
+    Msg(GERROR, "Unable to open file '%s'", name.c_str());
+    return 0;
+  }
+
+  med_int numMeshes = MEDnMaa(fid);
+
+  if(MEDfermer(fid) < 0){
+    Msg(GERROR, "Unable to close file '%s'", (char*)name.c_str());
+    return 0;
+  }
+
+  int ret;
+  for(int i = 0; i < numMeshes; i++){
+    GModel *m = new GModel;
+    ret = m->readMED(name, i);
+    if(!ret) return 0;
+  }
+  return ret;
+}
+
+int GModel::readMED(const std::string &name, int meshIndex)
 {
   med_idt fid = MEDouvrir((char*)name.c_str(), MED_LECTURE);
   if(fid < 0) {
@@ -68,22 +92,20 @@ int GModel::readMED(const std::string &name)
   }
   
   int numMeshes = MEDnMaa(fid);
-  if(!numMeshes){
-    Msg(INFO, "No mesh found in MED file");
+  if(meshIndex >= numMeshes){
+    Msg(INFO, "Could not find mesh %d in MED file", meshIndex);
     return 0;
   }
-  if(numMeshes > 1)
-    Msg(WARNING, "Reading mesh 1 of %d (ignoring the others)", numMeshes);
 
   // read mesh info
   char meshName[MED_TAILLE_NOM + 1], meshDesc[MED_TAILLE_DESC + 1];
   med_int meshDim;
   med_maillage meshType;
-  if(MEDmaaInfo(fid, 1, meshName, &meshDim, &meshType, meshDesc) < 0){
+  if(MEDmaaInfo(fid, meshIndex + 1, meshName, &meshDim, &meshType, meshDesc) < 0){
     Msg(GERROR, "Unable to read mesh information");
     return 0;
   }
-
+  setName(meshName);
   if(meshType == MED_NON_STRUCTURE){
     Msg(INFO, "Reading %d-D unstructured mesh <<%s>>", meshDim, meshName);
   }
@@ -168,8 +190,10 @@ int GModel::readMED(const std::string &name)
   _associateEntityWithMeshVertices();
   for(unsigned int i = 0; i < verts.size(); i++){
     GEntity *ge = verts[i]->onWhat();
+    // store vertices (except for points, which ok already)
     if(ge && ge->dim() > 0) ge->mesh_vertices.push_back(verts[i]);
-    if(!ge) delete verts[i]; // delete unused vertices
+    // delete unused vertices
+    if(!ge) delete verts[i];
   }
 
   // read family info
@@ -198,7 +222,7 @@ int GModel::readMED(const std::string &name)
 	Msg(GERROR, "Could not read info for MED family %d", i + 1);
       }
       else{
-	GEntity *ge; // family tags are unique (for all dims)
+	GEntity *ge; // family tags are unique (for all dimensions)
 	if((ge = getRegionByTag(-familyNum))){}
 	else if((ge = getFaceByTag(-familyNum))){}
 	else if((ge = getEdgeByTag(-familyNum))){}
@@ -265,7 +289,7 @@ int GModel::writeMED(const std::string &name, bool saveAll, double scalingFactor
 
   char *meshName = (char*)getName().c_str();
   
-  // create 3D unstructured mesh
+  // Gmsh always writes 3D unstructured meshes
   if(MEDmaaCr(fid, meshName, 3, MED_NON_STRUCTURE, (char*)"gmsh") < 0){
     Msg(GERROR, "Could not create MED mesh");
     return 0;
@@ -433,6 +457,13 @@ int GModel::readMED(const std::string &name)
   return 0;
 }
 
+int GModel::readMED(const std::string &name, int meshIndex)
+{
+  Msg(GERROR, "Gmsh has to be compiled with MED support to read '%s'",
+      name.c_str());
+  return 0;
+}
+
 int GModel::writeMED(const std::string &name, bool saveAll, double scalingFactor)
 {
   Msg(GERROR, "Gmsh has to be compiled with MED support to write '%s'",
diff --git a/Parser/OpenFile.cpp b/Parser/OpenFile.cpp
index 3d49e2731818676cb5b88e5517f0c2255c3d02f6..96ebe181fad132a74e9641165a4ebdc84a62b799 100644
--- a/Parser/OpenFile.cpp
+++ b/Parser/OpenFile.cpp
@@ -1,4 +1,4 @@
-// $Id: OpenFile.cpp,v 1.181 2008-03-29 10:19:42 geuzaine Exp $
+// $Id: OpenFile.cpp,v 1.182 2008-03-30 17:45:12 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -353,7 +353,7 @@ int MergeFile(const char *name, int warn_if_missing)
     status = m->readMESH(name);
   }
   else if(!strcmp(ext, ".med") || !strcmp(ext, ".MED")){
-    status = m->readMED(name);
+    status = GModel::readMED(name);
     if(status > 1) status = PView::readMED(name);
   }
   else if(!strcmp(ext, ".bdf") || !strcmp(ext, ".BDF") ||
diff --git a/Post/PViewDataGModelIO.cpp b/Post/PViewDataGModelIO.cpp
index 0b2de783b2dab39248925557f5a3d47f238de818..d3b1e71d94e6b660d39e02d91a212d26dcc07007 100644
--- a/Post/PViewDataGModelIO.cpp
+++ b/Post/PViewDataGModelIO.cpp
@@ -1,4 +1,4 @@
-// $Id: PViewDataGModelIO.cpp,v 1.24 2008-03-30 13:21:04 geuzaine Exp $
+// $Id: PViewDataGModelIO.cpp,v 1.25 2008-03-30 17:45:12 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -212,7 +212,6 @@ bool PViewDataGModel::readMED(std::string fileName, int fileIndex)
 	    ElementNodeData);
   }
 
-  _steps.clear();
   for(int step = 0; step < numSteps; step++){
     for(unsigned int pair = 0; pair < pairs.size(); pair++){
       // get step info
@@ -229,11 +228,15 @@ bool PViewDataGModel::readMED(std::string fileName, int fileIndex)
       }
 
       // create step data
-      // FIXME: search for model (i.e., mesh) by name and use the right one!
       if(!pair){
 	int numCompMsh = (numComp == 1) ? 1 : (numComp < 3) ? 3 : 9;
+	GModel *m = GModel::findByName(meshName);
+	if(!m){
+	  Msg(GERROR, "Could not find mesh <<%s>>", meshName);
+	  return false;
+	}
 	while(step >= (int)_steps.size())
-	  _steps.push_back(new stepData<double>(GModel::current(), numCompMsh));
+	  _steps.push_back(new stepData<double>(m, numCompMsh));
 	_steps[step]->setFileName(fileName);
 	_steps[step]->setFileIndex(fileIndex);
 	_steps[step]->setTime(dt);
@@ -318,6 +321,9 @@ bool PViewDataGModel::readMED(std::string fileName, int fileIndex)
 	  num = tags[profile[i] - 1];
 	}
 	double *d = _steps[step]->getData(num, true, mult);
+	// FIXME: for ElementNodeData, we need to reorder the data before
+	// storing it
+	// FIXME: what should we do with Gauss point data?
 	for(int j = 0; j < numComp * mult; j++)
 	  d[j] = val[numComp * i + j];
 	double s = ComputeScalarRep(_steps[step]->getNumComponents(), d);