diff --git a/examples/test_singlecellsolves.cpp b/examples/test_singlecellsolves.cpp
index 111f2c20e..e5a143d7b 100644
--- a/examples/test_singlecellsolves.cpp
+++ b/examples/test_singlecellsolves.cpp
@@ -211,6 +211,7 @@ main(int argc, char** argv)
         const double ff = s; // Simplified a lot...
         for (int conc = 0; conc < num_concs; ++conc) {
             const double c = poly_props.cMax()*double(conc)/double(num_concs - 1);
+            std::vector<double> polymer_inflow_c(num_cells, c);
             // std::cout << "(s, c) = (" << s << ", " << c << ")\n";
             transport_src[0] = src[0]*ff;
             // Resetting the state for next run.
@@ -223,8 +224,8 @@ main(int argc, char** argv)
             reorder_model.solve(&state.faceflux()[0],
                                 &porevol[0],
                                 &transport_src[0],
+                                &polymer_inflow_c[0],
                                 dt,
-                                c,
                                 state.saturation(),
                                 state.concentration(),
                                 state.maxconcentration());
diff --git a/opm/polymer/SimulatorCompressiblePolymer.cpp b/opm/polymer/SimulatorCompressiblePolymer.cpp
index 8e07c55a3..843155506 100644
--- a/opm/polymer/SimulatorCompressiblePolymer.cpp
+++ b/opm/polymer/SimulatorCompressiblePolymer.cpp
@@ -129,7 +129,7 @@ namespace Opm
         std::vector< std::vector<int> > columns_;
         // Misc. data
         std::vector<int> allcells_;
-        PolymerInflow poly_inflow_;
+        boost::scoped_ptr<PolymerInflowInterface> poly_inflow_;
     };
 
 
@@ -191,10 +191,7 @@ namespace Opm
           tsolver_(grid, props, poly_props, *rock_comp_props,
                    TransportModelCompressiblePolymer::Bracketing,
                    param.getDefault("nl_tolerance", 1e-9),
-                   param.getDefault("nl_maxiter", 30)),
-          poly_inflow_(param.getDefault("poly_start_days", 300.0)*Opm::unit::day,
-                       param.getDefault("poly_end_days", 800.0)*Opm::unit::day,
-                       param.getDefault("poly_amount", poly_props.cMax()))
+                   param.getDefault("nl_maxiter", 30))
     {
         // For output.
         output_ = param.getDefault("output", true);
@@ -230,6 +227,11 @@ namespace Opm
             extractColumn(grid_, columns_);
         }
 
+        // Polymer inflow control.
+        poly_inflow_.reset(new PolymerInflowBasic(param.getDefault("poly_start_days", 300.0)*Opm::unit::day,
+                                                  param.getDefault("poly_end_days", 800.0)*Opm::unit::day,
+                                                  param.getDefault("poly_amount", poly_props.cMax())));
+
         // Misc init.
         const int num_cells = grid.number_of_cells;
         allcells_.resize(num_cells);
@@ -244,7 +246,8 @@ namespace Opm
                                                             PolymerBlackoilState& state,
                                                             WellState& well_state)
     {
-        std::vector<double> transport_src;
+        std::vector<double> transport_src(grid_.number_of_cells);
+        std::vector<double> polymer_inflow_c(grid_.number_of_cells);
 
         // Initialisation.
         std::vector<double> initial_pressure;
@@ -320,12 +323,7 @@ namespace Opm
             // Find inflow rate.
             const double current_time = timer.currentTime();
             double stepsize = timer.currentStepLength();
-            const double inflowc0 = poly_inflow_(current_time + 1e-5*stepsize);
-            const double inflowc1 = poly_inflow_(current_time + (1.0 - 1e-5)*stepsize);
-            if (inflowc0 != inflowc1) {
-                std::cout << "**** Warning: polymer inflow rate changes during timestep. Using rate near start of step.";
-            }
-            const double inflow_c = inflowc0;
+            poly_inflow_->getInflowValues(current_time, current_time + stepsize, polymer_inflow_c);
 
             // Solve transport.
             transport_timer.start();
@@ -340,7 +338,7 @@ namespace Opm
             for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) {
                 tsolver_.solve(&state.faceflux()[0], initial_pressure,
                                state.pressure(), &initial_porevol[0], &porevol[0],
-                               &transport_src[0], stepsize, inflow_c,
+                               &transport_src[0], &polymer_inflow_c[0], stepsize,
                                state.saturation(), state.surfacevol(),
                                state.concentration(), state.maxconcentration());
                 double substep_injected[2] = { 0.0 };
@@ -350,7 +348,7 @@ namespace Opm
                 Opm::computeInjectedProduced(props_, poly_props_,
                                              state.pressure(), state.surfacevol(), state.saturation(),
                                              state.concentration(), state.maxconcentration(),
-                                             transport_src, stepsize, inflow_c,
+                                             transport_src, polymer_inflow_c, stepsize,
                                              substep_injected, substep_produced,
                                              substep_polyinj, substep_polyprod);
                 injected[0] += substep_injected[0];
diff --git a/opm/polymer/SimulatorPolymer.cpp b/opm/polymer/SimulatorPolymer.cpp
index 83e759d56..fa1f14e1a 100644
--- a/opm/polymer/SimulatorPolymer.cpp
+++ b/opm/polymer/SimulatorPolymer.cpp
@@ -126,7 +126,7 @@ namespace Opm
         std::vector< std::vector<int> > columns_;
         // Misc. data
         std::vector<int> allcells_;
-        PolymerInflow poly_inflow_;
+        boost::scoped_ptr<PolymerInflowInterface> poly_inflow_;
     };
 
 
@@ -187,10 +187,7 @@ namespace Opm
                    gravity, wells, src, bcs),
           tsolver_(grid, props, poly_props, TransportModelPolymer::Bracketing,
                    param.getDefault("nl_tolerance", 1e-9),
-                   param.getDefault("nl_maxiter", 30)),
-          poly_inflow_(param.getDefault("poly_start_days", 300.0)*Opm::unit::day,
-                       param.getDefault("poly_end_days", 800.0)*Opm::unit::day,
-                       param.getDefault("poly_amount", poly_props.cMax()))
+                   param.getDefault("nl_maxiter", 30))
     {
         // For output.
         output_ = param.getDefault("output", true);
@@ -226,6 +223,11 @@ namespace Opm
             extractColumn(grid_, columns_);
         }
 
+        // Polymer inflow control.
+        poly_inflow_.reset(new PolymerInflowBasic(param.getDefault("poly_start_days", 300.0)*Opm::unit::day,
+                                                  param.getDefault("poly_end_days", 800.0)*Opm::unit::day,
+                                                  param.getDefault("poly_amount", poly_props.cMax())));
+
         // Misc init.
         const int num_cells = grid.number_of_cells;
         allcells_.resize(num_cells);
@@ -240,7 +242,8 @@ namespace Opm
                                                 PolymerState& state,
                                                 WellState& well_state)
     {
-        std::vector<double> transport_src;
+        std::vector<double> transport_src(grid_.number_of_cells);
+        std::vector<double> polymer_inflow_c(grid_.number_of_cells);
 
         // Initialisation.
         std::vector<double> porevol;
@@ -316,12 +319,7 @@ namespace Opm
             // Find inflow rate.
             const double current_time = timer.currentTime();
             double stepsize = timer.currentStepLength();
-            const double inflowc0 = poly_inflow_(current_time + 1e-5*stepsize);
-            const double inflowc1 = poly_inflow_(current_time + (1.0 - 1e-5)*stepsize);
-            if (inflowc0 != inflowc1) {
-                std::cout << "**** Warning: polymer inflow rate changes during timestep. Using rate near start of step.";
-            }
-            const double inflow_c = inflowc0;
+            poly_inflow_->getInflowValues(current_time, current_time + stepsize, polymer_inflow_c);
 
             // Solve transport.
             transport_timer.start();
@@ -335,11 +333,11 @@ namespace Opm
             double substep_polyprod = 0.0;
             injected[0] = injected[1] = produced[0] = produced[1] = polyinj = polyprod = 0.0;
             for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) {
-                tsolver_.solve(&state.faceflux()[0], &porevol[0], &transport_src[0], stepsize, inflow_c,
+                tsolver_.solve(&state.faceflux()[0], &porevol[0], &transport_src[0], &polymer_inflow_c[0], stepsize,
                                state.saturation(), state.concentration(), state.maxconcentration());
                 Opm::computeInjectedProduced(props_, poly_props_,
                                              state.saturation(), state.concentration(), state.maxconcentration(),
-                                             transport_src, stepsize, inflow_c,
+                                             transport_src, polymer_inflow_c, stepsize,
                                              substep_injected, substep_produced, substep_polyinj, substep_polyprod);
                 injected[0] += substep_injected[0];
                 injected[1] += substep_injected[1];
diff --git a/opm/polymer/TransportModelCompressiblePolymer.cpp b/opm/polymer/TransportModelCompressiblePolymer.cpp
index 755e357e9..85f08d7fe 100644
--- a/opm/polymer/TransportModelCompressiblePolymer.cpp
+++ b/opm/polymer/TransportModelCompressiblePolymer.cpp
@@ -168,8 +168,8 @@ namespace Opm
           porevolume0_(0),
           porevolume_(0),
 	  source_(0),
+	  polymer_inflow_c_(0),
 	  dt_(0.0),
-	  inflow_c_(0.0),
 	  tol_(tol),
 	  maxit_(maxit),
 	  method_(method),
@@ -220,8 +220,8 @@ namespace Opm
                                                   const double* porevolume0,
                                                   const double* porevolume,
                                                   const double* source,
+                                                  const double* polymer_inflow_c,
                                                   const double dt,
-                                                  const double inflow_c,
                                                   std::vector<double>& saturation,
                                                   std::vector<double>& surfacevol,
                                                   std::vector<double>& concentration,
@@ -232,7 +232,7 @@ namespace Opm
         porevolume_ = porevolume;
 	source_ = source;
 	dt_ = dt;
-	inflow_c_ = inflow_c;
+	polymer_inflow_c_ = polymer_inflow_c;
         toWaterSat(saturation, saturation_);
 	concentration_ = &concentration[0];
 	cmax_ = &cmax[0];
@@ -363,7 +363,7 @@ namespace Opm
         bool src_is_inflow = src_flux < 0.0;
         B_cell0 = 1.0/tm.A0_[np*np*cell + 0];
         B_cell = 1.0/tm.A_[np*np*cell + 0];
-        // Not clear why we multiply by B_cell source terms.
+        // influx  =  src_is_inflow ? B_cell*src_flux : 0.0; // Use this after changing transport source.
         influx  =  src_is_inflow ? src_flux : 0.0;
         outflux = !src_is_inflow ? src_flux : 0.0;
         porevolume0 = tm.porevolume0_[cell];
@@ -376,7 +376,7 @@ namespace Opm
 	rhor = tm.polyprops_.rockDensity();
         tm.polyprops_.adsorption(c0, cmax0, ads0);
         double mc;
-        tm.computeMc(tm.inflow_c_, mc);
+        tm.computeMc(tm.polymer_inflow_c_[cell_index], mc);
 	influx_polymer = src_is_inflow ? src_flux*mc : 0.0;
 	for (int i = tm.grid_.cell_facepos[cell]; i < tm.grid_.cell_facepos[cell+1]; ++i) {
 	    int f = tm.grid_.cell_faces[i];
diff --git a/opm/polymer/TransportModelCompressiblePolymer.hpp b/opm/polymer/TransportModelCompressiblePolymer.hpp
index 063c2008d..469ea74cc 100644
--- a/opm/polymer/TransportModelCompressiblePolymer.hpp
+++ b/opm/polymer/TransportModelCompressiblePolymer.hpp
@@ -78,9 +78,14 @@ namespace Opm
 	/// \param[in] pressure            Array with pressure.
 	/// \param[in] porevolume0         Array with pore volume at start of timestep.
 	/// \param[in] porevolume          Array with pore volume.
-	/// \param[in] source              Transport source term.
+	/// \param[in] source              Transport source term, to be interpreted by sign:
+        ///                                 (+) Inflow, value is first phase flow (water)
+        ///                                     per second, in *surface* volumes (unlike the
+        ///                                     incompressible version).
+        ///                                 (-) Outflow, value is total flow of all phases
+        ///                                     per second, in reservoir volumes.
+	/// \param[in] polymer_inflow_c    Array of inflow polymer concentrations per cell.
 	/// \param[in] dt                  Time step.
-	/// \param[in] inflow_c            Inflow polymer.
 	/// \param[in, out] saturation     Phase saturations.
 	/// \param[in, out] surfacevol     Surface volumes.
 	/// \param[in, out] concentration  Polymer concentration.
@@ -91,8 +96,8 @@ namespace Opm
                    const double* porevolume0,
                    const double* porevolume,
 		   const double* source,
+		   const double* polymer_inflow_c,
 		   const double dt,
-		   const double inflow_c,
 		   std::vector<double>& saturation,
 		   std::vector<double>& surfacevol,
                    std::vector<double>& concentration,
@@ -134,8 +139,8 @@ namespace Opm
         const double* porevolume0_; // one volume per cell
         const double* porevolume_;  // one volume per cell
 	const double* source_;      // one source per cell
+	const double* polymer_inflow_c_;
 	double dt_;
-	double inflow_c_;
 	double tol_;
 	double maxit_;
 	SingleCellMethod method_;
diff --git a/opm/polymer/TransportModelPolymer.cpp b/opm/polymer/TransportModelPolymer.cpp
index 5227f5411..af79968d9 100644
--- a/opm/polymer/TransportModelPolymer.cpp
+++ b/opm/polymer/TransportModelPolymer.cpp
@@ -189,8 +189,8 @@ namespace Opm
 	  maxit_(maxit),
 	  darcyflux_(0),
 	  source_(0),
+	  polymer_inflow_c_(0),
 	  dt_(0.0),
-	  inflow_c_(0.0),
 	  concentration_(0),
 	  cmax_(0),
 	  fractionalflow_(grid.number_of_cells, -1.0),
@@ -232,8 +232,8 @@ namespace Opm
     void TransportModelPolymer::solve(const double* darcyflux,
                                       const double* porevolume,
 				      const double* source,
+				      const double* polymer_inflow_c,
 				      const double dt,
-				      const double inflow_c,
 				      std::vector<double>& saturation,
 				      std::vector<double>& concentration,
 				      std::vector<double>& cmax)
@@ -241,8 +241,8 @@ namespace Opm
 	darcyflux_ = darcyflux;
         porevolume_ = porevolume;
 	source_ = source;
+	polymer_inflow_c_ = polymer_inflow_c;
 	dt_ = dt;
-	inflow_c_ = inflow_c;
         toWaterSat(saturation, saturation_);
 	concentration_ = &concentration[0];
 	cmax_ = &cmax[0];
@@ -350,7 +350,7 @@ namespace Opm
 	bool src_is_inflow = dflux < 0.0;
 	influx  =  src_is_inflow ? dflux : 0.0;
         double mc;
-        tm.computeMc(tm.inflow_c_, mc);
+        tm.computeMc(tm.polymer_inflow_c_[cell_index], mc);
 	influx_polymer = src_is_inflow ? dflux*mc : 0.0;
 	outflux = !src_is_inflow ? dflux : 0.0;
 	comp_term = tm.source_[cell];   // Note: this assumes that all source flux is water.
diff --git a/opm/polymer/TransportModelPolymer.hpp b/opm/polymer/TransportModelPolymer.hpp
index da3ff2a63..c9d3fd5c7 100644
--- a/opm/polymer/TransportModelPolymer.hpp
+++ b/opm/polymer/TransportModelPolymer.hpp
@@ -67,17 +67,21 @@ namespace Opm
 	/// Using implicit Euler scheme, reordered.
 	/// \param[in] darcyflux           Array of signed face fluxes.
 	/// \param[in] porevolume          Array of pore volumes.
-	/// \param[in] source              Transport source term.
+	/// \param[in] source              Transport source term, to be interpreted by sign:
+        ///                                 (+) Inflow, value is first phase flow (water)
+        ///                                     per second, in reservoir volumes.
+        ///                                 (-) Outflow, value is total flow of all phases
+        ///                                     per second, in reservoir volumes.
+	/// \param[in] polymer_inflow_c    Array of inflow polymer concentrations per cell.
 	/// \param[in] dt                  Time step.
-	/// \param[in] inflow_c            Time step.
 	/// \param[in, out] saturation     Phase saturations.
 	/// \param[in, out] concentration  Polymer concentration.
 	/// \param[in, out] cmax           Highest concentration that has occured in a given cell.
 	void solve(const double* darcyflux,
                    const double* porevolume,
 		   const double* source,
+                   const double* polymer_inflow_c,
 		   const double dt,
-		   const double inflow_c,
 		   std::vector<double>& saturation,
                    std::vector<double>& concentration,
                    std::vector<double>& cmax);
@@ -149,8 +153,8 @@ namespace Opm
 
 	const double* darcyflux_;   // one flux per grid face
 	const double* source_;      // one source per cell
+	const double* polymer_inflow_c_;
 	double dt_;
-	double inflow_c_;
         std::vector<double> saturation_; // one per cell, only water saturation!
 	double* concentration_;
 	double* cmax_;
diff --git a/opm/polymer/polymerUtilities.cpp b/opm/polymer/polymerUtilities.cpp
index 910556a7c..3b103b35f 100644
--- a/opm/polymer/polymerUtilities.cpp
+++ b/opm/polymer/polymerUtilities.cpp
@@ -103,8 +103,8 @@ namespace Opm
     /// @param[in]  c         polymer concentration
     /// @param[in]  cmax      polymer maximum concentration
     /// @param[in]  src       if < 0: total outflow, if > 0: first phase inflow.
+    /// @param[in]  inj_c     injected concentration by cell
     /// @param[in]  dt        timestep used
-    /// @param[in]  inj_c     injected concentration
     /// @param[out] injected  must point to a valid array with P elements,
     ///                       where P = s.size()/src.size().
     /// @param[out] produced  must also point to a valid array with P elements.
@@ -116,8 +116,8 @@ namespace Opm
 				 const std::vector<double>& c,
 				 const std::vector<double>& cmax,
 				 const std::vector<double>& src,
+                                 const std::vector<double>& inj_c,
 				 const double dt,
-                                 const double inj_c,
 				 double* injected,
 				 double* produced,
                                  double& polyinj,
@@ -139,7 +139,7 @@ namespace Opm
         for (int cell = 0; cell < num_cells; ++cell) {
             if (src[cell] > 0.0) {
                 injected[0] += src[cell]*dt;
-                polyinj += src[cell]*dt*inj_c;
+                polyinj += src[cell]*dt*inj_c[cell];
             } else if (src[cell] < 0.0) {
                 const double flux = -src[cell]*dt;
                 const double* sat = &s[np*cell];
@@ -170,15 +170,13 @@ namespace Opm
     /// @param[in]  c         polymer concentration
     /// @param[in]  cmax      polymer maximum concentration
     /// @param[in]  src       if < 0: total outflow, if > 0: first phase inflow.
+    /// @param[in]  inj_c     injected concentration by cell
     /// @param[in]  dt        timestep used
-    /// @param[in]  inj_c     injected concentration
-    ///
     /// @param[out] injected  must point to a valid array with P elements,
     ///                       where P = s.size()/src.size().
     /// @param[out] produced  must also point to a valid array with P elements.
     /// @param[out] polyinj   injected mass of polymer
     /// @param[out] polyprod  produced mass of polymer
-
     void computeInjectedProduced(const BlackoilPropertiesInterface& props,
                                  const Opm::PolymerProperties& polyprops,
                                  const std::vector<double>& press,
@@ -187,8 +185,8 @@ namespace Opm
 				 const std::vector<double>& c,
 				 const std::vector<double>& cmax,
                                  const std::vector<double>& src,
+                                 const std::vector<double>& inj_c,
                                  const double dt,
-                                 const double inj_c,
                                  double* injected,
                                  double* produced,
                                  double& polyinj,
@@ -210,7 +208,7 @@ namespace Opm
         for (int cell = 0; cell < num_cells; ++cell) {
             if (src[cell] > 0.0) {
                 injected[0] += src[cell]*dt;
-                polyinj += src[cell]*dt*inj_c;
+                polyinj += src[cell]*dt*inj_c[cell];
             } else if (src[cell] < 0.0) {
                 const double flux = -src[cell]*dt;
                 const double* sat = &s[np*cell];
diff --git a/opm/polymer/polymerUtilities.hpp b/opm/polymer/polymerUtilities.hpp
index 66729af4c..2bde8c2a7 100644
--- a/opm/polymer/polymerUtilities.hpp
+++ b/opm/polymer/polymerUtilities.hpp
@@ -77,8 +77,8 @@ namespace Opm
     /// @param[in]  s         saturation values (for all P phases)
     /// @param[in]  c         polymer concentration
     /// @param[in]  src       if < 0: total outflow, if > 0: first phase inflow.
+    /// @param[in]  inj_c     injected concentration by cell
     /// @param[in]  dt        timestep used
-    /// @param[in]  inj_c     injected concentration
     /// @param[out] injected  must point to a valid array with P elements,
     ///                       where P = s.size()/src.size().
     /// @param[out] produced  must also point to a valid array with P elements.
@@ -90,8 +90,8 @@ namespace Opm
 				 const std::vector<double>& c,
                                  const std::vector<double>& cmax,
 				 const std::vector<double>& src,
+				 const std::vector<double>& inj_c,
 				 const double dt,
-                                 const double inj_c,
 				 double* injected,
 				 double* produced,
                                  double& polyinj,
@@ -111,8 +111,8 @@ namespace Opm
     /// @param[in]  c         polymer concentration
     /// @param[in]  cmax      polymer maximum concentration
     /// @param[in]  src       if < 0: total outflow, if > 0: first phase inflow.
+    /// @param[in]  inj_c     injected concentration by cell
     /// @param[in]  dt        timestep used
-    /// @param[in]  inj_c     injected concentration
     ///
     /// @param[out] injected  must point to a valid array with P elements,
     ///                       where P = s.size()/src.size().
@@ -127,9 +127,9 @@ namespace Opm
                                  const std::vector<double>& s,
 				 const std::vector<double>& c,
 				 const std::vector<double>& cmax,
-                                 const std::vector<double>& src,
-                                 const double dt,
-                                 const double inj_c,
+				 const std::vector<double>& src,
+				 const std::vector<double>& inj_c,
+				 const double dt,
                                  double* injected,
                                  double* produced,
                                  double& polyinj,
@@ -172,28 +172,44 @@ namespace Opm
                                   const RockCompressibility* rock_comp);
 
 
-    /// @brief Functor giving the injected amount of polymer as a function of time.
-    class PolymerInflow
+    class PolymerInflowInterface
+    {
+    public:
+        virtual ~PolymerInflowInterface() {}
+        virtual void getInflowValues(const double step_start,
+                                     const double step_end,
+                                     std::vector<double>& poly_inflow_c) = 0;
+    };
+
+
+
+    /// @brief Functor giving the injected amount of polymer per cell as a function of time.
+    class PolymerInflowBasic : public PolymerInflowInterface
     {
     public:
         /// Constructor.
         /// @param[in]  starttime  Start time of injection in seconds.
         /// @param[in]  endtime    End time of injection in seconds.
         /// @param[in]  amount     Amount to be injected per second.
-        PolymerInflow(const double starttime,
-                      const double endtime,
-                      const double amount)
+        PolymerInflowBasic(const double starttime,
+                           const double endtime,
+                           const double amount)
             : stime_(starttime), etime_(endtime), amount_(amount)
         {
         }
-        /// Get the current injection rate.
-        /// @param[in]  time   Current time in seconds.
-        double operator()(double time)
+
+        virtual void getInflowValues(const double step_start,
+                                     const double step_end,
+                                     std::vector<double>& poly_inflow_c)
         {
-            if (time >= stime_ && time < etime_) {
-                return amount_;
+            const double eps = 1e-5*(step_end - step_start);
+            if (step_start + eps >= stime_ && step_end - eps <= etime_) {
+                std::fill(poly_inflow_c.begin(), poly_inflow_c.end(), amount_);
+            } else if (step_start + eps <= etime_ && step_end - eps >= stime_) {
+                MESSAGE("Warning: polymer injection set to change inside timestep. Using value at start of step.");
+                std::fill(poly_inflow_c.begin(), poly_inflow_c.end(), amount_);
             } else {
-                return 0.0;
+                std::fill(poly_inflow_c.begin(), poly_inflow_c.end(), 0.0);
             }
         }
     private: