diff --git a/Geo/gmshLevelset.cpp b/Geo/gmshLevelset.cpp
index a7dbd96eaf55305cc7f463767dbd3a9cf477a445..5fdadd4934fea33ac744e122b3f4cc0e4cce0a31 100644
--- a/Geo/gmshLevelset.cpp
+++ b/Geo/gmshLevelset.cpp
@@ -1157,3 +1157,43 @@ gLevelsetConrod::gLevelsetConrod(const double *pt, const double *dir1,
 
 gLevelsetConrod::gLevelsetConrod(const gLevelsetConrod &lv)
   : gLevelsetImproved(lv){}
+
+
+double gLevelsetNACA00::operator() (const double x, const double y, const double z) const
+{
+
+  static const int maxIter = 100;
+  static const double tol = 1.e-8;
+
+  const double tolr = tol/_c;                 // Tolerance (scaled bu chord)
+  double distSq;                              // Square of dist. between point and boundary
+  bool in = false;                            // Whether the point is inside the airfoil
+
+  const double xt = x-_x0, yt = fabs(y-_y0);  // Point translated according to airfoil origin and symmetry
+
+  if (xt-_c > 1.16925*_t*yt)                  // Behind line normal to airfoil at trailing edge, closest
+    distSq = (xt-_c)*(xt-_c)+yt*yt;           // boundary point is trailing edge
+  else {                                      // Otherwise Newton-Raphson to find closest boundary point
+    const double fact = 5.*_t*_c;
+    double xb = std::max(xt,tolr);
+    for (int it=0; it<maxIter; it++) {
+      const double xbr = xb/_c, sxbr = sqrt(xbr), xbr32 = xbr*sxbr,
+                    xbr2 = xbr*xbr, xbr3 = xbr2*xbr, xbr4 = xbr2*xbr2;
+      double yb = fact*(0.2969*sxbr-0.1260*xbr-0.3516*xbr2+0.2843*xbr3-0.1015*xbr4);
+      const double dyb = fact*(0.14845/sxbr-0.406*xbr3+0.8529*xbr2-0.7032*xbr-0.126)/_c;
+      const double ddyb = fact*(-0.074225/xbr32-1.218*xbr2+1.7058*xbr-0.7032)/(_c*_c);
+      const double xx = xt-xb, yy = yt-yb;
+      in = (xt > 0) && (yy < 0);
+      distSq = xx*xx+yy*yy;
+      const double dDistSq = -2.*(xx+dyb*yy);
+      const double ddDistSq = 2.*(1.-ddyb*yy+dyb*dyb);
+      const double mIncr = dDistSq/ddDistSq;
+      if (fabs(mIncr) < tolr) break;
+      else xb -= mIncr;
+      if (xb < tolr) xb = tolr;
+    }
+  }
+
+  return in ? -sqrt(distSq) : sqrt(distSq);
+
+}
diff --git a/Geo/gmshLevelset.h b/Geo/gmshLevelset.h
index 476992e55f52e6b161a35426485e0c4a9d6faa7c..d4cb0f2d33640eaa2b533a3895db0f0b2877256e 100644
--- a/Geo/gmshLevelset.h
+++ b/Geo/gmshLevelset.h
@@ -600,4 +600,14 @@ public:
 /*   int type2() const {return DISK;} */
 /* }; */
 
+class gLevelsetNACA00: public gLevelsetPrimitive
+{
+  double _x0, _y0, _c, _t;
+public:
+  gLevelsetNACA00(double x0, double y0, double c, double t) : _x0(x0), _y0(y0), _c(c), _t(t) {}
+  ~gLevelsetNACA00() {}
+  double operator () (const double x, const double y, const double z) const;
+  int type() const { return UNKNOWN; }
+};
+
 #endif