diff --git a/CEP/BB/BBSKernel/include/BBSKernel/MNS/MeqParm.h b/CEP/BB/BBSKernel/include/BBSKernel/MNS/MeqParm.h
index 17700d2a9af2a7abdb528f5149ebce51ca697dae..3539ecf886fa3f032c56857b14b21987c8264140 100644
--- a/CEP/BB/BBSKernel/include/BBSKernel/MNS/MeqParm.h
+++ b/CEP/BB/BBSKernel/include/BBSKernel/MNS/MeqParm.h
@@ -79,14 +79,19 @@ public:
                  const MeqDomain&);
 
   // Initialize the solvable parameter for the given domain.
-  virtual int initDomain (const vector<MeqDomain>&, int& pertIndex,
-              vector<int>& scidIndex);
+  virtual void initDomain (const vector<MeqDomain>&);
 
   // Make parameter solvable, thus perturbed values have to be calculated.
   // spidIndex is the index of the first spid of this parm.
   // It returns the number of spids in this parm.
-  void setSolvable (bool solvable)
-    { itsIsSolvable = solvable; }
+  virtual int setSolvable(int)
+  {
+    itsIsSolvable = true;
+    return 0;
+  }
+  
+  virtual void unsetSolvable()
+  { itsIsSolvable = false; }
 
   bool isSolvable() const
     { return itsIsSolvable; }
@@ -151,11 +156,12 @@ public:
   void fillFunklets (const std::map<std::string,LOFAR::ParmDB::ParmValueSet>& parmSet,
              const MeqDomain& domain)
     { itsParmPtr->fillFunklets (parmSet, domain); }
-  int initDomain (const vector<MeqDomain>& solveDomains, int& pertIndex,
-              vector<int>& scidIndex)
-    { return itsParmPtr->initDomain (solveDomains, pertIndex, scidIndex); }
-  void setSolvable (bool solvable)
-    { itsParmPtr->setSolvable (solvable); }
+  void initDomain (const vector<MeqDomain>& solveDomains)
+    { itsParmPtr->initDomain (solveDomains); }
+  int setSolvable(int spid)
+  { return itsParmPtr->setSolvable(spid); }
+  void unsetSolvable()
+  { itsParmPtr->unsetSolvable(); }
   bool isSolvable() const
     { return itsParmPtr->isSolvable(); }
   const vector<MeqFunklet*>& getFunklets() const
diff --git a/CEP/BB/BBSKernel/include/BBSKernel/MNS/MeqParmFunklet.h b/CEP/BB/BBSKernel/include/BBSKernel/MNS/MeqParmFunklet.h
index 47bdfa2ae4047eeea57d1b60e00cf99719c53d0c..6ecc790b21d6d2cb68612807f243448a562e4d46 100644
--- a/CEP/BB/BBSKernel/include/BBSKernel/MNS/MeqParmFunklet.h
+++ b/CEP/BB/BBSKernel/include/BBSKernel/MNS/MeqParmFunklet.h
@@ -83,14 +83,16 @@ public:
   // Add a funklet.
   void add (const MeqFunklet& funklet);
 
+  virtual int setSolvable(int spid);
+  virtual void unsetSolvable();
+
   // Initialize the solvable parameter for the given solve domain size.
   // FillFunklets must have been done before.
   // The only important thing about the solve domain is its size, not the
   // start and end values. However, it is checked if the start of the
   // solve domain matches the start of the work domain given to fillFunklets.
   // It returns the number of spids used for this parm.
-  virtual int initDomain (const vector<MeqDomain>&, int& pertIndex,
-			  vector<int>& scidIndex);
+  virtual void initDomain (const vector<MeqDomain>&);
 
   // Get the current funklets.
   virtual const vector<MeqFunklet*>& getFunklets() const;
diff --git a/CEP/BB/BBSKernel/include/BBSKernel/Prediffer.h b/CEP/BB/BBSKernel/include/BBSKernel/Prediffer.h
index b01d5616a03622834c8ca005f3cfb9f7d7e9fd1e..c612249d0333e5cb28f4588ac985266d2590d756 100644
--- a/CEP/BB/BBSKernel/include/BBSKernel/Prediffer.h
+++ b/CEP/BB/BBSKernel/include/BBSKernel/Prediffer.h
@@ -121,9 +121,9 @@ public:
     
     // (Distributed) solving.
     // <group>
-    bool setParameterSelection(const vector<string> &include,
+    bool setSolvableParameters(const vector<string> &include,
         const vector<string> &exclude);        
-    void clearParameterSelection();
+    void resetSolvableParameters();
 
     // Set the global cell grid.
     bool setCellGrid(const Grid &cellGrid);
@@ -263,7 +263,7 @@ private:
     Location                                    itsStartCell, itsEndCell;
 
     //# Parameters selected for solving.
-    vector<MeqPExpr>                            itsParameterSelection;
+    vector<MeqPExpr>                            itsSolvableParameters;
     //# Local coefficient index.
     CoeffIndex                                  itsCoeffIndex;
     //# Mapping from global to local coefficients.
diff --git a/CEP/BB/BBSKernel/src/MNS/MeqParm.cc b/CEP/BB/BBSKernel/src/MNS/MeqParm.cc
index 7d3e8d618e05bcfac7a104fb9d4d83129301f26f..f216a135b19e59ec95614d955abcdd4243678a15 100644
--- a/CEP/BB/BBSKernel/src/MNS/MeqParm.cc
+++ b/CEP/BB/BBSKernel/src/MNS/MeqParm.cc
@@ -57,10 +57,8 @@ const vector<MeqFunklet*>& MeqParm::getFunklets() const
   return vec;
 }
 
-int MeqParm::initDomain (const vector<MeqDomain>&, int&, vector<int>&)
+void MeqParm::initDomain (const vector<MeqDomain>&)
 {
-  ASSERT (! isSolvable());
-  return 0;
 }
 
 void MeqParm::save()
diff --git a/CEP/BB/BBSKernel/src/MNS/MeqParmFunklet.cc b/CEP/BB/BBSKernel/src/MNS/MeqParmFunklet.cc
index 8f8cc9cd9878bb11f82850276540d93362fd5695..73fd4b2d857ef6b285d7a07a3fe3c17e585d429a 100644
--- a/CEP/BB/BBSKernel/src/MNS/MeqParmFunklet.cc
+++ b/CEP/BB/BBSKernel/src/MNS/MeqParmFunklet.cc
@@ -289,82 +289,100 @@ void MeqParmFunklet::fillFunklets
   itsWorkDomain = domain;
 }
 
-int MeqParmFunklet::initDomain (const vector<MeqDomain>& solveDomains,
-				int& pertIndex, vector<int>& scidIndex)
+
+int MeqParmFunklet::setSolvable(int spid)
 {
-  itsNrPert  = 0;
-  itsPertInx = -1;
-  if (isSolvable()) {
-    // See if the values of the previous domain can be used.
-    // We only do that if the default was used in a single previous domain.
-    ///    if (itsDefUsed  &&  polcs.size() == 1) {
-    ///      polcs[0].setDomain (domain);
-    ///      polcs[0].setNewParm();
-    ///    } else {
+    ASSERT(itsNrPert == 0 && itsPertInx == -1);
+    ASSERT(!itsFunklets.empty());
+
+    // Call the base class.
+    MeqParm::setSolvable(spid);
+    
+    // Set coefficient index and determine number of coefficients.
+    itsPertInx = spid;
+    itsNrPert = itsFunklets[0]->makeSolvable(itsPertInx);
+
+    for(uint i = 1; i < itsFunklets.size(); ++i)
+    {
+        int nrs = itsFunklets[i]->makeSolvable(itsPertInx);
+        ASSERT(nrs == itsNrPert);
+    }
+    
+    return itsNrPert;    
+}
+
+
+void MeqParmFunklet::unsetSolvable()
+{
+    // Call the base class.
+    MeqParm::unsetSolvable();
+    
+    itsPertInx = -1;
+    itsNrPert = 0;
+    for (uint i=0; i<itsFunklets.size(); i++)
+    {
+        itsFunklets[i]->clearSolvable();
+    }
+}
+
 
+void MeqParmFunklet::initDomain(const vector<MeqDomain>& solveDomains)
+{
     uint nDomain = solveDomains.size();
-    ASSERTSTR (nDomain > 0 &&
-	       solveDomains[0].endX()>itsWorkDomain.startX() &&
-	       solveDomains[0].endY()>itsWorkDomain.startY() &&
-	       solveDomains[nDomain-1].startX()<itsWorkDomain.endX() &&
-	       solveDomains[nDomain-1].startY()<itsWorkDomain.endY(),
-	       "MeqParmFunklet::initDomain of parm " << getName() << " - "
-	       "solvedomains " << solveDomains[0] << solveDomains[nDomain-1]
-	       << " mismatch workdomain " << itsWorkDomain
-	       << " given in last fillFunklets");
+
+    ASSERTSTR(nDomain > 0 &&
+        solveDomains[0].endX()>itsWorkDomain.startX() &&
+        solveDomains[0].endY()>itsWorkDomain.startY() &&
+        solveDomains[nDomain-1].startX()<itsWorkDomain.endX() &&
+        solveDomains[nDomain-1].startY()<itsWorkDomain.endY(),
+        "MeqParmFunklet::initDomain of parm " << getName() << " - "
+        "solvedomains " << solveDomains[0] << solveDomains[nDomain-1]
+        << " mismatch workdomain " << itsWorkDomain
+        << " given in last fillFunklets");
+
     // The nr of solve domains should match the nr of funklets.
     // However, if a default value is used, only one funklet is in use.
     // In that case more funklets are created if needed.
-//    if (itsDefUsed  &&  nDomain > 1  &&  itsFunklets.size() != nDomain) {
-    if (itsDefUsed) {
-      ASSERT (itsFunklets.size() == 1);
-      itsFunklets[0]->setDomain (solveDomains[0]);
-      itsFunklets.resize (nDomain);
-      for (uint i=1; i<nDomain; ++i) {
-        itsFunklets[i] = itsFunklets[0]->clone();
-	    itsFunklets[i]->setDomain (solveDomains[i]);
-      }
-      itsDefUsed = false;
-    }
+    if(itsDefUsed)
+    {
+        ASSERT(itsFunklets.size() == 1);
 
-    ASSERTSTR (itsFunklets.size() == nDomain,
-	       "Solvable parameter " << getName() << " has "
-	       << itsFunklets.size() << " funklet domains mismatching the "
-	       << nDomain << " solve domains for work domain "
-	       << itsWorkDomain);
-    for (uint i=0; i<nDomain; ++i) {
-      const MeqDomain& polDom = itsFunklets[i]->domain();
-      const MeqDomain& solDom = solveDomains[i];
-      ASSERTSTR (near(solDom.startX(), polDom.startX())  &&
-		 near(solDom.endX(), polDom.endX())  &&
-		 near(solDom.startY(), polDom.startY())  &&
-		 near(solDom.endY(), polDom.endY()),
-		 "Solvable parameter " << getName() <<
-		 " has a partially instead of fully matching entry for "
-		 "solve domain " << solDom << endl
-		 << '(' << polDom << ')');
+        itsFunklets[0]->setDomain(solveDomains[0]);
+        itsFunklets.resize(nDomain);
+        for(uint i=1; i<nDomain; ++i)
+        {
+            itsFunklets[i] = itsFunklets[0]->clone();
+            itsFunklets[i]->setDomain(solveDomains[i]);
+        }
+        
+        itsDefUsed = false;
     }
-    // Determine the solvable coeff ids for each funklet.
-    // The highest nr of scids is the nr of perturbed values.
-    itsPertInx = pertIndex;
-    for (uint i=0; i<nDomain; ++i) {
-//      int nrs = itsFunklets[i]->makeSolvable (scidIndex[i]);
-      int nrs = itsFunklets[i]->makeSolvable (itsPertInx);
-      scidIndex[i] += nrs;
-      if (nrs > itsNrPert) {
-      	itsNrPert = nrs;
-      }
-    }
-//    itsPertInx = pertIndex;
-    pertIndex += itsNrPert;
-  } else {
-    for (uint i=0; i<itsFunklets.size(); i++) {
-      itsFunklets[i]->clearSolvable();
-    }
-  }    
-  return itsNrPert;
+    else
+    {
+        ASSERTSTR (itsFunklets.size() == nDomain,
+            "Solvable parameter " << getName() << " has "
+            << itsFunklets.size() << " funklet domains mismatching the "
+            << nDomain << " solve domains for work domain "
+            << itsWorkDomain);
+        
+        for(uint i=0; i<nDomain; ++i)
+        {
+            const MeqDomain& polDom = itsFunklets[i]->domain();
+            const MeqDomain& solDom = solveDomains[i];
+
+            ASSERTSTR (near(solDom.startX(), polDom.startX())  &&
+                near(solDom.endX(), polDom.endX())  &&
+                near(solDom.startY(), polDom.startY())  &&
+                near(solDom.endY(), polDom.endY()),
+                "Solvable parameter " << getName() <<
+                " has a partially instead of fully matching entry for "
+                "solve domain " << solDom << endl
+                << '(' << polDom << ')');
+        }
+    }        
 }
 
+
 const vector<MeqFunklet*>& MeqParmFunklet::getFunklets() const
 {
   return itsFunklets;
diff --git a/CEP/BB/BBSKernel/src/Prediffer.cc b/CEP/BB/BBSKernel/src/Prediffer.cc
index 1953c4e0a25bc1152d0645cfb9b4adfa325cdb04..bdf7ace799c4edb0caa2d4091c65ca22a5d009c2 100644
--- a/CEP/BB/BBSKernel/src/Prediffer.cc
+++ b/CEP/BB/BBSKernel/src/Prediffer.cc
@@ -179,12 +179,13 @@ void Prediffer::resetState()
     itsStartCell = Location();
     itsEndCell = Location();
     itsCoeffIndex.clear();
-    clearParameterSelection();
     itsCoeffMapping.clear();    
 
+    // Reset solvable state of all parameters.
+    resetSolvableParameters();
+
     // Reset model equations.
     itsModel->clearEquations();
-//    itsParameters.clear();
 }
 
 
@@ -401,7 +402,7 @@ void Prediffer::setModelConfig(OperationType operation,
 }   
 
 
-bool Prediffer::setParameterSelection(const vector<string> &include,
+bool Prediffer::setSolvableParameters(const vector<string> &include,
         const vector<string> &exclude)
 {
     ASSERT(itsOperation == CONSTRUCT);
@@ -427,7 +428,7 @@ bool Prediffer::setParameterSelection(const vector<string> &include,
 
     // Clear the solvable flag of any parameter marked for solving and empty
     // the parameter selection.
-    clearParameterSelection();
+    resetSolvableParameters();
     
     // Find all parameters matching context.unknowns, exclude those that
     // match context.excludedUnknowns.
@@ -461,28 +462,27 @@ bool Prediffer::setParameterSelection(const vector<string> &include,
 
                 if(includeFlag)
                 {
-                    parm_it->second.setSolvable(true);
-                    itsParameterSelection.push_back(parm_it->second);
+                    itsSolvableParameters.push_back(parm_it->second);
                 }
                 
                 break;
             }
         }
     }
-    
-    return true;
+
+    return !itsSolvableParameters.empty();
 }
 
 
-void Prediffer::clearParameterSelection()
+void Prediffer::resetSolvableParameters()
 {
-    itsParameterSelection.clear();
+    itsSolvableParameters.clear();
 
     for (MeqParmGroup::iterator it = itsParameters.begin();
         it != itsParameters.end();
         ++it)
     {
-        it->second.setSolvable(false);
+        it->second.unsetSolvable();
     }
 }
 
@@ -579,7 +579,8 @@ bool Prediffer::setCellGrid(const Grid &cellGrid)
     LOG_DEBUG(oss.str());
     // DEBUG DEBUG
 
-    // Initialize model parameters marked for solving.
+    // Initialize model parameters selected for solving and construct
+    // coefficient index.
     vector<MeqDomain> domains;
     for(size_t t = itsStartCell.second; t <= itsEndCell.second; ++t)
     {
@@ -591,29 +592,20 @@ bool Prediffer::setCellGrid(const Grid &cellGrid)
         }
     }
 
-    // Clear coefficient index.
     itsCoeffIndex.clear();
-
-    // Temporary variables (required by MeqParmFunklet::initDomain() but not
-    // needed here).
-    int nCoeff = 0;
-    vector<int> tmp(nCells, 0);
-    for(MeqParmGroup::iterator it = itsParameters.begin();
-        it != itsParameters.end();
+    for(vector<MeqPExpr>::iterator it = itsSolvableParameters.begin();
+        it != itsSolvableParameters.end();
         ++it)
     {
-        // Determine maximal number of coefficients over all solution
-        // domains for this parameter.
-        int count = it->second.initDomain(domains, nCoeff, tmp);
-        
-        if(it->second.isSolvable())
-        {
-            itsCoeffIndex.insert(it->first, count);
-        }
+        // NOTE: initDomains() should come before setSolvable(), because
+        // it can modify the number of funklets of a Parm.
+        it->initDomain(domains);
+
+        int nCoeff = it->setSolvable(itsCoeffIndex.getCoeffCount());
+        itsCoeffIndex.insert(it->getName(), nCoeff);
     }
-    ASSERT(itsCoeffIndex.getCoeffCount() == static_cast<size_t>(nCoeff));
 
-    LOG_DEBUG_STR("nCoeff: " << itsCoeffIndex.getCoeffCount());
+//    LOG_DEBUG_STR("CoeffIndex: " << endl << itsCoeffIndex);
     return (itsCoeffIndex.getCoeffCount() > 0);
 }
 
@@ -629,10 +621,10 @@ void Prediffer::setCoeffIndex(const CoeffIndex &global)
 {
     ASSERT(itsOperation == CONSTRUCT);
 
-    itsCoeffMapping.resize(itsParameterSelection.size());
-    for(size_t i = 0; i < itsParameterSelection.size(); ++i)
+    itsCoeffMapping.resize(itsSolvableParameters.size());
+    for(size_t i = 0; i < itsSolvableParameters.size(); ++i)
     {
-        const MeqPExpr &parm = itsParameterSelection[i];
+        const MeqPExpr &parm = itsSolvableParameters[i];
         
         // Note: log(N) look-up.
         CoeffIndex::const_iterator it = global.find(parm.getName());
@@ -653,9 +645,9 @@ void Prediffer::getCoeff(const Location &cell, vector<double> &coeff) const
 
     // Copy the coefficient values for every selected parameter.
     coeff.reserve(itsCoeffIndex.getCoeffCount());
-    for(size_t i = 0; i < itsParameterSelection.size(); ++i)
+    for(size_t i = 0; i < itsSolvableParameters.size(); ++i)
     {
-        const MeqPExpr &parm = itsParameterSelection[i];
+        const MeqPExpr &parm = itsSolvableParameters[i];
         ASSERT(localId < parm.getFunklets().size());
 
         const MeqFunklet *funklet = parm.getFunklets()[localId];
@@ -732,9 +724,9 @@ void Prediffer::setCoeff(const Location &cell, const vector<double> &coeff)
     const size_t localId = getLocalCellId(cell);
     LOG_DEBUG_STR("Local cell id: " << localId);
 
-    for(size_t j = 0; j < itsParameterSelection.size(); ++j)
+    for(size_t j = 0; j < itsSolvableParameters.size(); ++j)
     {
-        MeqPExpr &parm = itsParameterSelection[j];
+        MeqPExpr &parm = itsSolvableParameters[j];
         parm.update(localId, coeff, itsCoeffMapping[j].start);
     }
 }
@@ -1703,6 +1695,14 @@ void Prediffer::initPolarizationMap()
 
 void Prediffer::loadParameterValues()
 {
+    // Remove the funklets from all parameters.
+    for(MeqParmGroup::iterator it = itsParameters.begin();
+        it != itsParameters.end();
+        ++it)
+    {
+        it->second.removeFunklets();
+    }
+
     // Clear all cached parameter values.
     itsParameterValues.clear();
 
@@ -1714,13 +1714,6 @@ void Prediffer::loadParameterValues()
 
     itsInstrumentDb.getValues(itsParameterValues, vector<string>(), domain);
     itsSkyDb.getValues(itsParameterValues, vector<string>(), domain);
-
-    for(MeqParmGroup::iterator it = itsParameters.begin();
-        it != itsParameters.end();
-        ++it)
-    {
-        it->second.removeFunklets();
-    }
 }
 
 
@@ -1731,9 +1724,9 @@ void Prediffer::storeParameterValues()
         itsSkyDb.lock(true);
         itsInstrumentDb.lock(true);
         
-        for(size_t i = 0; i < itsParameterSelection.size(); ++i)
+        for(size_t i = 0; i < itsSolvableParameters.size(); ++i)
         {
-            itsParameterSelection[i].save();
+            itsSolvableParameters[i].save();
         }
 
         itsSkyDb.unlock();