diff --git a/examples/polymer_reorder.cpp b/examples/polymer_reorder.cpp
index 7d06ee801..65e26d2c9 100644
--- a/examples/polymer_reorder.cpp
+++ b/examples/polymer_reorder.cpp
@@ -63,14 +63,16 @@
 
 
 
-class ReservoirState {
+class ReservoirState
+{
 public:
     ReservoirState(const UnstructuredGrid* g, const int num_phases = 2)
         : press_ (g->number_of_cells, 0.0),
           fpress_(g->number_of_faces, 0.0),
           flux_  (g->number_of_faces, 0.0),
           sat_   (num_phases * g->number_of_cells, 0.0),
-	  concentration_(g->number_of_cells, 0.0)
+	  concentration_(g->number_of_cells, 0.0),
+	  cmax_(g->number_of_cells, 0.0)
     {
 	for (int cell = 0; cell < g->number_of_cells; ++cell) {
 	    sat_[num_phases*cell + num_phases - 1] = 1.0;
@@ -79,24 +81,27 @@ public:
 
     int numPhases() const { return sat_.size()/press_.size(); }
 
-    ::std::vector<double>& pressure    () { return press_ ; }
-    ::std::vector<double>& facepressure() { return fpress_; }
-    ::std::vector<double>& faceflux    () { return flux_  ; }
-    ::std::vector<double>& saturation  () { return sat_   ; }
-    ::std::vector<double>& concentration() { return concentration_; }
+    std::vector<double>& pressure    () { return press_ ; }
+    std::vector<double>& facepressure() { return fpress_; }
+    std::vector<double>& faceflux    () { return flux_  ; }
+    std::vector<double>& saturation  () { return sat_   ; }
+    std::vector<double>& concentration() { return concentration_; }
+    std::vector<double>& cmax() { return cmax_; }
 
-    const ::std::vector<double>& pressure    () const { return press_ ; }
-    const ::std::vector<double>& facepressure() const { return fpress_; }
-    const ::std::vector<double>& faceflux    () const { return flux_  ; }
-    const ::std::vector<double>& saturation  () const { return sat_   ; }
-    const ::std::vector<double>& concentration() const { return concentration_; }
+    const std::vector<double>& pressure    () const { return press_ ; }
+    const std::vector<double>& facepressure() const { return fpress_; }
+    const std::vector<double>& faceflux    () const { return flux_  ; }
+    const std::vector<double>& saturation  () const { return sat_   ; }
+    const std::vector<double>& concentration() const { return concentration_; }
+    const std::vector<double>& cmax() const { return cmax_; }
 
 private:
-    ::std::vector<double> press_ ;
-    ::std::vector<double> fpress_;
-    ::std::vector<double> flux_  ;
-    ::std::vector<double> sat_   ;
-    ::std::vector<double> concentration_   ;
+    std::vector<double> press_ ;
+    std::vector<double> fpress_;
+    std::vector<double> flux_  ;
+    std::vector<double> sat_   ;
+    std::vector<double> concentration_;
+    std::vector<double> cmax_;
 };
 
 
@@ -169,12 +174,18 @@ main(int argc, char** argv)
 	props.reset(new Opm::IncompPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells));
 	polydata.c_max_limit = param.getDefault("c_max_limit", 1.0);
 	polydata.omega = param.getDefault("omega", 1.0);
-	polydata.c_vals.resize(2);
-	polydata.c_vals[0] = 0.0;
-	polydata.c_vals[0] = polydata.c_max_limit;
+	polydata.rhor = param.getDefault("rock_density", 1000.0);
+	polydata.dps = param.getDefault("dead_pore_space", 0.15);
+	polydata.c_vals_visc.resize(2);
+	polydata.c_vals_visc[0] = 0.0;
+	polydata.c_vals_visc[0] = polydata.c_max_limit;
 	polydata.visc_mult_vals.resize(2);
 	polydata.visc_mult_vals[0] = 1.0;
 	polydata.visc_mult_vals[1] = param.getDefault("c_max_viscmult", 30.0);
+	polydata.c_vals_ads = polydata.c_vals_visc;
+	polydata.ads_vals.resize(2);
+	polydata.ads_vals[0] = 1.0;
+	polydata.ads_vals[1] = param.getDefault("c_max_ads", 30.0);
     }
 
     // Extra rock init.
@@ -252,6 +263,7 @@ main(int argc, char** argv)
 	// source term following the same convention.
 	transport_timer.start();
 	polymertransport(&porevol[0],
+			 props->porosity(),
 			 &reorder_src[0],
 			 stepsize,
 			 const_cast<UnstructuredGrid*>(grid->c_grid()),
@@ -259,7 +271,8 @@ main(int argc, char** argv)
 			 &polydata,
 			 &state.faceflux()[0],
 			 &reorder_sat[0],
-			 &state.concentration()[0]);
+			 &state.concentration()[0],
+			 &state.cmax()[0]);
 	Opm::toBothSat(reorder_sat, state.saturation());
 	transport_timer.stop();
 	double tt = transport_timer.secsSinceStart();
diff --git a/opm/polymer/polymermodel.cpp b/opm/polymer/polymermodel.cpp
index 436154a7b..5a395113b 100644
--- a/opm/polymer/polymermodel.cpp
+++ b/opm/polymer/polymermodel.cpp
@@ -19,14 +19,9 @@ struct ParametersSRes
 {
     double c;
     double s0;
-    // double c0;
     double dtpv;    /* dt/pv(i)                  */
     double influx;  /* sum_j min(v_ij, 0)*f(s_j) */
-    // double influx_polymer;
     double outflux; /* sum_j max(v_ij, 0)        */
-    // double dps;
-    // double rhor;
-    // double phi;
     int cell;
     const Opm::IncompPropertiesInterface* props;
     const PolymerData* polydata;
@@ -34,10 +29,19 @@ struct ParametersSRes
 
 struct ParametersCRes
 {
+    double s0;
+    double c0;
+    double cmax0;
+    double dtpv;    /* dt/pv(i)                  */
+    double influx;  /* sum_j min(v_ij, 0)*f(s_j) */
+    double influx_polymer;
+    double outflux; /* sum_j max(v_ij, 0)        */
+    double porosity;
     PolymerSolverData* psdata;
     int cell;
     NonlinearSolverCtrl* ctrl;
     double s;
+    const PolymerData* polydata;
 };
 
 
@@ -51,6 +55,9 @@ static double fluxfun_props(double s,
 			    int cell,
 			    const Opm::IncompPropertiesInterface* props, 
 			    const PolymerData* polydata);
+static double compute_mc(double c,
+			 const Opm::IncompPropertiesInterface* props, 
+			 const PolymerData* polydata);
 
 void
 destroy_solverdata(struct PolymerSolverData *d)
@@ -58,6 +65,7 @@ destroy_solverdata(struct PolymerSolverData *d)
     if (d!=NULL)
     {
         free(d->fractionalflow);
+        free(d->mc);
     }
     free(d);
 }
@@ -68,10 +76,12 @@ init_solverdata(struct UnstructuredGrid *grid,
 		const PolymerData* polydata,
 		const double *darcyflux,
                 const double *porevolume,
+                const double *porosity,
 		const double *source,
                 const double dt,
 		double *saturation,
-		double *concentration)
+		double *concentration,
+		double *cmax)
 {
     int i;
     struct PolymerSolverData *d = (struct PolymerSolverData*) malloc(sizeof *d);
@@ -83,14 +93,18 @@ init_solverdata(struct UnstructuredGrid *grid,
 	d->polydata   = polydata;
         d->darcyflux  = darcyflux;
         d->porevolume = porevolume;
+        d->porosity   = porosity;
         d->source     = source;
         d->dt         = dt;
 
         d->saturation     = saturation;
 	d->concentration  = concentration;
+	d->cmax           = cmax;
         d->fractionalflow = (double*) malloc(grid->number_of_cells *
 					     sizeof *d->fractionalflow);
-        if (d->fractionalflow == NULL)
+        d->mc             = (double*) malloc(grid->number_of_cells *
+					     sizeof *d->mc);
+        if (d->fractionalflow == NULL || d->mc == NULL)
         {
             destroy_solverdata(d);
             d = NULL;
@@ -110,8 +124,10 @@ void polymer_solvecell(void *data, struct NonlinearSolverCtrl *ctrl, int cell)
     struct ParametersCRes   prm = get_parameters_c(d, cell, ctrl);
 
     d->concentration[cell] = find_zero(residual_c, &prm, ctrl);
+    d->cmax[cell] = std::max(d->cmax[cell], d->concentration[cell]);
     d->saturation[cell] = prm.s;
     d->fractionalflow[cell] = fluxfun_props(d->saturation[cell], d->concentration[cell], cell, d->props, d->polydata);
+    d->mc[cell] = compute_mc(d->concentration[cell], d->props, d->polydata);
 }
 
 
@@ -134,13 +150,25 @@ residual_s(double s, void *data)
 static double
 residual_c(double c, void *data)
 {
-    struct ParametersCRes *prm_c = (struct ParametersCRes*) data;
-    int cell = prm_c->cell;
-    struct ParametersSRes prm_s = get_parameters_s(prm_c->psdata, cell);
+    struct ParametersCRes *p = (struct ParametersCRes*) data;
+    int cell = p->cell;
+    struct ParametersSRes prm_s = get_parameters_s(p->psdata, cell);
     prm_s.c = c;
-    double s = find_zero(residual_s, &prm_s, prm_c->ctrl);
-    prm_c->s = s;
-    return  c - 0.0;
+    double s = find_zero(residual_s, &prm_s, p->ctrl);
+    p->s = s;
+    double ff = fluxfun_props(s, c, p->cell, prm_s.props, p->polydata);
+    double mc = compute_mc(c, prm_s.props, p->polydata);
+    double dps = p->polydata->dps;
+    double s0 = p->s0;
+    double c0 = p->c0;
+    double rhor = p->polydata->rhor;
+    double porosity = p->porosity;
+    double cmax0 = p->cmax0;
+    double ads0 = p->polydata->adsorbtion(std::max(c0, cmax0));
+    double ads = p->polydata->adsorbtion(std::max(c, cmax0));
+    return  (s - dps)*c - (s0 - dps)*c0
+	+ rhor*((1.0 - porosity)/porosity)*(ads - ads0)
+	+ p->dtpv*(p->outflux*ff*mc + p->influx_polymer);
 }
 
 static struct ParametersSRes
@@ -190,12 +218,49 @@ get_parameters_s(struct PolymerSolverData *d, int cell)
 static struct ParametersCRes
 get_parameters_c(struct PolymerSolverData *d, int cell, NonlinearSolverCtrl* ctrl)
 {
-    ParametersCRes prm;
-    prm.psdata = d;
-    prm.cell = cell;
-    prm.ctrl = ctrl;
-    prm.s = -1e100;
-    return prm;
+    int i;
+    struct UnstructuredGrid *g  = d->grid;
+    double flux;
+    int f, other;
+
+    ParametersCRes p;
+    p.c0      = d->concentration[cell];
+    p.cmax0   = d->cmax[cell];
+    p.s0      = d->saturation[cell];
+    p.dtpv    = d->dt/d->porevolume[cell];
+    p.influx  = d->source[cell] >  0 ? -d->source[cell] : 0.0;
+    p.influx_polymer  = d->source[cell] >  0 ? -d->source[cell] : 0.0; // TODO. Wrong if nonzero, mult by mc.
+    for (i=g->cell_facepos[cell]; i<g->cell_facepos[cell+1]; ++i) {
+        f = g->cell_faces[i];
+
+        /* Compute cell flux*/
+        if (cell == g->face_cells[2*f]) {
+            flux  = d->darcyflux[f];
+            other = g->face_cells[2*f+1];
+        }
+        else {
+            flux  =-d->darcyflux[f];
+            other = g->face_cells[2*f];
+        }
+
+        if (other != -1) {
+            if (flux < 0.0) {
+                p.influx  += flux*d->fractionalflow[other];
+		p.influx_polymer += flux*d->fractionalflow[other]*d->mc[other];
+            }
+            else {
+                p.outflux += flux;
+            }
+        }
+    }
+    p.outflux = d->source[cell] <= 0 ? -d->source[cell] : 0.0;
+    p.porosity = d->porosity[cell];
+    p.psdata = d;
+    p.cell = cell;
+    p.ctrl = ctrl;
+    p.s = -1e100;
+    p.polydata = d->polydata;
+    return p;
 }
 
 
@@ -223,6 +288,24 @@ static double fluxfun_props(double s, double c, int cell,
     return mob[0]/(mob[0] + mob[1]);
 }
 
+static double compute_mc(double c,
+			 const Opm::IncompPropertiesInterface* props,
+			 const PolymerData* pd)
+{
+    const double* visc = props->viscosity();
+    double c_max_limit = pd->c_max_limit;
+    double cbar = c/c_max_limit;
+    double mu_w = visc[0];
+    double mu_m = pd->viscMult(c)*mu_w;
+    double mu_p = pd->viscMult(pd->c_max_limit)*mu_w;
+    double omega = pd->omega;
+    double mu_m_omega = std::pow(mu_m, omega);
+    double mu_w_e   = mu_m_omega*std::pow(mu_w, 1.0 - omega);
+    double mu_p_eff = mu_m_omega*std::pow(mu_p, 1.0 - omega);
+    double inv_mu_w_eff = (1.0 - cbar)/mu_w_e + cbar/mu_p_eff;
+    return c/(inv_mu_w_eff*mu_p_eff);
+}
+
 
 /* Local Variables:    */
 /* c-basic-offset:4    */
diff --git a/opm/polymer/polymermodel.hpp b/opm/polymer/polymermodel.hpp
index cf507913b..3a89f6f7b 100644
--- a/opm/polymer/polymermodel.hpp
+++ b/opm/polymer/polymermodel.hpp
@@ -18,11 +18,19 @@ struct PolymerData
     double omega;
     double viscMult(double c) const
     {
-	return Opm::linearInterpolation(c_vals, visc_mult_vals, c);
+	return Opm::linearInterpolation(c_vals_visc, visc_mult_vals, c);
+    }
+    double rhor;
+    double dps;
+    double adsorbtion(double c) const
+    {
+	return Opm::linearInterpolation(c_vals_ads, ads_vals, c);
     }
 
-    std::vector<double> c_vals;
+    std::vector<double> c_vals_visc;
     std::vector<double> visc_mult_vals;
+    std::vector<double> c_vals_ads;
+    std::vector<double> ads_vals;
 };
 
 
@@ -32,11 +40,14 @@ struct PolymerSolverData  {
     const PolymerData* polydata;
     const double            *darcyflux;   /* one flux per face  in cdata::grid*/
     const double            *porevolume;  /* one volume per cell */
+    const double            *porosity;
     const double            *source;      /* one source per cell */
     double                   dt;
     double                  *saturation;      /* one per cell */
     double                  *concentration;   /* one per cell */
+    double                  *cmax;
     double                  *fractionalflow;  /* one per cell */
+    double                  *mc;
 };
 
 struct NonlinearSolverCtrl;
@@ -54,10 +65,12 @@ init_solverdata(struct UnstructuredGrid *grid,
 		const PolymerData* polydata,
 		const double *darcyflux,
                 const double *porevolume,
+		const double *porosity,
 		const double *source,
                 const double dt,
 		double *saturation,
-		double *concentration);
+		double *concentration,
+		double *cmax);
 
 
 #endif /* POLYMER_H_INCLUDED */
diff --git a/opm/polymer/polymertransport.cpp b/opm/polymer/polymertransport.cpp
index e4e9aca45..890bf4532 100644
--- a/opm/polymer/polymertransport.cpp
+++ b/opm/polymer/polymertransport.cpp
@@ -16,6 +16,7 @@
 
 void polymertransport(
     const double            *porevolume,
+    const double            *porosity,
     const double            *source,
     double                   dt,
     struct UnstructuredGrid *grid,
@@ -23,7 +24,8 @@ void polymertransport(
     const PolymerData* polydata,
     const double            *darcyflux,
     double                  *saturation,
-    double                  *concentration)
+    double                  *concentration,
+    double                  *cmax)
 {
     int    i;
 
@@ -32,8 +34,8 @@ void polymertransport(
     int    *components;
 
     PolymerSolverData *data = init_solverdata(grid, props, polydata, darcyflux,
-					      porevolume, source, dt, saturation,
-					      concentration);
+					      porevolume, porosity, source, dt, saturation,
+					      concentration, cmax);
 
 
     struct NonlinearSolverCtrl ctrl;
diff --git a/opm/polymer/polymertransport.hpp b/opm/polymer/polymertransport.hpp
index 6f8f6d74d..8e6f0a731 100644
--- a/opm/polymer/polymertransport.hpp
+++ b/opm/polymer/polymertransport.hpp
@@ -14,6 +14,7 @@ struct PolymerData;
 
 void polymertransport(
     const double *porevolume,
+    const double *porosity,
     const double *source,
     double dt,
     struct UnstructuredGrid *grid,
@@ -21,7 +22,8 @@ void polymertransport(
     const PolymerData* polydata,
     const double *darcyflux,
     double *saturation,
-    double *concentration);
+    double *concentration,
+    double *cmax);
 
 
 #endif /* POLYMERTRANSPORT_HPP_INCLUDED */