diff --git a/examples/Makefile.am b/examples/Makefile.am
index 5515aea..b2ed501 100644
--- a/examples/Makefile.am
+++ b/examples/Makefile.am
@@ -6,8 +6,9 @@ noinst_PROGRAMS = upscale_perm \
                   upscale_relperm \
                   upscale_avg  \
                   upscale_cond \
-		  cpchop \
-                  exp_variogram
+                  cpchop \
+                  exp_variogram \
+                  cp_regularize
 
 upscale_perm_SOURCES = upscale_perm.cpp
 
@@ -21,6 +22,8 @@ cpchop_SOURCES = cpchop.cpp
 
 exp_variogram_SOURCES = exp_variogram.cpp
 
+cp_regularize_SOURCES = cp_regularize.cpp
+
 AM_CPPFLAGS += $(DUNEMPICPPFLAGS) $(BOOST_CPPFLAGS) $(SUPERLU_CPPFLAGS)
 AM_LDFLAGS  += $(DUNEMPILDFLAGS) $(BOOST_LDFLAGS) $(SUPERLU_LDFLAGS)
 LDADD = $(DUNEMPILIBS) $(BOOST_UNIT_TEST_FRAMEWORK_LIB) \
diff --git a/examples/cp_regularize.cpp b/examples/cp_regularize.cpp
index 7fc5113..9a4fb39 100644
--- a/examples/cp_regularize.cpp
+++ b/examples/cp_regularize.cpp
@@ -74,7 +74,105 @@ int main(int argc, char** argv)
     std::vector<double> permy;
     std::vector<double> permz;
 
-    // Run through the new regular grid to find its properties
 
+    // Construct mapping from coarse i and j indices to fine
+    std::vector<int> iidx_f, jidx_f;
+    int finesprcoarse_i = floor(dims[0] / ires);
+    int remainder_i = dims[0] - ires*finesprcoarse;
+    for (int iidx_c=0; iidx_c < remainder_i; ++iidx_c) {
+        iidx_f.push_back(finesprcoarse_i + 1);  
+    }
+    for (int iidx_c=remainder_i; iidx_c < ires; ++iidx_c) {
+        iidx_f.push_back(finesprcoarse_i);
+    }
+    int finesprcoarse_j = floor(dims[1] / jres);
+    int remainder_j = dims[1] - jres*finesprcoarse_j;
+    for (int jidx_c=0; iidx_c < remainder; ++iidx_c) {
+        jidx_f.push_back(finesprcoarse_j + 1);
+    }
+    for (int jidx_c=remainder; jidx_c < jres; ++jidx_c) {
+        jidx_f.push_back(finesprcoarse_j);
+    }
+
+
+
+    // Run through the new regular grid to find its properties
+    for (zidx_c=0; zidx_c < zres; ++zidx_c) {
+        for (jidx_c=0; jidx_c < jres; ++jidx_c) {
+            for (iidx_c=0; iidx_c < ires; ++iidx_c) {
+                cp.chop(iidx_f(iidx_c), iidx_f(iidx_c+1),
+                    jidx_f(jidx_c), jidx_f(jidx_c+1),
+                    zcorn_c(zidc_c), zcorn_c(zidx_c+1),
+                    false);
+                Dune::EclipseGridParser subparser = ch.subparser();
+                Dune::SinglePhaseUpscaler upscaler;
+                upscaler.init(subparser, Dune::SinglePhaseUpscaler::Fixed, 0.0, z_tolerance,
+                          residual_tolerance, linsolver_verbosity, linsolver_type, false);
+            
+                Dune::SinglePhaseUpscaler::permtensor_t upscaled_K = upscaler.upscaleSinglePhase();
+                upscaled_K *= (1.0/(Dune::prefix::milli*Dune::unit::darcy));
+                poro.push_back(upscaler.upscalePorosity());
+                permx.push_back(upscaled_K(0,0));
+                permy.push_back(upscaled_K(1,1));
+                permz.push_back(upscaled_K(2,2)); 
+            }
+        }
+    }
     // Write regularized grid to outputfile
+    std::ofstream out(resultgrid.c_str());
+    if (!out) {
+        std::cerr << "Could not open file " << filename << "\n";
+        throw std::runtime_error("Could not open output file.");
+    }
+    out << "SPECGRID\n" << ires << ' ' << jres << ' ' << zres 
+                << " 1 F\n/\n\n";
+
+    out << "COORD\n";
+    for (int j = 0; j < jres; ++j) {
+        for (int i = 0; i < ires; ++i) {
+            // Calculate these on the fly instead!
+            out << x(i) << " " << y(j) << " " << zmin << " "
+                << x(i) << " " << y(j) << " " << zmax << "\n";
+        }
+    }
+
+    /*
+     Write ZCORN, that is the Z-coordinates along the pillars, specifying
+     the eight corners of each cell. Each corner is specified for each
+     cell, even though it is the same corner that is used in other
+     cells. 
+
+     We loop over corners in each grid cell, directions: z, y, x (x innermost).
+     The code here *IS* redundant, but the grid is also very redundant
+     for a grid that is really uniform..
+   */ 
+  out << "ZCORN\n";
+ //  for zidx=1:numel(z)-1,
+ //    for yidx=1:numel(y)-1,
+  //       for xidx=1:numel(x)-1,
+   //         fprintf(outfile, ' %g %g', ...
+    //                z(zidx), z(zidx));
+     //    end
+      //   fprintf(outfile, '\n');
+       //  for xidx=1:numel(x)-1,
+        //    fprintf(outfile, ' %g %g', ...
+         //           z(zidx), z(zidx));
+   //      end
+    //  end
+     // for yidx=1:numel(y)-1,
+      //   for xidx=1:numel(x)-1,
+       //     fprintf(outfile, ' %g %g', ...
+        //            z(zidx+1), z(zidx+1));
+      //   end
+       //  fprintf(outfile, '\n');
+        // for xidx=1:numel(x)-1,
+         //   fprintf(outfile, ' %g %g', ...
+          //          z(zidx+1), z(zidx+1));
+     //    end
+     // end
+
+  // end
+   //fprintf(outfile, '\n/\n\n');
+    out.close();
 }
+