diff --git a/Common/Context.h b/Common/Context.h
index b24360a542fd5d21a0a3bd19f35df45a76443ec5..cf87afb2421f982da44e67ada65ddc68f613c31d 100644
--- a/Common/Context.h
+++ b/Common/Context.h
@@ -55,6 +55,7 @@ struct contextMeshOptions {
   std::map<int,int> curvature_control_per_face;
   int bunin;
   int ignorePartBound;
+  int NewtonConvergenceTestXYZ;
 };
 
 struct contextGeometryOptions {
diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index 721b818f1685d16791ba140247b8d485c9a58cb9..5164ddbae842f3f0f437d0ad2bb46925fae8847c 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -963,6 +963,8 @@ StringXNumber MeshOptions_Number[] = {
 
   { F|O, "FlexibleTransfinite" , opt_mesh_flexible_transfinite , 0 ,
     "Allow transfinite contraints to be modified for Blossom or by global mesh size factor" },
+  { F|O, "NewtonConvergenceTestXYZ" , opt_mesh_newton_convergence_test_xyz , 0. ,
+    "Force inverse surface mapping algorithm (Newton-Raphson) to converge in real coordinates (experimental)" },
   { F|O, "Format" , opt_mesh_file_format , FORMAT_AUTO ,
     "Mesh output format (1=msh, 2=unv, 10=automatic, 19=vrml, 27=stl, 30=mesh, 31=bdf, "
     "32=cgns, 33=med, 40=ply2)" },
diff --git a/Common/Options.cpp b/Common/Options.cpp
index 67773b4e9f5b36babd45d793bf0b1a528c2c9480..efbcfbb9c5b3557c841fbd45d074d94147daecf3 100644
--- a/Common/Options.cpp
+++ b/Common/Options.cpp
@@ -5250,6 +5250,13 @@ double opt_mesh_file_format(OPT_ARGS_NUM)
   return CTX::instance()->mesh.fileFormat;
 }
 
+double opt_mesh_newton_convergence_test_xyz(OPT_ARGS_NUM)
+{
+  if(action & GMSH_SET)
+    CTX::instance()->mesh.NewtonConvergenceTestXYZ = (int)val;
+  return CTX::instance()->mesh.NewtonConvergenceTestXYZ;
+}
+
 double opt_mesh_msh_file_version(OPT_ARGS_NUM)
 {
   if(action & GMSH_SET)
diff --git a/Common/Options.h b/Common/Options.h
index c6a5c5ce2bfaa115eb1b23a8314d64fcbdcbbcba..fccd0d8562a2070b4f621062343103bfefd5433c 100644
--- a/Common/Options.h
+++ b/Common/Options.h
@@ -415,6 +415,7 @@ double opt_mesh_light(OPT_ARGS_NUM);
 double opt_mesh_light_lines(OPT_ARGS_NUM);
 double opt_mesh_light_two_side(OPT_ARGS_NUM);
 double opt_mesh_file_format(OPT_ARGS_NUM);
+double opt_mesh_newton_convergence_test_xyz(OPT_ARGS_NUM);
 double opt_mesh_msh_file_version(OPT_ARGS_NUM);
 double opt_mesh_msh_file_partitioned(OPT_ARGS_NUM);
 double opt_mesh_partition_hex_weight(OPT_ARGS_NUM);
diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 841977ea5614687fba3a240b2175c1e9f55a1a38..8efeb53b2fc842db4e57c80e29c749d8cd7e98e6 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -877,11 +877,22 @@ void GFace::XYZtoUV(double X, double Y, double Z, double &U, double &V,
       if(iter < MaxIter && err <= tol &&
          Unew <= umax && Vnew <= vmax &&
          Unew >= umin && Vnew >= vmin){
-        if (onSurface && err2 > 1.e-4 * CTX::instance()->lc)
+
+        if (onSurface && err2 > 1.e-4 * CTX::instance()->lc &&
+            !CTX::instance()->mesh.NewtonConvergenceTestXYZ){
           Msg::Warning("Converged for i=%d j=%d (err=%g iter=%d) BUT "
                        "xyz error = %g in point (%e,%e,%e) on surface %d",
                        i, j, err, iter, err2, X, Y, Z, tag());
-        return;
+        }
+
+        if(onSurface && err2 > 1.e-4 * CTX::instance()->lc &&
+           CTX::instance()->mesh.NewtonConvergenceTestXYZ){
+          // not converged in XYZ coordinates
+        }
+        else{
+          return;
+        }
+
       }
     }
   }