diff --git a/pystare/PySTARE.cpp b/pystare/PySTARE.cpp
index 4d280bc..1d7440b 100644
--- a/pystare/PySTARE.cpp
+++ b/pystare/PySTARE.cpp
@@ -12,9 +12,12 @@
 
 #include <algorithm>
 
+
+#include <omp.h>
 #undef FMTX
 #define FMTX(x) " 0x" << setfill('0') << setw(16) << hex << x << dec
 
+
 // Info
 const char* stare_version() { return STARE_version(); }
 
@@ -25,34 +28,36 @@ void _from_latlon(double* lat, int len_lat, double * lon, int len_lon, int64_t*
     }
 }
 
-void _from_latlon2D(double* lat, int lalen1, int lalen2, 
-		    double* lon, int lolen1, int lolen2, 
-		    int64_t* indices, int len1, int len2, int level, bool adapt_resolution,
-		    bool fill_value_enabled,
-		    double fill_value_in,
-		    int fill_value_out // would prefer int64_t, but it doesn't work.
-		    ) {
+
+void _from_latlon2D(double* lat, int lalen1, int lalen2,
+		            double* lon, int lolen1, int lolen2,
+		            int64_t* indices, int len1, int len2, int level, bool adapt_resolution,
+		            bool fill_value_enabled,
+		            double fill_value_in,
+		            int fill_value_out // would prefer int64_t, but it doesn't work.
+        		    ) {
     int n;   
     static EmbeddedLevelNameEncoding lj;       // Use this to get the mask
     int lvl;
 
     if (!fill_value_enabled) {
-      for (int i=0; i<lalen1; i++) {        
-        for (int j=0; j<lalen2; j++) {   
-	  n = i * lalen2 + j;
-	  indices[n] = stare.ValueFromLatLonDegrees(lat[n], lon[n], level);   
-	  if (adapt_resolution) {
-	    if (j==0) {        
-	      indices[n+1] = stare.ValueFromLatLonDegrees(lat[n+1], lon[n+1], level);
-	      lvl = stare.cmpSpatialResolutionEstimateI(indices[n], indices[n+1]);
-	      indices[n] = (indices[n]& ~lj.levelMaskSciDB) | lvl;
-	    } else {                    
-	      lvl = stare.cmpSpatialResolutionEstimateI(indices[n-1], indices[n]);
-	      indices[n] = (indices[n]& ~lj.levelMaskSciDB) | lvl;
-	    }
-	  }
+        for (int i=0; i<lalen1; i++) {
+        #pragma omp parallel for
+            for (int j=0; j<lalen2; j++) {
+	            n = i * lalen2 + j;
+	            indices[n] = stare.ValueFromLatLonDegrees(lat[n], lon[n], level);
+	            if (adapt_resolution) {
+	                if (j==0) {
+	                    indices[n+1] = stare.ValueFromLatLonDegrees(lat[n+1], lon[n+1], level);
+	                    lvl = stare.cmpSpatialResolutionEstimateI(indices[n], indices[n+1]);
+	                    indices[n] = (indices[n]& ~lj.levelMaskSciDB) | lvl;
+	                } else {
+	                    lvl = stare.cmpSpatialResolutionEstimateI(indices[n-1], indices[n]);
+	                    indices[n] = (indices[n]& ~lj.levelMaskSciDB) | lvl;
+	                }
+	            }
+            }
         }
-      }
     } else {
       for (int i=0; i<lalen1; i++) {
 	bool skipped = true;