diff --git a/doc/source/_templates/layout.html b/doc/source/_templates/layout.html
index 5eb0f6f9ea..d180ef1e1e 100644
--- a/doc/source/_templates/layout.html
+++ b/doc/source/_templates/layout.html
@@ -104,11 +104,12 @@
      style="height:1in;"> A Python toolkit for computational methods and simulation</p>
     <div class="btn-group" role="group">
       <a class="btn btn-success"
-          href="https://github.com/erdc-cm/proteus/archive/1.4.1.tar.gz">
+          href="https://github.com/erdc-cm/proteus/archive/1.4.2.tar.gz">
             <span class="glyphicon glyphicon-download"></span> Download
       </a>
       <button type="button" class="btn btn-success dropdown-toggle" data-toggle="dropdown"><span class="caret"></span></button>
       <ul class="dropdown-menu">
+        <li><a https://github.com/erdc-cm/proteus/archive/1.4.2.tar.gz">Version 1.4.2</a></li>
         <li><a https://github.com/erdc-cm/proteus/archive/1.4.1.tar.gz">Version 1.4.1</a></li>
         <li><a https://github.com/erdc-cm/proteus/archive/1.4.0.tar.gz">Version 1.4.0</a></li>
         <li><a https://github.com/erdc-cm/proteus/archive/1.3.3.tar.gz">Version 1.3.3</a></li>
diff --git a/proteus/ModelFactory.h b/proteus/ModelFactory.h
index 209bfe8258..4ec98ca606 100644
--- a/proteus/ModelFactory.h
+++ b/proteus/ModelFactory.h
@@ -2,6 +2,20 @@
 #define MODELFACTORY_H
 #include <iostream>
 
+#define NO_INSTANCE std::cout<<"Constructing model object from template class:"<<std::endl \
+  <<"return static_cast<Model_Base*>(new ModelTemplate<CompKernelTemplate<" \
+  <<nSpaceIn<<","							\
+  <<nDOF_mesh_trial_elementIn<<"," \
+  <<nDOF_trial_elementIn<<"," \
+  <<nDOF_test_elementIn<<">," \
+  <<nSpaceIn<<"," \
+  <<nQuadraturePoints_elementIn<<"," \
+  <<nDOF_mesh_trial_elementIn<<"," \
+  <<nDOF_trial_elementIn<<"," \
+  <<nDOF_test_elementIn<<"," \
+  <<nQuadraturePoints_elementBoundaryIn<<">());" \
+  <<std::endl<<std::flush;
+
 namespace proteus
 {
   template<class Model_Base, 
@@ -26,22 +40,6 @@ namespace proteus
 						int nQuadraturePoints_elementBoundaryIn,
 						int CompKernelFlag)//0=Parametric
     {
-      std::cout<<"Constructing model object from template class:"<<std::endl
-	       <<"return static_cast<Model_Base*>(new ModelTemplate<CompKernelTemplate<"
-	       <<nSpaceIn<<","
-	       <<nDOF_mesh_trial_elementIn<<","
-	       <<nDOF_trial_elementIn<<","
-	       <<nDOF_test_elementIn<<">,"
-	       <<nSpaceIn<<","
-	       <<nQuadraturePoints_elementIn<<","
-	       <<nDOF_mesh_trial_elementIn<<","
-	       <<nDOF_trial_elementIn<<","
-	       <<nDOF_test_elementIn<<","
-	       <<nQuadraturePoints_elementBoundaryIn<<">());"
-	       <<std::endl<<std::flush;
-
-      std::cout<<"this is chooseAndAllocateDiscretization"<<std::endl;
-      
       if (CompKernelFlag == 0)
 	{
 	  if (nSpaceIn == 3)
@@ -55,7 +53,7 @@ namespace proteus
 	  	      else if (nQuadraturePoints_elementIn == 4)
 	  		return static_cast<Model_Base*>(new ModelTemplate<CompKernelTemplate<3,4,4,4>,3,4,4,4,4,3>());
 	  	      else
-	  		abort();
+	  		NO_INSTANCE
 	  	    }
 	  	  else if (nDOF_mesh_trial_elementIn == 8)
 	  	    {
@@ -65,24 +63,24 @@ namespace proteus
 	  		else if ( nQuadraturePoints_elementBoundaryIn == 9)
 	  		  return static_cast<Model_Base*>(new ModelTemplate<CompKernelTemplate<3,8,8,8>,3,8,8,8,8,9>());
 	  		else
-	  		  abort();
+	  		  NO_INSTANCE
 	  	      else if (nQuadraturePoints_elementIn == 27)
 	  		if ( nQuadraturePoints_elementBoundaryIn == 4)
 	  		  return static_cast<Model_Base*>(new ModelTemplate<CompKernelTemplate<3,8,8,8>,3,27,8,8,8,4>());
 	  		else if ( nQuadraturePoints_elementBoundaryIn == 9)
 	  		  return static_cast<Model_Base*>(new ModelTemplate<CompKernelTemplate<3,8,8,8>,3,27,8,8,8,9>());
 	  		else
-	  		  abort();
+	  		  NO_INSTANCE
 	  	      else
 			
-	  		abort();
+	  		NO_INSTANCE
 	  	    }
 	  	  else if (nDOF_mesh_trial_elementIn == 10)
 	  	    {
 	  	      if (nQuadraturePoints_elementIn == 15)
 	  		return static_cast<Model_Base*>(new ModelTemplate<CompKernelTemplate<3,10,10,10>,3,15,10,10,10,7>());
 	  	      else
-	  		abort();
+	  		NO_INSTANCE
 	  	    }
 	  	  else if (nDOF_mesh_trial_elementIn == 27)
 	  	    {
@@ -92,19 +90,19 @@ namespace proteus
 	  		else if ( nQuadraturePoints_elementBoundaryIn == 9)
 	  		  return static_cast<Model_Base*>(new ModelTemplate<CompKernelTemplate<3,27,27,27>,3,8,27,27,27,9>());
 	  		else
-	  		  abort();
+	  		  NO_INSTANCE
 	  	      else if (nQuadraturePoints_elementIn == 27)
 	  		if ( nQuadraturePoints_elementBoundaryIn == 4)
 	  		         return static_cast<Model_Base*>(new ModelTemplate<CompKernelTemplate<3,27,27,27>,3,27,27,27,27,4>());
 	  		else if ( nQuadraturePoints_elementBoundaryIn == 9)
 	  		  return static_cast<Model_Base*>(new ModelTemplate<CompKernelTemplate<3,27,27,27>,3,27,27,27,27,9>());
 	  		else
-	  		  abort();
+	  		  NO_INSTANCE
 	  	      else
-	  		abort();
+	  		NO_INSTANCE
 	  	    }
 	  	  else
-	  	    abort();
+	  	    NO_INSTANCE
 	  	}
 	      else if (nDOF_mesh_trial_elementIn == 4)//sub-parametric tets
 	  	{
@@ -136,11 +134,11 @@ namespace proteus
 	  		}
 	  	      else
 	  		{
-	  		  abort();
+	  		  NO_INSTANCE
 	  		}
 	  	    }
 	  	  else
-	  	    abort();
+	  	    NO_INSTANCE
 	  	}
 	      else if (nDOF_mesh_trial_elementIn == 8)//sub-parametric hexes
 	  	{
@@ -149,13 +147,13 @@ namespace proteus
 	  	      return static_cast<Model_Base*>(new ModelTemplate<CompKernelTemplate<3,8,27,27>,3,125,8,27,27,25>());
 	  	}
 	      else
-	  	abort();
+	  	NO_INSTANCE
 	    }
 	  else
-	    abort();
+	    NO_INSTANCE
 	}
       else
-	abort();	       
+	NO_INSTANCE	       
       return NULL;   
     }	  		      
   template<class Model_Base, 
@@ -180,130 +178,78 @@ namespace proteus
 						  int nQuadraturePoints_elementBoundaryIn,
 						  int CompKernelFlag)//0=Parametric
     {
-      std::cout<<"Constructing model object from template class:"<<std::endl 
-	       <<"return static_cast<Model_Base*>(new ModelTemplate<CompKernelTemplate<" 
-	       <<nSpaceIn<<"," 
-	       <<nDOF_mesh_trial_elementIn<<"," 
-	       <<nDOF_trial_elementIn<<"," 
-	       <<nDOF_test_elementIn<<">," 
-	       <<nSpaceIn<<"," 
-	       <<nQuadraturePoints_elementIn<<"," 
-	       <<nDOF_mesh_trial_elementIn<<"," 
-	       <<nDOF_trial_elementIn<<"," 
-	       <<nDOF_test_elementIn<<"," 
-	       <<nQuadraturePoints_elementBoundaryIn<<">());" 
-	       <<std::endl<<std::flush;
-    
-      std::cout<<"this is chooseAndAllocateDiscretization2D"
-          <<std::endl<<std::flush;
-    
-
       if (CompKernelFlag == 0)	  
 	{
 	  if (nSpaceIn == 2)
 	    {
 	      if (nDOF_mesh_trial_elementIn == nDOF_trial_elementIn)//iso-parametric
 		{
-                        std::cout<<"I am here"
-                            <<std::endl<<std::flush;
-        
 		  if (nDOF_mesh_trial_elementIn == 3)
 		    {
-                std::cout<<"3333333"
-                            <<std::endl<<std::flush;
-
-                
 		      if (nQuadraturePoints_elementIn == 1)
 		      	{
-                        std::cout<<"111111"
-                            <<std::endl<<std::flush;
-                    
 		      	  if (nQuadraturePoints_elementBoundaryIn == 1)
 		      	    return static_cast<Model_Base*>(new ModelTemplate<CompKernelTemplate<2,3,3,3>,2,1,3,3,3,1>());
 		      	  else
-		      	    abort();
+		      	    NO_INSTANCE
 		      	}
 		      else if (nQuadraturePoints_elementIn == 3)
 		      	{
-                        std::cout<<"3333333"
-                            <<std::endl<<std::flush;
-                    
 		      	  if (nQuadraturePoints_elementBoundaryIn == 2)
 		      	    return static_cast<Model_Base*>(new ModelTemplate<CompKernelTemplate<2,3,3,3>,2,3,3,3,3,2>());
 		      	  else
-		      	    abort();
+		      	    NO_INSTANCE
 		      	}
 		      else if (nQuadraturePoints_elementIn == 4)
 		      	{
-                        std::cout<<"444444444"
-                            <<std::endl<<std::flush;
-                    
 		      	  if (nQuadraturePoints_elementBoundaryIn == 3)
 		      	    return static_cast<Model_Base*>(new ModelTemplate<CompKernelTemplate<2,3,3,3>,2,4,3,3,3,3>());
 		      	  else
-		      	    abort();
+		      	    NO_INSTANCE
 		      	}
 		      else if (nQuadraturePoints_elementIn == 6)
 		      	{
-                        std::cout<<"666666666"
-                            <<std::endl<<std::flush;
-                    
 		      	  if (nQuadraturePoints_elementBoundaryIn == 4)
 		      	    return static_cast<Model_Base*>(new ModelTemplate<CompKernelTemplate<2,3,3,3>,2,6,3,3,3,4>());
 		      	  else
-		      	    abort();
+		      	    NO_INSTANCE
 		      	}
 		      else if (nQuadraturePoints_elementIn == 7)
 		      	{
-                        std::cout<<"777777777"
-                            <<std::endl<<std::flush;
-                    
 		      	  if (nQuadraturePoints_elementBoundaryIn == 5)
 		      	    return static_cast<Model_Base*>(new ModelTemplate<CompKernelTemplate<2,3,3,3>,2,7,3,3,3,5>());
 		      	  else
-		      	    abort();
+		      	    NO_INSTANCE
 		      	}
-              else if (nQuadraturePoints_elementIn == 16)//hk=0.5
-                {
-                        std::cout<<"16-16-16-16-16-16-16"
-                            <<std::endl<<std::flush;
-                    
-                  if (nQuadraturePoints_elementBoundaryIn == 3)
-                    {
-                        std::cout<<"I am here"
-                            <<std::endl<<std::flush;
-		      	    return static_cast<Model_Base*>(new ModelTemplate<CompKernelTemplate<2,3,3,3>,2,16,3,3,3,3>());
-                    }else
-		      	        abort();
-                }
-              else if (nQuadraturePoints_elementIn == 36)//hk=0.3
-                {
-                        std::cout<<"36-36-36-36-36-36-36"
-                            <<std::endl<<std::flush;
-                    
-                  if (nQuadraturePoints_elementBoundaryIn == 3)
-                    {
-                        std::cout<<"I am here"
-                            <<std::endl<<std::flush;
-		      	    return static_cast<Model_Base*>(new ModelTemplate<CompKernelTemplate<2,3,3,3>,2,36,3,3,3,3>());
-                    }else
-		      	        abort();
-                }
-              else if (nQuadraturePoints_elementIn == 100)//hk=0.2
-                {
-                        std::cout<<"100-100-100-100-100"
-                            <<std::endl<<std::flush;
-                    
-                  if (nQuadraturePoints_elementBoundaryIn == 3)
-                    {
-                        std::cout<<"I am here"
-                            <<std::endl<<std::flush;
-		      	    return static_cast<Model_Base*>(new ModelTemplate<CompKernelTemplate<2,3,3,3>,2,100,3,3,3,3>());
-                    }else
-		      	        abort();
-                }
+		      else if (nQuadraturePoints_elementIn == 16)//hk=0.5
+			{
+			  if (nQuadraturePoints_elementBoundaryIn == 3)
+			    {
+			      return static_cast<Model_Base*>(new ModelTemplate<CompKernelTemplate<2,3,3,3>,2,16,3,3,3,3>());
+			    }
+			  else
+			    NO_INSTANCE
+			}
+		      else if (nQuadraturePoints_elementIn == 36)//hk=0.3
+			{
+			  if (nQuadraturePoints_elementBoundaryIn == 3)
+			    {
+			      return static_cast<Model_Base*>(new ModelTemplate<CompKernelTemplate<2,3,3,3>,2,36,3,3,3,3>());
+			    }
+			  else
+			    NO_INSTANCE
+			}
+		      else if (nQuadraturePoints_elementIn == 100)//hk=0.2
+			{
+			  if (nQuadraturePoints_elementBoundaryIn == 3)
+			    {
+			      return static_cast<Model_Base*>(new ModelTemplate<CompKernelTemplate<2,3,3,3>,2,100,3,3,3,3>());
+			    }
+			  else
+			    NO_INSTANCE
+			}
 		      else
-			        abort();
+			NO_INSTANCE
 		    }
                   else  if(nDOF_mesh_trial_elementIn == 4)
                     {
@@ -312,13 +258,13 @@ namespace proteus
                           if (nQuadraturePoints_elementBoundaryIn == 2)
                             return static_cast<Model_Base*>(new ModelTemplate<CompKernelTemplate<2,4,4,4>,2,4,4,4,4,2>());
                           else
-                            abort();
+                            NO_INSTANCE
                         }
                       else
-                        abort();
+                        NO_INSTANCE
                     }
 		  else
-		    abort();
+		    NO_INSTANCE
 		}
 	      else if (nDOF_mesh_trial_elementIn == 3)
 		{
@@ -329,28 +275,28 @@ namespace proteus
 		      	  if (nQuadraturePoints_elementBoundaryIn == 1)
 		      	    return static_cast<Model_Base*>(new ModelTemplate<CompKernelTemplate<2,3,6,6>,2,1,3,6,6,1>());
 		      	  else
-		      	    abort();
+		      	    NO_INSTANCE
 		      	}
 		      else if (nQuadraturePoints_elementIn == 3)
 		      	{
 		      	  if (nQuadraturePoints_elementBoundaryIn == 2)
 		      	    return static_cast<Model_Base*>(new ModelTemplate<CompKernelTemplate<2,3,6,6>,2,3,3,6,6,2>());
 		      	  else
-		      	    abort();
+		      	    NO_INSTANCE
 		      	}
 		      else if (nQuadraturePoints_elementIn == 4)
 		      	{
 		      	  if (nQuadraturePoints_elementBoundaryIn == 3)
 		      	    return static_cast<Model_Base*>(new ModelTemplate<CompKernelTemplate<2,3,6,6>,2,4,3,6,6,3>());
 		      	  else
-		      	    abort();
+		      	    NO_INSTANCE
 		      	}
 		      else if (nQuadraturePoints_elementIn == 6)
 		      	{
 		      	  if (nQuadraturePoints_elementBoundaryIn == 4)
 		      	    return static_cast<Model_Base*>(new ModelTemplate<CompKernelTemplate<2,3,6,6>,2,6,3,6,6,4>());
 		      	  else
-		      	    abort();
+		      	    NO_INSTANCE
 		      	}
 		      else if (nQuadraturePoints_elementIn == 7)
 		      	{
@@ -359,21 +305,21 @@ namespace proteus
 		      	      return static_cast<Model_Base*>(new ModelTemplate<CompKernelTemplate<2,3,6,6>,2,7,3,6,6,5>());
 		      	    }
 		      	  else
-		      	    abort();
+		      	    NO_INSTANCE
 		      	}
 		      else
-			abort();
+			NO_INSTANCE
 		    }
 		  else
-		    abort();
+		    NO_INSTANCE
 		}
 	      else
-		abort();
+		NO_INSTANCE
 	    }	  		      
 	}
       else
         {
-          abort();
+          NO_INSTANCE
         }
       return NULL;
     }
diff --git a/proteus/SimTools.py b/proteus/SimTools.py
index baae4a7f1f..2334c38116 100644
--- a/proteus/SimTools.py
+++ b/proteus/SimTools.py
@@ -43,6 +43,7 @@ class SimulationProcessor:
                     'dataDir' :'.',              #where file is located
                     'appendResults':False,       #append to existing data files?
                     'echo':False,                #print to screen
+                    'echoRelativeErrors':False,  #print to screen also relative errors
                     'components':[0],            #list of components to monitor
                     'errorQuantities':[None],    #quantities in which to estimate error
                     'errorNorms':[None],         #norms to use for error calc
@@ -779,8 +780,11 @@ def uOfXT(self,X,T):
                             self.errorData[ci][il][kerr].append(err)
                             self.errorData[ci][il][kexa].append(exa)
                             if self.flags['echo']:
-                                logEvent("""\nt= %g; %s[%d][%d]= %g;""" % (tsim,kerr,ci,il,err),level=0)
-                            #end if
+                                if self.flags['echoRelativeErrors']:
+                                    logEvent("""\nt= %g; %s[%d][%d]= %g; relative_error= %g;""" % (tsim,kerr,ci,il,err,err/(exa+1E-15)),level=0)
+                                else:
+                                    logEvent("""\nt= %g; %s[%d][%d]= %g;""" % (tsim,kerr,ci,il,err),level=0)
+                                #end if
                         #if calcL2u
                         calcL1u = ('L1' in self.flags['errorNorms'] and
                                    'u' in self.flags['errorQuantities'])
@@ -804,13 +808,15 @@ def uOfXT(self,X,T):
                                 err = Norms.L1errorSFEM(mFine.q[('dV_u',ci)][0:mFine.mesh.subdomainMesh.nElements_owned],udense[0:mFine.mesh.subdomainMesh.nElements_owned],uproj[ci][0:mFine.mesh.subdomainMesh.nElements_owned])
                                 exa = Norms.L1errorSFEM(mFine.q[('dV_u',ci)][0:mFine.mesh.subdomainMesh.nElements_owned],udense[0:mFine.mesh.subdomainMesh.nElements_owned],
                                                         numpy.zeros(udense[0:mFine.mesh.subdomainMesh.nElements_owned].shape,'d'))
-
                             kerr = 'error_'+'u'+'_'+'L1'
                             kexa = 'exact_'+'u'+'_'+'L1'
                             self.errorData[ci][il][kerr].append(err)
                             self.errorData[ci][il][kexa].append(exa)
                             if self.flags['echo']:
-                                logEvent("""\nt= %g; %s[%d][%d]= %g;""" % (tsim,kerr,ci,il,err),level=0)
+                                if self.flags['echoRelativeErrors']:
+                                    logEvent("""\nt= %g; %s[%d][%d]= %g; relative_error= %g;""" % (tsim,kerr,ci,il,err,err/(exa+1E-15)),level=0)
+                                else:
+                                    logEvent("""\nt= %g; %s[%d][%d]= %g;""" % (tsim,kerr,ci,il,err),level=0)                                
                             #end if
                         #if calcL1u
                         calcLIu = ('LI' in self.flags['errorNorms'] and
@@ -837,7 +843,10 @@ def uOfXT(self,X,T):
                             self.errorData[ci][il][kerr].append(err)
                             self.errorData[ci][il][kexa].append(exa)
                             if self.flags['echo']:
-                                logEvent("""\nt= %g; %s[%d][%d]= %g;""" % (tsim,kerr,ci,il,err),level=0)
+                                if self.flags['echoRelativeErrors']:
+                                    logEvent("""\nt= %g; %s[%d][%d]= %g; relative_error= %g;""" % (tsim,kerr,ci,il,err,err/(exa+1E-15)),level=0)
+                                else:
+                                    logEvent("""\nt= %g; %s[%d][%d]= %g;""" % (tsim,kerr,ci,il,err),level=0)
                             #end if
                             #
                         #calcLIu
@@ -894,7 +903,10 @@ def uOfXT(self,X,T):
                             self.errorData[ci][il][kerr].append(err)
                             self.errorData[ci][il][kexa].append(exa)
                             if self.flags['echo']:
-                                logEvent("""\nt= %g; %s[%d][%d]= %g;""" % (tsim,kerr,ci,il,err),level=0)
+                                if self.flags['echoRelativeErrors']:
+                                    logEvent("""\nt= %g; %s[%d][%d]= %g; relative_error= %g;""" % (tsim,kerr,ci,il,err,err/(exa+1E-15)),level=0)
+                                else:
+                                    logEvent("""\nt= %g; %s[%d][%d]= %g;""" % (tsim,kerr,ci,il,err),level=0)
                             #end if
                         #if calcH1u
                         calcH1semiU = ('H1semi' in self.flags['errorNorms'] and
@@ -928,7 +940,10 @@ def uOfXT(self,X,T):
                             self.errorData[ci][il][kerr].append(err)
                             self.errorData[ci][il][kexa].append(exa)
                             if self.flags['echo']:
-                                logEvent("""\nt= %g; %s[%d][%d]= %g;""" % (tsim,kerr,ci,il,err),level=0)
+                                if self.flags['echoRelativeErrors']:
+                                    logEvent("""\nt= %g; %s[%d][%d]= %g; relative_error= %g;""" % (tsim,kerr,ci,il,err,err/(exa+1E-15)),level=0)
+                                else:
+                                    logEvent("""\nt= %g; %s[%d][%d]= %g;""" % (tsim,kerr,ci,il,err),level=0)
                             #end if
                         #if calcH1semiu
                         calcW11u = ('W11' in self.flags['errorNorms'] and
@@ -983,7 +998,10 @@ def uOfXT(self,X,T):
                             self.errorData[ci][il][kerr].append(err)
                             self.errorData[ci][il][kexa].append(exa)
                             if self.flags['echo']:
-                                logEvent("""\nt= %g; %s[%d][%d]= %g;""" % (tsim,kerr,ci,il,err),level=0)
+                                if self.flags['echoRelativeErrors']:
+                                    logEvent("""\nt= %g; %s[%d][%d]= %g; relative_error= %g;""" % (tsim,kerr,ci,il,err,err/(exa+1E-15)),level=0)
+                                else:
+                                    logEvent("""\nt= %g; %s[%d][%d]= %g;""" % (tsim,kerr,ci,il,err),level=0)
                             #end if
                         #if calcW11u
                         calcW11semiU = ('W11semi' in self.flags['errorNorms'] and
@@ -1016,7 +1034,10 @@ def uOfXT(self,X,T):
                             self.errorData[ci][il][kerr].append(err)
                             self.errorData[ci][il][kexa].append(exa)
                             if self.flags['echo']:
-                                logEvent("""\nt= %g; %s[%d][%d]= %g;""" % (tsim,kerr,ci,il,err),level=0)
+                                if self.flags['echoRelativeErrors']:
+                                    logEvent("""\nt= %g; %s[%d][%d]= %g; relative_error= %g;""" % (tsim,kerr,ci,il,err,err/(exa+1E-15)),level=0)
+                                else:
+                                    logEvent("""\nt= %g; %s[%d][%d]= %g;""" % (tsim,kerr,ci,il,err),level=0)
                             #end if
                         #if calcH1semiu
                         calcTVu = ('TV' in self.flags['errorNorms'] and
@@ -1053,7 +1074,10 @@ def uOfXT(self,X,T):
                             self.errorData[ci][il][kerr].append(err)
                             self.errorData[ci][il][kexa].append(exa)
                             if self.flags['echo']:
-                                logEvent("""\nt= %g; %s[%d][%d]= %g;""" % (tsim,kerr,ci,il,err),level=0)
+                                if self.flags['echoRelativeErrors']:
+                                    logEvent("""\nt= %g; %s[%d][%d]= %g; relative_error= %g;""" % (tsim,kerr,ci,il,err,err/(exa+1E-15)),level=0)
+                                else:
+                                    logEvent("""\nt= %g; %s[%d][%d]= %g;""" % (tsim,kerr,ci,il,err),level=0)
                             #end if
                         #end calcTV
                         ############### velocity specific calculations ###############
@@ -1090,7 +1114,10 @@ def uOfXT(self,X,T):
                             self.errorData[ci][il][kerr].append(err)
                             self.errorData[ci][il][kexa].append(exa)
                             if self.flags['echo']:
-                                logEvent("""\nt= %g; %s[%d][%d]= %g;""" % (tsim,kerr,ci,il,err),level=0)
+                                if self.flags['echoRelativeErrors']:
+                                    logEvent("""\nt= %g; %s[%d][%d]= %g; relative_error= %g;""" % (tsim,kerr,ci,il,err,err/(exa+1E-15)),level=0)
+                                else:
+                                    logEvent("""\nt= %g; %s[%d][%d]= %g;""" % (tsim,kerr,ci,il,err),level=0)
                             #end if
                         #if calcL2vel
                         calcL1vel = ('L1' in self.flags['errorNorms'] and
@@ -1124,7 +1151,10 @@ def uOfXT(self,X,T):
                             self.errorData[ci][il][kerr].append(err)
                             self.errorData[ci][il][kexa].append(exa)
                             if self.flags['echo']:
-                                logEvent("""\nt= %g; %s[%d][%d]= %g;""" % (tsim,kerr,ci,il,err),level=0)
+                                if self.flags['echoRelativeErrors']:
+                                    logEvent("""\nt= %g; %s[%d][%d]= %g; relative_error= %g;""" % (tsim,kerr,ci,il,err,err/(exa+1E-15)),level=0)
+                                else:
+                                    logEvent("""\nt= %g; %s[%d][%d]= %g;""" % (tsim,kerr,ci,il,err),level=0)
                             #end if
                         #if calcL2vel
 
@@ -1156,7 +1186,10 @@ def uOfXT(self,X,T):
                             self.errorData[ci][il][kerr].append(err)
                             self.errorData[ci][il][kexa].append(exa)
                             if self.flags['echo']:
-                                logEvent("""\nt= %g; %s[%d][%d]= %g;""" % (tsim,kerr,ci,il,err),level=0)
+                                if self.flags['echoRelativeErrors']:
+                                    logEvent("""\nt= %g; %s[%d][%d]= %g; relative_error= %g;""" % (tsim,kerr,ci,il,err,err/(exa+1E-15)),level=0)
+                                else:
+                                    logEvent("""\nt= %g; %s[%d][%d]= %g;""" % (tsim,kerr,ci,il,err),level=0)
                             #end if
                         #if calcLIvel
 
diff --git a/proteus/__init__.py b/proteus/__init__.py
index 6728068d14..5c51e7e7b2 100644
--- a/proteus/__init__.py
+++ b/proteus/__init__.py
@@ -9,7 +9,7 @@
     import pkgutil
     __path__ = pkgutil.extend_path(__path__, __name__)
 
-__version__ = '1.4.1.dev0'
+__version__ = '1.4.2'
 
 __all__ = ["Archiver",
            "Domain",
diff --git a/proteus/mprans/SW2DCV.h b/proteus/mprans/SW2DCV.h
index 2d29cdd1d2..3d38b4f575 100644
--- a/proteus/mprans/SW2DCV.h
+++ b/proteus/mprans/SW2DCV.h
@@ -13,7 +13,6 @@
 //5. Try other choices of variables h,hu,hv, Bova-Carey symmetrization?
 
 #define POWER_SMOOTHNESS_INDICATOR 2
-#define LINEAR_FRICTION 0
 #define VEL_FIX_POWER 2.
 #define REESTIMATE_MAX_EDGE_BASED_CFL 1
 
@@ -54,9 +53,9 @@ namespace proteus
 			    int NNZ, //number on non-zero entries on sparsity pattern
 			    int numDOFs, //number of DOFs
 			    double* lumped_mass_matrix, //lumped mass matrix (as vector)
-			    double* h_lstage, //DOFs of solution at last stage
-			    double* hu_lstage, 
-			    double* hv_lstage, 
+			    double* h_old, //DOFs of solution at last stage
+			    double* hu_old, 
+			    double* hv_old, 
 			    double* b_dof,
 			    double* high_order_hnp1, //DOFs of high order solution at tnp1
 			    double* high_order_hunp1, 
@@ -79,9 +78,9 @@ namespace proteus
       virtual double calculateEdgeBasedCFL(double g, 
 					   int numDOFsPerEqn, //number of DOFs
 					   double* lumped_mass_matrix, //lumped mass matrix (as vector)
-					   double* h_lstage, //DOFs of solution at last stage
-					   double* hu_lstage, 
-					   double* hv_lstage, 
+					   double* h_old, //DOFs of solution at last stage
+					   double* hu_old, 
+					   double* hv_old, 
 					   double* b_dof,
 					   int* csrRowIndeces_DofLoops, //csr row indeces 
 					   int* csrColumnOffsets_DofLoops, //csr column offsets 
@@ -138,12 +137,6 @@ namespace proteus
 				   double* h_dof_old,
  				   double* hu_dof_old, 
 				   double* hv_dof_old,
-				   double* h_dof_old_old,
- 				   double* hu_dof_old_old, 
-				   double* hv_dof_old_old,
-				   double* h_dof_lstage,
- 				   double* hu_dof_lstage, 
-				   double* hv_dof_lstage,
 				   double* b_dof,
 				   double* h_dof, 
 				   double* hu_dof, 
@@ -235,6 +228,7 @@ namespace proteus
 				   double cE, 
 				   int LUMPED_MASS_MATRIX, 
 				   double dt,
+				   int LINEAR_FRICTION,
 				   double mannings, 
 				   double* quantDOFs,
 				   int SECOND_CALL_CALCULATE_RESIDUAL,
@@ -289,12 +283,6 @@ namespace proteus
 				   double* h_dof_old,
 				   double* hu_dof_old, 
 				   double* hv_dof_old,
-				   double* h_dof_old_old, 
-				   double* hu_dof_old_old, 
-				   double* hv_dof_old_old,
-				   double* h_dof_lstage, 
-				   double* hu_dof_lstage, 
-				   double* hv_dof_lstage,
 				   double* b_dof,
 				   double* h_dof, 
 				   double* hu_dof, 
@@ -386,6 +374,7 @@ namespace proteus
 				   double cE,
 				   int LUMPED_MASS_MATRIX,
 				   double dt,
+				   int LINEAR_FRICTION,
 				   double mannings,
 				   // Quant of interests
 				   double* quantDOFs,
@@ -1800,9 +1789,9 @@ namespace proteus
 		    int NNZ, //number on non-zero entries on sparsity pattern
 		    int numDOFs, //number of DOFs
 		    double* lumped_mass_matrix, //lumped mass matrix (as vector))
-		    double* h_lstage, //DOFs of solution at last stage
-		    double* hu_lstage, 
-		    double* hv_lstage,
+		    double* h_old, //DOFs of solution at last stage
+		    double* hu_old, 
+		    double* hv_old,
 		    double* b_dof,
 		    double* high_order_hnp1, //DOFs of high order solution at tnp1
 		    double* high_order_hunp1, 
@@ -1832,7 +1821,7 @@ namespace proteus
 	{
 	  //read some vectors 
 	  double high_order_hnp1i  = high_order_hnp1[i];
-	  double hni = h_lstage[i];
+	  double hni = h_old[i];
 	  double Zi = b_dof[i];
 	  double mi = lumped_mass_matrix[i];
 
@@ -1843,7 +1832,7 @@ namespace proteus
 	    {
 	      int j = csrColumnOffsets_DofLoops[offset];
 	      // read some vectors
-	      double hnj = h_lstage[j];
+	      double hnj = h_old[j];
 	      double Zj = b_dof[j];
 
 	      // COMPUTE STAR SOLUTION // hStar, huStar and hvStar
@@ -1887,9 +1876,9 @@ namespace proteus
 	  double high_order_hnp1i  = high_order_hnp1[i];
 	  double high_order_hunp1i = high_order_hunp1[i];
 	  double high_order_hvnp1i = high_order_hvnp1[i];
-	  double hni = h_lstage[i];
-	  double huni = hu_lstage[i];
-	  double hvni = hv_lstage[i];
+	  double hni = h_old[i];
+	  double huni = hu_old[i];
+	  double hvni = hv_old[i];
 	  double Zi = b_dof[i];
 	  double mi = lumped_mass_matrix[i];
 	  double one_over_hiReg = 2*hni/(hni*hni+std::pow(fmax(hni,hEps),2)); //hEps
@@ -1903,9 +1892,9 @@ namespace proteus
 	    {
 	      int j = csrColumnOffsets_DofLoops[offset];
 	      // read some vectors
-	      double hnj = h_lstage[j];
-	      double hunj = hu_lstage[j];
-	      double hvnj = hv_lstage[j];
+	      double hnj = h_old[j];
+	      double hunj = hu_old[j];
+	      double hvnj = hv_old[j];
 	      double Zj = b_dof[j];
 	      double one_over_hjReg = 2*hnj/(hnj*hnj+std::pow(fmax(hnj,hEps),2)); //hEps
 
@@ -1971,9 +1960,9 @@ namespace proteus
     double calculateEdgeBasedCFL(double g, 
 				 int numDOFsPerEqn, //number of DOFs
 				 double* lumped_mass_matrix, //lumped mass matrix (as vector))
-				 double* h_dof_lstage, //DOFs of solution at last stage
-				 double* hu_dof_lstage, 
-				 double* hv_dof_lstage,
+				 double* h_dof_old, //DOFs of solution at last stage
+				 double* hu_dof_old, 
+				 double* hv_dof_old,
 				 double* b_dof,
 				 int* csrRowIndeces_DofLoops, //csr row indeces 
 				 int* csrColumnOffsets_DofLoops, //csr column offsets 
@@ -1992,9 +1981,9 @@ namespace proteus
       int ij=0;
       for (int i=0; i<numDOFsPerEqn; i++)
 	{
-	  double hi = h_dof_lstage[i]; // solution at time tn for the ith DOF
-	  double hui = hu_dof_lstage[i]; 
-	  double hvi = hv_dof_lstage[i]; 
+	  double hi = h_dof_old[i]; // solution at time tn for the ith DOF
+	  double hui = hu_dof_old[i]; 
+	  double hvi = hv_dof_old[i]; 
 	  double dLowii = 0.;
 
 	  double alphai;
@@ -2003,9 +1992,9 @@ namespace proteus
 	  for (int offset=csrRowIndeces_DofLoops[i]; offset<csrRowIndeces_DofLoops[i+1]; offset++)
 	    { //loop in j (sparsity pattern)
 	      int j = csrColumnOffsets_DofLoops[offset];
-	      double hj = h_dof_lstage[j]; // solution at time tn for the jth DOF
-	      double huj = hu_dof_lstage[j];
-	      double hvj = hv_dof_lstage[j];
+	      double hj = h_dof_old[j]; // solution at time tn for the jth DOF
+	      double huj = hu_dof_old[j];
+	      double hvj = hv_dof_old[j];
 
 	      if (i != j) 
 		{
@@ -2039,7 +2028,7 @@ namespace proteus
 	  // CALCULATE EDGE BASED CFL //
 	  //////////////////////////////
 	  double mi = lumped_mass_matrix[i];
-	  edge_based_cfl[i] = fabs(dLowii)/mi;
+	  edge_based_cfl[i] = 2*fabs(dLowii)/mi;
 	  max_edge_based_cfl = fmax(max_edge_based_cfl,edge_based_cfl[i]);
 
 	  //////////////////////////////////
@@ -2067,9 +2056,9 @@ namespace proteus
 	  ij=0;
 	  for (int i=0; i<numDOFsPerEqn; i++)
 	    {
-	      double hi = h_dof_lstage[i]; // solution at time tn for the ith DOF
-	      double hui = hu_dof_lstage[i]; 
-	      double hvi = hv_dof_lstage[i]; 
+	      double hi = h_dof_old[i]; // solution at time tn for the ith DOF
+	      double hui = hu_dof_old[i]; 
+	      double hvi = hv_dof_old[i]; 
 	      double Zi = b_dof[i];
 	      double one_over_hiReg = 2*hi/(hi*hi+std::pow(fmax(hi,hEps),2)); // hEps
 	      double ui = hui*one_over_hiReg;
@@ -2082,9 +2071,9 @@ namespace proteus
 	      for (int offset=csrRowIndeces_DofLoops[i]; offset<csrRowIndeces_DofLoops[i+1]; offset++)
 		{
 		  int j = csrColumnOffsets_DofLoops[offset];
-		  double hj = h_dof_lstage[j]; // solution at time tn for the jth DOF
-		  double huj = hu_dof_lstage[j];
-		  double hvj = hv_dof_lstage[j];
+		  double hj = h_dof_old[j]; // solution at time tn for the jth DOF
+		  double huj = hu_dof_old[j];
+		  double hvj = hv_dof_old[j];
 		  double Zj = b_dof[j];	      
 		  double one_over_hjReg = 2*hj/(hj*hj+std::pow(fmax(hj,hEps),2)); //hEps
 		  double uj = huj*one_over_hjReg;
@@ -2173,12 +2162,6 @@ namespace proteus
 			   double* h_dof_old,
 			   double* hu_dof_old, 
 			   double* hv_dof_old,
-			   double* h_dof_old_old, 
-			   double* hu_dof_old_old, 
-			   double* hv_dof_old_old, 
-			   double* h_dof_lstage, 
-			   double* hu_dof_lstage, 
-			   double* hv_dof_lstage, 
 			   double* b_dof, 
 			   double* h_dof, 
 			   double* hu_dof, 
@@ -2269,7 +2252,8 @@ namespace proteus
 			   double* muH_minus_muL,
 			   double cE,
 			   int LUMPED_MASS_MATRIX,
-			   double dt, 
+			   double dt,
+			   int LINEAR_FRICTION,
 			   double mannings,
 			   // Quant of interests
 			   double* quantDOFs, 
@@ -3158,12 +3142,6 @@ namespace proteus
 			   double* h_dof_old,
 			   double* hu_dof_old, 
 			   double* hv_dof_old,
-			   double* h_dof_old_old, 
-			   double* hu_dof_old_old, 
-			   double* hv_dof_old_old, 
-			   double* h_dof_lstage, 
-			   double* hu_dof_lstage, 
-			   double* hv_dof_lstage, 
 			   double* b_dof, 
 			   double* h_dof, 
 			   double* hu_dof, 
@@ -3255,6 +3233,7 @@ namespace proteus
 			   double cE,
 			   int LUMPED_MASS_MATRIX,
 			   double dt,
+			   int LINEAR_FRICTION,
 			   double mannings,
 			   // Quant of interests
 			   double* quantDOFs,
@@ -3269,8 +3248,10 @@ namespace proteus
       //FOR FRICTION//
       double n2 = std::pow(mannings,2.);
       double gamma=4./3;
-      double xi=2.;
-
+      // mql. This parameter relaxes the cfl restriction.
+      // It is now conservative. I might change it after the local bounds are implemented 
+      double xi=10.; 
+      
       //////////////////////////////////////
       // ********** CELL LOOPS ********** //
       //////////////////////////////////////
@@ -3303,7 +3284,7 @@ namespace proteus
 		eN_nDOF_trial_element = eN*nDOF_trial_element;
 	      register double 
 		h=0.0,hu=0.0,hv=0.0, // solution at current time
-		h_lstage=0.0,hu_lstage=0.0,hv_lstage=0.0, // solution at lstage
+		h_old=0.0,hu_old=0.0,hv_old=0.0, // solution at lstage
 		jac[nSpace*nSpace], jacDet, jacInv[nSpace*nSpace],
 		h_test_dV[nDOF_trial_element],
 		dV,x,y,xt,yt;
@@ -3325,15 +3306,15 @@ namespace proteus
 	      ck.valFromDOF(hu_dof,&vel_l2g[eN_nDOF_trial_element],&vel_trial_ref[k*nDOF_trial_element],hu);
 	      ck.valFromDOF(hv_dof,&vel_l2g[eN_nDOF_trial_element],&vel_trial_ref[k*nDOF_trial_element],hv);
 	      // get the solution at the lstage
-	      ck.valFromDOF(h_dof_lstage,&h_l2g[eN_nDOF_trial_element],&h_trial_ref[k*nDOF_trial_element],h_lstage);
-	      ck.valFromDOF(hu_dof_lstage,&vel_l2g[eN_nDOF_trial_element],&vel_trial_ref[k*nDOF_trial_element],hu_lstage);
-	      ck.valFromDOF(hv_dof_lstage,&vel_l2g[eN_nDOF_trial_element],&vel_trial_ref[k*nDOF_trial_element],hv_lstage);
+	      ck.valFromDOF(h_dof_old,&h_l2g[eN_nDOF_trial_element],&h_trial_ref[k*nDOF_trial_element],h_old);
+	      ck.valFromDOF(hu_dof_old,&vel_l2g[eN_nDOF_trial_element],&vel_trial_ref[k*nDOF_trial_element],hu_old);
+	      ck.valFromDOF(hv_dof_old,&vel_l2g[eN_nDOF_trial_element],&vel_trial_ref[k*nDOF_trial_element],hv_old);
 	      // calculate cell based CFL to keep a reference
 	      calculateCFL(elementDiameter[eN],
 			   g,
-			   h_lstage,
-			   hu_lstage,
-			   hv_lstage,
+			   h_old,
+			   hu_old,
+			   hv_old,
 			   hEps,
 			   q_cfl[eN_k]);
 	      //precalculate test function products with integration weights
@@ -3349,9 +3330,9 @@ namespace proteus
 	      for(int i=0;i<nDOF_test_element;i++)
 		{
 		  // compute time derivative part of global residual. NOTE: no lumping
-		  elementResidual_h[i]  += (h  - h_lstage)*h_test_dV[i];
-		  elementResidual_hu[i] += (hu - hu_lstage)*h_test_dV[i];
-		  elementResidual_hv[i] += (hv - hv_lstage)*h_test_dV[i];
+		  elementResidual_h[i]  += (h  - h_old)*h_test_dV[i];
+		  elementResidual_hu[i] += (hu - hu_old)*h_test_dV[i];
+		  elementResidual_hv[i] += (hv - hv_old)*h_test_dV[i];
 		}
 	    }
 	  // distribute
@@ -3379,9 +3360,9 @@ namespace proteus
 	  for (int i=0; i<numDOFsPerEqn; i++)
 	    {
 	      // COMPUTE ENTROPY. NOTE: WE CONSIDER A FLAT BOTTOM 
-	      double hni = h_dof_lstage[i]; 
+	      double hni = h_dof_old[i]; 
 	      double one_over_hniReg = 2*hni/(hni*hni+std::pow(fmax(hni,hEps),2)); //hEps
-	      eta[i] = ENTROPY(g,hni,hu_dof_lstage[i],hv_dof_lstage[i],0.,one_over_hniReg);
+	      eta[i] = ENTROPY(g,hni,hu_dof_old[i],hv_dof_old[i],0.,one_over_hniReg);
 	    }
 	  // ********** END OF COMPUTING ENTROPY ********** //
 	  
@@ -3395,9 +3376,9 @@ namespace proteus
 	  register double psi[numDOFsPerEqn], etaMax[numDOFsPerEqn], etaMin[numDOFsPerEqn];
 	  for (int i=0; i<numDOFsPerEqn; i++)
 	    {	      
-	      double hi = h_dof_lstage[i]; // solution at time tn for the ith DOF
-	      double hui = hu_dof_lstage[i]; 
-	      double hvi = hv_dof_lstage[i]; 
+	      double hi = h_dof_old[i]; // solution at time tn for the ith DOF
+	      double hui = hu_dof_old[i]; 
+	      double hvi = hv_dof_old[i]; 
 	      double one_over_hiReg = 2*hi/(hi*hi+std::pow(fmax(hi,hEps),2)); //hEps
 
 	      // For eta min and max
@@ -3418,9 +3399,9 @@ namespace proteus
 	      for (int offset=csrRowIndeces_DofLoops[i]; offset<csrRowIndeces_DofLoops[i+1]; offset++)
 		{ //loop in j (sparsity pattern)
 		  int j = csrColumnOffsets_DofLoops[offset];
-		  double hj = h_dof_lstage[j]; // solution at time tn for the jth DOF
-		  double huj = hu_dof_lstage[j];
-		  double hvj = hv_dof_lstage[j];
+		  double hj = h_dof_old[j]; // solution at time tn for the jth DOF
+		  double huj = hu_dof_old[j];
+		  double hvj = hv_dof_old[j];
 		  double one_over_hjReg = 2*hj/(hj*hj+std::pow(fmax(hj,hEps),2)); //hEps
 		  double uj = huj*one_over_hjReg;
 		  double vj = hvj*one_over_hjReg;
@@ -3491,9 +3472,9 @@ namespace proteus
 	  ij = 0;
 	  for (int i=0; i<numDOFsPerEqn; i++)
 	    {
-	      double hi = h_dof_lstage[i];
-	      double hui = hu_dof_lstage[i];
-	      double hvi = hv_dof_lstage[i];	  
+	      double hi = h_dof_old[i];
+	      double hui = hu_dof_old[i];
+	      double hvi = hv_dof_old[i];	  
 	      double Zi = b_dof[i];
 	      double mi = lumped_mass_matrix[i];
 	      double one_over_hiReg = 2*hi/(hi*hi+std::pow(fmax(hi,hEps),2)); // hEps
@@ -3531,15 +3512,11 @@ namespace proteus
 	      double veli_norm = std::sqrt(ui*ui+vi*vi);
 	      double hi_to_the_gamma = std::pow(hi,gamma);
 	      double friction_aux = 
-		veli_norm == 0. ? 0. : 2*g*n2*veli_norm*mi/(hi_to_the_gamma+fmax(hi_to_the_gamma,10*xi*g*n2*dt*veli_norm));
+		veli_norm == 0. ? 0.: (2*g*n2*veli_norm*mi/
+				       (hi_to_the_gamma+fmax(hi_to_the_gamma,xi*g*n2*dt*veli_norm)));
 	      double ith_friction_term2 = friction_aux*hui;
 	      double ith_friction_term3 = friction_aux*hvi;
 
-	      //double ith_friction_term2 =  
-	      //veli_norm==0. ? 0. : 2*g*n2*hui*veli_norm*mi/(hi_to_the_gamma+fmax(hi_to_the_gamma,xi*g*n2*dt*veli_norm));
-	      //double ith_friction_term3 =  
-	      //veli_norm==0. ? 0. : 2*g*n2*hvi*veli_norm*mi/(hi_to_the_gamma+fmax(hi_to_the_gamma,xi*g*n2*dt*veli_norm));
-	      
 	      if (LINEAR_FRICTION==1)
 		{
 		  ith_friction_term2 = mannings*hui*mi;
@@ -3550,9 +3527,9 @@ namespace proteus
 	      for (int offset=csrRowIndeces_DofLoops[i]; offset<csrRowIndeces_DofLoops[i+1]; offset++)
 		{
 		  int j = csrColumnOffsets_DofLoops[offset];
-		  double hj = h_dof_lstage[j];
-		  double huj = hu_dof_lstage[j];
-		  double hvj = hv_dof_lstage[j];
+		  double hj = h_dof_old[j];
+		  double huj = hu_dof_old[j];
+		  double hvj = hv_dof_old[j];
 		  double Zj = b_dof[j];
 		  double one_over_hjReg = 2*hj/(hj*hj+std::pow(fmax(hj,hEps),2)); //hEps
 		  double uj = huj*one_over_hjReg;
@@ -3599,7 +3576,7 @@ namespace proteus
 								      hEps,hEps,false)*cji_norm); //hEps
 			}
 		      dLij = dLowij*fmax(psi[i],psi[j]); // enhance the order to 2nd order. No EV
-		      
+
 		      ///////////////////////////////////////
 		      // WELL BALANCING DISSIPATIVE MATRIX //
 		      ///////////////////////////////////////
diff --git a/proteus/mprans/SW2DCV.py b/proteus/mprans/SW2DCV.py
index 83a8766c31..7710e46174 100644
--- a/proteus/mprans/SW2DCV.py
+++ b/proteus/mprans/SW2DCV.py
@@ -75,7 +75,7 @@ def updateShockCapturingHistory(self):
                                                                     globalMax(self.numDiff_last[1].max()),
                                                                     globalMax(self.numDiff_last[2].max())))
 
-class RKEV(proteus.TimeIntegration.SSP33):
+class RKEV(proteus.TimeIntegration.SSP):
     from proteus import TimeIntegration
     """
     Wrapper for SSPRK time integration using EV
@@ -87,54 +87,38 @@ def __init__(self, transport, timeOrder=1, runCFL=0.1, integrateInterpolationPoi
         self.trasport = transport
         self.runCFL=runCFL
         self.dtLast=None
+        self.dtRatioMax = 2.0
         self.isAdaptive=True
         # About the cfl 
         assert hasattr(transport,'edge_based_cfl'), "No edge based cfl defined"
         self.edge_based_cfl = transport.edge_based_cfl
-        self.cell_based_cfl = transport.q[('cfl',0)]
-        # Stuff particular for SSP33
+        # Stuff particular for SSP
         self.timeOrder = timeOrder  #order of approximation
         self.nStages = timeOrder  #number of stages total
         self.lstage = 0  #last stage completed
-        # storage vectors
-        # previous time step mass and solution dof per component
-        self.m_last = {}
-        #temporarily use this to stash previous solution since m_last used
-        #in EV transport models for previous solution value
-        self.m_last_save = {}
+        # storage vectors (at old time step)
         self.u_dof_last = {}
-        # per component stage values, list with array at each stage
-        self.m_stage = {}
-        self.u_dof_stage = {}
+        # per component lstage values
+        self.u_dof_lstage = {}
         for ci in range(self.nc):
-             if transport.q.has_key(('m',ci)):
-                self.m_last[ci] = transport.q[('m',ci)].copy()
-                self.m_last_save[ci] = transport.q[('m',ci)].copy()
-
-                self.u_dof_last[ci] = transport.u[ci].dof.copy()
-                self.m_stage[ci] = []
-                self.u_dof_stage[ci] = []
-                for k in range(self.nStages+1):                    
-                    self.m_stage[ci].append(transport.q[('m',ci)].copy())
-                    self.u_dof_stage[ci].append(transport.u[ci].dof.copy())
+            self.u_dof_last[ci] = transport.u[ci].dof.copy()
+            self.u_dof_lstage[ci] = transport.u[ci].dof.copy()
         
-    #def set_dt(self, DTSET):
-    #    self.dt = DTSET #  don't update t
     def choose_dt(self):
         maxCFL = 1.0e-6
         # COMPUTE edge_based_cfl        
         rowptr_cMatrix, colind_cMatrix, Cx = self.transport.cterm_global[0].getCSRrepresentation()
-        rowptr_cMatrix, colind_cMatrix, Cy = self.transport.cterm_global[1].getCSRrepresentation()        
+        rowptr_cMatrix, colind_cMatrix, Cy = self.transport.cterm_global[1].getCSRrepresentation()
         rowptr_cMatrix, colind_cMatrix, CTx = self.transport.cterm_global_transpose[0].getCSRrepresentation()
         rowptr_cMatrix, colind_cMatrix, CTy = self.transport.cterm_global_transpose[1].getCSRrepresentation()
         numDOFsPerEqn = self.transport.u[0].dof.size
-        
+
         adjusted_maxCFL = self.transport.sw2d.calculateEdgeBasedCFL(
             self.transport.coefficients.g, 
             numDOFsPerEqn,
             self.transport.ML,
-            self.transport.u[0].dof, 
-            self.transport.u[1].dof, 
+            self.transport.u[0].dof,
+            self.transport.u[1].dof,
             self.transport.u[2].dof,
             self.transport.coefficients.b.dof,
             rowptr_cMatrix, 
@@ -150,14 +134,15 @@ def choose_dt(self):
             self.transport.edge_based_cfl)
 
         maxCFL = max(maxCFL,max(adjusted_maxCFL, globalMax(self.edge_based_cfl.max())))
-        # maxCFL = max(maxCFL,globalMax(self.cell_based_cfl.max()))
-        self.dt = self.runCFL/maxCFL            
+        self.dt = self.runCFL/maxCFL
+
         if self.dtLast is None:
             self.dtLast = self.dt
+        if self.dt/self.dtLast  > self.dtRatioMax:
+            self.dt = self.dtLast*self.dtRatioMax
         self.t = self.tLast + self.dt
-        # mwf debug
-        #print "RKEv max cfl component ci dt dtLast {0} {1} {2} {3}".format(maxCFL,ci,self.dt,self.dtLast)
         self.substeps = [self.t for i in range(self.nStages)] #Manuel is ignoring different time step levels for now
+
     def initialize_dt(self,t0,tOut,q):
         """
         Modify self.dt
@@ -171,93 +156,90 @@ def setCoefficients(self):
         beta are all 1's here
         mwf not used right now
         """
-        self.alpha = numpy.zeros((self.nStages, self.nStages),'d')
-        self.dcoefs = numpy.zeros((self.nStages),'d')
+        # Not needed for an implementation when alpha and beta are not used
         
     def updateStage(self):
         """
         Need to switch to use coefficients
         """
-        #mwf debug
-        #import pdb
-        #pdb.set_trace()
         self.lstage += 1
-        assert self.timeOrder in [1,3]
+        assert self.timeOrder in [1,2,3]
         assert self.lstage > 0 and self.lstage <= self.timeOrder
-        if self.timeOrder == 3:
+        #print "within update stage...: ", self.lstage
+        if self.timeOrder == 3:            
             if self.lstage == 1:
+                logEvent("First stage of SSP33 method finished",level=4)
                 for ci in range(self.nc):
-                    self.u_dof_stage[ci][self.lstage][:] = numpy.copy(self.transport.u[ci].dof) #no need for .copy?
-                    self.m_stage[ci][self.lstage][:] = numpy.copy(self.transport.q[('m',ci)])
-                    #needs to be updated for non-scalar equations
-                    #this as used as last stage value in EV Transport model
-                    #mwf TODO, get rid of m_last here
-                    self.m_last[ci] = numpy.copy(self.transport.q[('m',ci)])
-
+                    self.u_dof_lstage[ci][:] = self.transport.u[ci].dof
+                # update u_dof_old
+                self.transport.h_dof_old[:] = self.u_dof_lstage[0]
+                self.transport.hu_dof_old[:] = self.u_dof_lstage[1]
+                self.transport.hv_dof_old[:] = self.u_dof_lstage[2]
             elif self.lstage == 2:
+                logEvent("Second stage of SSP33 method finished",level=4)
+                for ci in range(self.nc):
+                    self.u_dof_lstage[ci][:] = self.transport.u[ci].dof
+                    self.u_dof_lstage[ci] *= 1./4.
+                    self.u_dof_lstage[ci] += 3./4.*self.u_dof_last[ci]
+                # update u_dof_old
+                self.transport.h_dof_old[:] = self.u_dof_lstage[0]
+                self.transport.hu_dof_old[:] = self.u_dof_lstage[1]
+                self.transport.hv_dof_old[:] = self.u_dof_lstage[2]
+            else:
+                logEvent("Third stage of SSP33 method finished",level=4)
                 for ci in range(self.nc):
-                    self.u_dof_stage[ci][self.lstage][:] = numpy.copy(self.transport.u[ci].dof) 
-                    self.u_dof_stage[ci][self.lstage] *= 1./4.
-                    self.u_dof_stage[ci][self.lstage] += 3./4.*self.u_dof_last[ci]
-                    self.m_stage[ci][self.lstage][:] = numpy.copy(self.transport.q[('m',ci)])
-                    self.m_stage[ci][self.lstage] *= 1./4.
-                    #mwf this has to be fixed
-                    #previous stage updated m_last to the stage value
-                    #either have another temporary here or have the VOF code use m_stage
-                    #instead of m_last
-                    self.m_stage[ci][self.lstage] += 3./4.*self.m_last_save[ci] 
-                    #mwf TODO get rid of this
-                    self.m_last[ci] = numpy.copy(self.m_stage[ci][self.lstage])
-            elif self.lstage == 3:
+                    self.u_dof_lstage[ci][:] = self.transport.u[ci].dof
+                    self.u_dof_lstage[ci][:] *= 2.0/3.0
+                    self.u_dof_lstage[ci][:] += 1.0/3.0*self.u_dof_last[ci]
+                    # update solution to u[0].dof
+                    self.transport.u[ci].dof[:] = self.u_dof_lstage[ci]
+                # update u_dof_old
+                self.transport.h_dof_old[:] = self.u_dof_last[0]
+                self.transport.hu_dof_old[:] = self.u_dof_last[1]
+                self.transport.hv_dof_old[:] = self.u_dof_last[2]
+        elif self.timeOrder == 2:
+            if self.lstage == 1:
+                logEvent("First stage of SSP22 method finished",level=4)
                 for ci in range(self.nc):
-                    self.u_dof_stage[ci][self.lstage][:] = numpy.copy(self.transport.u[ci].dof)
-                    self.u_dof_stage[ci][self.lstage][:] *= 2.0/3.0
-                    self.u_dof_stage[ci][self.lstage][:] += 1.0/3.0*self.u_dof_last[ci]
-                    #switch  time history back
-                    #mwf TODO this needs to be fixed for multipcomponent
-                    self.transport.u[ci].dof[:] = numpy.copy(self.u_dof_stage[ci][self.lstage])
-                    self.m_last[ci] = numpy.copy(self.m_last_save[ci])
-                    tmp_dof_stage = self.u_dof_stage[ci][self.lstage].copy()
-                    self.u_dof_stage[ci][self.lstage][:] = self.u_dof_last[ci]
-                    #globalResidualDummy = numpy.zeros(self.u.shape,'d')
-                    #self.transport.getResidual(tmp_dof_stage,
-                    #                           globalResidualDummy)
-                    self.u_dof_stage[ci][self.lstage][:] = tmp_dof_stage
-                    #self.setUnkowns(self.u)
-
+                    self.u_dof_lstage[ci][:] = self.transport.u[ci].dof
+                # Update u_dof_old
+                self.transport.h_dof_old[:] = self.u_dof_lstage[0]
+                self.transport.hu_dof_old[:] = self.u_dof_lstage[1]
+                self.transport.hv_dof_old[:] = self.u_dof_lstage[2]
+            else:
+                logEvent("Second stage of SSP22 method finished",level=4)
+                for ci in range(self.nc):
+                    self.u_dof_lstage[ci][:] = self.transport.u[ci].dof
+                    self.u_dof_lstage[ci][:] *= 1./2.
+                    self.u_dof_lstage[ci][:] += 1./2.*self.u_dof_last[ci]
+                    # update solution to u[0].dof
+                    self.transport.u[ci].dof[:] = self.u_dof_lstage[ci]
+                # Update u_dof_old
+                self.transport.h_dof_old[:] = self.u_dof_last[0] #HHHEEEEEREEEEE!!!
+                self.transport.hu_dof_old[:] = self.u_dof_last[1]
+                self.transport.hv_dof_old[:] = self.u_dof_last[2]                                
         else:
             assert self.timeOrder == 1
-            for ci in range(self.nc):
-                self.m_stage[ci][self.lstage][:]=self.transport.q[('m',ci)][:]
-                self.u_dof_stage[ci][self.lstage][:] = self.transport.u[ci].dof[:]
-                            
+            logEvent("FE method finished",level=4)
+            
     def initializeTimeHistory(self,resetFromDOF=True):
         """
         Push necessary information into time history arrays
         """
         for ci in range(self.nc):
-            self.m_last[ci][:] = self.transport.q[('m',ci)][:]
             self.u_dof_last[ci][:] = self.transport.u[ci].dof[:]
-            self.m_last_save[ci][:] = self.transport.q[('m',ci)][:]
-            for k in range(self.nStages):
-                self.m_stage[ci][k][:] = self.transport.q[('m',ci)][:]
-                self.u_dof_stage[ci][k][:] = self.transport.u[ci].dof[:]
+
     def updateTimeHistory(self,resetFromDOF=False):
         """
         assumes successful step has been taken
         """
-        
         self.t = self.tLast + self.dt
         for ci in range(self.nc):
-            self.m_last[ci][:] = self.transport.q[('m',ci)][:]
-            self.m_last_save[ci][:] = self.transport.q[('m',ci)][:]
             self.u_dof_last[ci][:] = self.transport.u[ci].dof[:]
-            for k in range(self.nStages):
-                self.m_stage[ci][k][:]=self.transport.q[('m',ci)][:]
-                self.u_dof_stage[ci][k][:] = self.transport.u[ci].dof[:]
         self.lstage=0
         self.dtLast = self.dt
         self.tLast = self.t
+
     def generateSubsteps(self,tList):
         """
         create list of substeps over time values given in tList. These correspond to stages
@@ -278,16 +260,11 @@ def resetOrder(self,order):
         self.lstage = 0  #last stage completed
         # storage vectors
         # per component stage values, list with array at each stage
-        self.m_stage = {}
-        self.u_dof_stage = {}
+        self.u_dof_lstage = {}
         for ci in range(self.nc):
-             if self.transport.q.has_key(('m',ci)):
-                self.m_stage[ci] = []
-                self.u_dof_stage[ci] = []
-                for k in range(self.nStages+1):                    
-                    self.m_stage[ci].append(self.transport.q[('m',ci)].copy())
-                    self.u_dof_stage[ci].append(self.transport.u[ci].dof.copy())
+            self.u_dof_lstage[ci] = self.transport.u[ci].dof.copy()
         self.substeps = [self.t for i in range(self.nStages)]            
+
     def setFromOptions(self,nOptions):
         """
         allow classes to set various numerical parameters
@@ -317,8 +294,9 @@ def __init__(self,
 		 useMetrics=0.0,
                  modelIndex=0, 
                  cE=1.0,
-                 LUMPED_MASS_MATRIX=1, 
-                 mannings=0.):
+                 LUMPED_MASS_MATRIX=1,
+                 LINEAR_FRICTION=0,
+                 mannings=0.):        
         self.bathymetry = bathymetry 
         self.useRBLES=useRBLES
         self.useMetrics=useMetrics
@@ -328,6 +306,7 @@ def __init__(self,
         self.nd=nd
         self.cE=cE
         self.LUMPED_MASS_MATRIX=LUMPED_MASS_MATRIX
+        self.LINEAR_FRICTION=LINEAR_FRICTION
         self.mannings=mannings
         self.modelIndex=modelIndex
         mass={}
@@ -399,12 +378,9 @@ def updateToMovingDomain(self,t,c):
     def evaluate(self,t,c):
         pass
     def preStep(self,t,firstStep=False):
-        self.model.h_dof_old_old = numpy.copy(self.model.h_dof_old)
-        self.model.hu_dof_old_old = numpy.copy(self.model.hu_dof_old)
-        self.model.hv_dof_old_old = numpy.copy(self.model.hv_dof_old)
-        self.model.h_dof_old = numpy.copy(self.model.u[0].dof)
-        self.model.hu_dof_old = numpy.copy(self.model.u[1].dof)
-        self.model.hv_dof_old = numpy.copy(self.model.u[2].dof)
+        self.model.h_dof_old[:] = self.model.u[0].dof
+        self.model.hu_dof_old[:] = self.model.u[1].dof
+        self.model.hv_dof_old[:] = self.model.u[2].dof
     def postStep(self,t,firstStep=False):
         pass
 
@@ -439,7 +415,6 @@ def __init__(self,
                  bdyNullSpace=False):
         self.bdyNullSpace=bdyNullSpace
         self.inf_norm_hu=[] #To test 1D well balancing
-        self.firstCalculateResidualCall=True
         self.secondCallCalculateResidual=0
         self.postProcessing = False#this is a hack to test the effect of post-processing
         #
@@ -645,12 +620,9 @@ def __init__(self,
         #NOTE (Mql): It is important to link h_dof_old by reference with u[0].dof (and so on).
         # This is because  I need the initial condition to be passed to them as well (before calling calculateResidual). 
         # During preStep I change this and copy the values instead of keeping the reference. 
-        self.h_dof_old_old = self.u[0].dof
-        self.hu_dof_old_old = self.u[1].dof
-        self.hv_dof_old_old = self.u[2].dof
-        self.h_dof_old = self.u[0].dof
-        self.hu_dof_old = self.u[1].dof
-        self.hv_dof_old = self.u[2].dof
+        self.h_dof_old = None
+        self.hu_dof_old = None
+        self.hv_dof_old = None
 
         #Vector for mass matrix
         self.check_positivity_water_height=True
@@ -922,9 +894,9 @@ def FCTStep(self):
                           self.nnz, #number of non zero entries 
                           len(rowptr)-1, #number of DOFs
                           self.ML, #Lumped mass matrix
-                          self.timeIntegration.u_dof_stage[0][self.timeIntegration.lstage], #hn
-                          self.timeIntegration.u_dof_stage[1][self.timeIntegration.lstage], #hun
-                          self.timeIntegration.u_dof_stage[2][self.timeIntegration.lstage], #hvn
+                          self.h_dof_old,
+                          self.hu_dof_old,
+                          self.hv_dof_old,
                           self.coefficients.b.dof,
                           self.timeIntegration.u[hIndex], #high order solution 
                           self.timeIntegration.u[huIndex],
@@ -954,6 +926,10 @@ def getResidual(self,u,r):
         Calculate the element residuals and add in to the global residual
         """
 
+        if self.h_dof_old is None:
+            self.h_dof_old = numpy.copy(self.u[0].dof)
+            self.hu_dof_old = numpy.copy(self.u[1].dof)
+            self.hv_dof_old = numpy.copy(self.u[2].dof)                        
         #COMPUTE hEps
         if self.hEps is None: 
             eps=1E-14
@@ -1210,12 +1186,6 @@ def getResidual(self,u,r):
         if (self.check_positivity_water_height==True):
             assert self.u[0].dof.min() >= 0, ("Negative water height: ", self.u[0].dof.min())
 
-        if (self.firstCalculateResidualCall):
-            self.timeIntegration.u_dof_stage[0][self.timeIntegration.lstage][:] = self.u[0].dof
-            self.timeIntegration.u_dof_stage[1][self.timeIntegration.lstage][:] = self.u[1].dof
-            self.timeIntegration.u_dof_stage[2][self.timeIntegration.lstage][:] = self.u[2].dof
-            self.firstCalculateResidualCall = False
-
         self.calculateResidual(
             self.u[0].femSpace.elementMaps.psi,
             self.u[0].femSpace.elementMaps.grad_psi,
@@ -1259,12 +1229,6 @@ def getResidual(self,u,r):
             self.h_dof_old,
             self.hu_dof_old,
             self.hv_dof_old,
-            self.h_dof_old_old,
-            self.hu_dof_old_old,
-            self.hv_dof_old_old,
-            self.timeIntegration.u_dof_stage[0][self.timeIntegration.lstage],
-            self.timeIntegration.u_dof_stage[1][self.timeIntegration.lstage],
-            self.timeIntegration.u_dof_stage[2][self.timeIntegration.lstage],
             self.coefficients.b.dof,
             self.u[0].dof,
             self.u[1].dof,
@@ -1349,7 +1313,8 @@ def getResidual(self,u,r):
             self.muH_minus_muL,
             self.coefficients.cE, 
             self.coefficients.LUMPED_MASS_MATRIX, 
-            self.timeIntegration.dt, 
+            self.timeIntegration.dt,
+            self.coefficients.LINEAR_FRICTION,
             self.coefficients.mannings,
             self.quantDOFs, 
             self.secondCallCalculateResidual, 
diff --git a/proteus/mprans/VOF.py b/proteus/mprans/VOF.py
index 57790d4905..e28349a33e 100644
--- a/proteus/mprans/VOF.py
+++ b/proteus/mprans/VOF.py
@@ -32,7 +32,7 @@ def updateShockCapturingHistory(self):
         if self.lag:
             for ci in range(self.nc):
                 self.numDiff_last[ci][:] = self.numDiff[ci]
-        if self.lag == False and self.nStepsToDelay != None and self.nSteps > self.nStepsToDelay:
+        if self.lag == False and self.nStepsToDelay is not None and self.nSteps > self.nStepsToDelay:
             logEvent("VOF.ShockCapturing: switched to lagged shock capturing")
             self.lag = True
             self.numDiff_last=[]
@@ -395,7 +395,7 @@ def attachModels(self,modelList):
             self.q_porosity = numpy.ones(modelList[self.modelIndex].q[('u',0)].shape,'d')
             self.porosity_dof = numpy.ones(modelList[self.modelIndex].u[0].dof.shape,'d')
 
-            if self.setParamsFunc != None:
+            if self.setParamsFunc is not None:
                 self.setParamsFunc(modelList[self.modelIndex].q['x'],self.q_porosity)
             #
         #
@@ -404,7 +404,7 @@ def attachModels(self,modelList):
         elif modelList[self.modelIndex].ebq.has_key(('u',0)):
             self.ebq_porosity = numpy.ones(modelList[self.modelIndex].ebq[('u',0)].shape,
                                            'd')
-            if self.setParamsFunc != None:
+            if self.setParamsFunc is not None:
                 self.setParamsFunc(modelList[self.modelIndex].ebq['x'],self.ebq_porosity)
             #
         #
@@ -413,7 +413,7 @@ def attachModels(self,modelList):
         else:
             self.ebqe_porosity = numpy.ones(modelList[self.LS_modelIndex].ebqe[('u',0)].shape,
                                             'd')
-            if self.setParamsFunc != None:
+            if self.setParamsFunc is not None:
                 self.setParamsFunc(modelList[self.LS_modelIndex].ebqe['x'],self.ebqe_porosity)
             #
         #
@@ -485,7 +485,7 @@ def evaluate(self,t,c):
             v = self.ebqe_v
             phi = self.ebqe_phi
             porosity  = self.ebq_porosity
-        elif ((self.ebq_v != None and self.ebq_phi != None) and c[('f',0)].shape == self.ebq_v.shape):
+        elif ((self.ebq_v is not None and self.ebq_phi is not None) and c[('f',0)].shape == self.ebq_v.shape):
             v = self.ebq_v
             phi = self.ebq_phi
             porosity  = self.ebq_porosity
@@ -493,7 +493,7 @@ def evaluate(self,t,c):
             v=None
             phi=None
             porosity=None
-        if v != None:
+        if v is not None:
             # self.VOFCoefficientsEvaluate(self.eps,
             #                              v,
             #                              phi,
@@ -591,7 +591,7 @@ def __init__(self,
         #determine whether  the stabilization term is nonlinear
         self.stabilizationIsNonlinear = False
         #cek come back
-	if self.stabilization != None:
+	if self.stabilization is not None:
 	    for ci in range(self.nc):
 		if coefficients.mass.has_key(ci):
 		    for flag in coefficients.mass[ci].values():
@@ -621,8 +621,8 @@ def __init__(self,
         #determine if we need element boundary storage
         self.elementBoundaryIntegrals = {}
         for ci  in range(self.nc):
-            self.elementBoundaryIntegrals[ci] = ((self.conservativeFlux != None) or
-                                                 (numericalFluxType != None) or
+            self.elementBoundaryIntegrals[ci] = ((self.conservativeFlux is not None) or
+                                                 (numericalFluxType is not None) or
                                                  (self.fluxBoundaryConditions[ci] == 'outFlow') or
                                                  (self.fluxBoundaryConditions[ci] == 'mixedFlow') or
                                                  (self.fluxBoundaryConditions[ci] == 'setFlow'))
@@ -655,7 +655,7 @@ def __init__(self,
         else:
             for I in self.coefficients.elementIntegralKeys:
                 elementQuadratureDict[I] = elementQuadrature
-        if self.stabilization != None:
+        if self.stabilization is not None:
             for I in self.coefficients.elementIntegralKeys:
                 if elemQuadIsDict:
                     if elementQuadrature.has_key(I):
@@ -664,7 +664,7 @@ def __init__(self,
                         elementQuadratureDict[('stab',)+I[1:]] = elementQuadrature['default']
                 else:
                     elementQuadratureDict[('stab',)+I[1:]] = elementQuadrature
-        if self.shockCapturing != None:
+        if self.shockCapturing is not None:
             for ci in self.shockCapturing.components:
                 if elemQuadIsDict:
                     if elementQuadrature.has_key(('numDiff',ci,ci)):
@@ -805,7 +805,7 @@ def __init__(self,
         else:
              self.timeIntegration = TimeIntegrationClass(self)
 
-        if options != None:
+        if options is not None:
             self.timeIntegration.setFromOptions(options)
         logEvent(memory("TimeIntegration","OneLevelTransport"),level=4)
         logEvent("Calculating numerical quadrature formulas",2)
@@ -844,10 +844,10 @@ def __init__(self,
         comm = Comm.get()
         self.comm=comm
         if comm.size() > 1:
-            assert numericalFluxType != None and numericalFluxType.useWeakDirichletConditions,"You must use a numerical flux to apply weak boundary conditions for parallel runs"
+            assert numericalFluxType is not None and numericalFluxType.useWeakDirichletConditions,"You must use a numerical flux to apply weak boundary conditions for parallel runs"
 
         logEvent(memory("stride+offset","OneLevelTransport"),level=4)
-        if numericalFluxType != None:
+        if numericalFluxType is not None:
             if options is None or options.periodicDirichletConditions is None:
                 self.numericalFlux = numericalFluxType(self,
                                                        dofBoundaryConditionsSetterDict,
@@ -966,7 +966,7 @@ def updateVelocityFieldAsFunction(self):
             self.coefficients.ebqe_v[...,2] = self.velocityFieldAsFunction[2](ebqe_X,t)
 
     def calculateElementResidual(self):
-        if self.globalResidualDummy != None:
+        if self.globalResidualDummy is not None:
             self.getResidual(self.u[0].dof,self.globalResidualDummy)
 
     def getResidual(self,u,r):
@@ -1412,10 +1412,10 @@ def calculateElementQuadrature(self):
         self.u[0].femSpace.getBasisValuesRef(self.elementQuadraturePoints)
         self.u[0].femSpace.getBasisGradientValuesRef(self.elementQuadraturePoints)
         self.coefficients.initializeElementQuadrature(self.timeIntegration.t,self.q)
-        if self.stabilization != None:
+        if self.stabilization is not None:
             self.stabilization.initializeElementQuadrature(self.mesh,self.timeIntegration.t,self.q)
             self.stabilization.initializeTimeIntegration(self.timeIntegration)
-        if self.shockCapturing != None:
+        if self.shockCapturing is not None:
             self.shockCapturing.initializeElementQuadrature(self.mesh,self.timeIntegration.t,self.q)
     def calculateElementBoundaryQuadrature(self):
         pass
diff --git a/proteus/mprans/cSW2DCV.pyx b/proteus/mprans/cSW2DCV.pyx
index 126f80bea9..0c3a5e7243 100644
--- a/proteus/mprans/cSW2DCV.pyx
+++ b/proteus/mprans/cSW2DCV.pyx
@@ -10,9 +10,9 @@ cdef extern from "SW2DCV.h" namespace "proteus":
 	             int NNZ,
 		     int numDOFs,
 		     double* lumped_mass_matrix, 
-		     double* h_lstage,
-		     double* hu_lstage,
-		     double* hv_lstage, 
+		     double* h_old,
+		     double* hu_old,
+		     double* hv_old, 
 		     double* b_dof,
 		     double* high_order_hnp1, 
 		     double* high_order_hunp1, 
@@ -35,9 +35,9 @@ cdef extern from "SW2DCV.h" namespace "proteus":
 		     double g,
 		     int numDOFsPerEqn,
 		     double* lumped_mass_matrix, 
-		     double* h_lstage,
-		     double* hu_lstage,
-		     double* hv_lstage, 
+		     double* h_old,
+		     double* hu_old,
+		     double* hv_old, 
 		     double* b_dof,
 		     int* csrRowIndeces_DofLoops, 
 		     int* csrColumnOffsets_DofLoops, 
@@ -90,12 +90,6 @@ cdef extern from "SW2DCV.h" namespace "proteus":
                                double* h_dof_old, 
                                double* hu_dof_old, 
                                double* hv_dof_old,
-                               double* h_dof_old_old, 
-                               double* hu_dof_old_old, 
-                               double* hv_dof_old_old,
-                               double* h_dof_lstage, 
-                               double* hu_dof_lstage, 
-                               double* hv_dof_lstage,
                                double* b_dof,
                                double* h_dof, 
                                double* hu_dof, 
@@ -180,7 +174,8 @@ cdef extern from "SW2DCV.h" namespace "proteus":
 			       double* muH_minus_muL,
 			       double cE, 
 			       int LUMPED_MASS_MATRIX, 
-			       double dt, 
+			       double dt,
+			       int LINEAR_FRICTION,
 			       double mannings,
 			       double* quantDOFs, 			       
 			       int SECOND_CALL_CALCULATE_RESIDUAL,
@@ -229,12 +224,6 @@ cdef extern from "SW2DCV.h" namespace "proteus":
                                double* h_dof_old, 
                                double* hu_dof_old, 
                                double* hv_dof_old,
-                               double* h_dof_old_old, 
-                               double* hu_dof_old_old, 
-                               double* hv_dof_old_old,
-                               double* h_dof_lstage, 
-                               double* hu_dof_lstage, 
-                               double* hv_dof_lstage,
                                double* b_dof,
                                double* h_dof, 
                                double* hu_dof, 
@@ -319,7 +308,8 @@ cdef extern from "SW2DCV.h" namespace "proteus":
 			       double* muH_minus_muL,
 			       double cE,
 			       int LUMPED_MASS_MATRIX, 
-			       double dt, 
+			       double dt,
+			       int LINEAR_FRICTION,
 			       double mannings,						
 			       double* quantDOFs, 			       
 			       int SECOND_CALL_CALCULATE_RESIDUAL,
@@ -690,9 +680,9 @@ cdef class cSW2DCV_base:
 	       int NNZ,
 	       int numDOFs,
 	       numpy.ndarray lumped_mass_matrix, 
-	       numpy.ndarray h_lstage, 
-	       numpy.ndarray hu_lstage, 
-	       numpy.ndarray hv_lstage, 
+	       numpy.ndarray h_old, 
+	       numpy.ndarray hu_old, 
+	       numpy.ndarray hv_old, 
 	       numpy.ndarray b_dof, 
 	       numpy.ndarray high_order_hnp1,
 	       numpy.ndarray high_order_hunp1,
@@ -715,9 +705,9 @@ cdef class cSW2DCV_base:
        			    NNZ,
 	                    numDOFs,
 			    <double*> lumped_mass_matrix.data, 
-			    <double*> h_lstage.data, 
-			    <double*> hu_lstage.data, 
-			    <double*> hv_lstage.data, 
+			    <double*> h_old.data, 
+			    <double*> hu_old.data, 
+			    <double*> hv_old.data, 
 			    <double*> b_dof.data, 
 			    <double*> high_order_hnp1.data,
 			    <double*> high_order_hunp1.data,
@@ -740,9 +730,9 @@ cdef class cSW2DCV_base:
        	       double g,
 	       int numDOFsPerEqn,
 	       numpy.ndarray lumped_mass_matrix, 
-	       numpy.ndarray h_lstage, 
-	       numpy.ndarray hu_lstage, 
-	       numpy.ndarray hv_lstage, 
+	       numpy.ndarray h_old, 
+	       numpy.ndarray hu_old, 
+	       numpy.ndarray hv_old, 
 	       numpy.ndarray b_dof, 
 	       numpy.ndarray csrRowIndeces_DofLoops, 
 	       numpy.ndarray csrColumnOffsets_DofLoops, 
@@ -759,9 +749,9 @@ cdef class cSW2DCV_base:
 			    g, 
 			    numDOFsPerEqn,
 			    <double*> lumped_mass_matrix.data, 
-			    <double*> h_lstage.data, 
-			    <double*> hu_lstage.data, 
-			    <double*> hv_lstage.data, 
+			    <double*> h_old.data, 
+			    <double*> hu_old.data, 
+			    <double*> hv_old.data, 
 			    <double*> b_dof.data,   
 			    <int*> csrRowIndeces_DofLoops.data,
 			    <int*> csrColumnOffsets_DofLoops.data,
@@ -815,12 +805,6 @@ cdef class cSW2DCV_base:
                          numpy.ndarray h_dof_old, 
                          numpy.ndarray hu_dof_old, 
                          numpy.ndarray hv_dof_old, 
-                         numpy.ndarray h_dof_old_old, 
-                         numpy.ndarray hu_dof_old_old, 
-                         numpy.ndarray hv_dof_old_old, 
-                         numpy.ndarray h_dof_lstage, 
-                         numpy.ndarray hu_dof_lstage, 
-                         numpy.ndarray hv_dof_lstage, 
                          numpy.ndarray b_dof, 
                          numpy.ndarray h_dof, 
                          numpy.ndarray hu_dof, 
@@ -905,7 +889,8 @@ cdef class cSW2DCV_base:
 			 numpy.ndarray muH_minus_muL,
 			 double cE, 
 			 int LUMPED_MASS_MATRIX,	
-			 double dt, 
+			 double dt,
+			 int LINEAR_FRICTION,
 			 double mannings,
 			 numpy.ndarray quantDOFs,	
 			 int SECOND_CALL_CALCULATE_RESIDUAL, 
@@ -957,12 +942,6 @@ cdef class cSW2DCV_base:
             <double*> h_dof_old.data, 
             <double*> hu_dof_old.data, 
             <double*> hv_dof_old.data, 
-            <double*> h_dof_old_old.data, 
-            <double*> hu_dof_old_old.data, 
-            <double*> hv_dof_old_old.data, 
-            <double*> h_dof_lstage.data, 
-            <double*> hu_dof_lstage.data, 
-            <double*> hv_dof_lstage.data, 
             <double*> b_dof.data, 
             <double*> h_dof.data, 
             <double*> hu_dof.data, 
@@ -1047,7 +1026,8 @@ cdef class cSW2DCV_base:
 	    <double*> muH_minus_muL.data,
 	    cE, 
 	    LUMPED_MASS_MATRIX, 
-	    dt, 
+	    dt,
+	    LINEAR_FRICTION,
 	    mannings,
 	    <double*> quantDOFs.data,
 	    SECOND_CALL_CALCULATE_RESIDUAL,
@@ -1097,12 +1077,6 @@ cdef class cSW2DCV_base:
                          numpy.ndarray h_dof_old, 
                          numpy.ndarray hu_dof_old, 
                          numpy.ndarray hv_dof_old, 
-                         numpy.ndarray h_dof_old_old, 
-                         numpy.ndarray hu_dof_old_old, 
-                         numpy.ndarray hv_dof_old_old, 
-                         numpy.ndarray h_dof_lstage, 
-                         numpy.ndarray hu_dof_lstage, 
-                         numpy.ndarray hv_dof_lstage, 
                          numpy.ndarray b_dof, 
                          numpy.ndarray h_dof, 
                          numpy.ndarray hu_dof, 
@@ -1187,7 +1161,8 @@ cdef class cSW2DCV_base:
 			 numpy.ndarray muH_minus_muL,
 			 double cE,
 			 int LUMPED_MASS_MATRIX, 
-			 double dt, 
+			 double dt,
+			 int LINEAR_FRICTION,
 			 double mannings,
 			 numpy.ndarray quantDOFs, 
 			 int SECOND_CALL_CALCULATE_RESIDUAL,
@@ -1239,12 +1214,6 @@ cdef class cSW2DCV_base:
             <double*> h_dof_old.data, 
             <double*> hu_dof_old.data, 
             <double*> hv_dof_old.data, 
-            <double*> h_dof_old_old.data, 
-            <double*> hu_dof_old_old.data, 
-            <double*> hv_dof_old_old.data, 
-            <double*> h_dof_lstage.data, 
-            <double*> hu_dof_lstage.data, 
-            <double*> hv_dof_lstage.data, 
             <double*> b_dof.data, 
             <double*> h_dof.data, 
             <double*> hu_dof.data, 
@@ -1329,7 +1298,8 @@ cdef class cSW2DCV_base:
 	    <double*> muH_minus_muL.data,
 	    cE,
 	    LUMPED_MASS_MATRIX, 
-	    dt, 
+	    dt,
+	    LINEAR_FRICTION,
 	    mannings,
 	    <double*> quantDOFs.data,
 	    SECOND_CALL_CALCULATE_RESIDUAL,
diff --git a/proteus/tests/SWEs/1D_dambreak_flat_bottom/L2_batch.py b/proteus/tests/SWEs/1D_dambreak_flat_bottom/L2_batch.py
index c67e46d04b..32f713e750 100644
--- a/proteus/tests/SWEs/1D_dambreak_flat_bottom/L2_batch.py
+++ b/proteus/tests/SWEs/1D_dambreak_flat_bottom/L2_batch.py
@@ -11,9 +11,10 @@
 
 simFlagsList[0]['errorQuantities']=['u']
 simFlagsList[0]['errorTypes']= ['numericalSolution'] #compute error in soln and glob. mass bal
-simFlagsList[0]['errorNorms']= ['L1','L2'] #compute L2 norm in space or H0 or ...
+simFlagsList[0]['errorNorms']= ['L1'] #compute L2 norm in space or H0 or ...
 simFlagsList[0]['errorTimes']= ['Last'] #'All', 'Last'
 simFlagsList[0]['echo']=True
+simFlagsList[0]['echoRelativeErrors']=True
 
 #
 start
diff --git a/proteus/tests/SWEs/1D_dambreak_flat_bottom/sw_hump_2d_c0p1_n.py b/proteus/tests/SWEs/1D_dambreak_flat_bottom/sw_hump_2d_c0p1_n.py
index db82f0cf49..9b350e657c 100644
--- a/proteus/tests/SWEs/1D_dambreak_flat_bottom/sw_hump_2d_c0p1_n.py
+++ b/proteus/tests/SWEs/1D_dambreak_flat_bottom/sw_hump_2d_c0p1_n.py
@@ -4,8 +4,8 @@
 
 refinement=3
 runCFL=0.5
-#timeIntegration_sw2d = "SSP33"
-timeIntegration_sw2d = "FE"
+timeIntegration_sw2d = "SSP33"
+#timeIntegration_sw2d = "FE"
 
 multilevelNonlinearSolver  = Newton
 if (LUMPED_MASS_MATRIX==1):
diff --git a/proteus/tests/SWEs/dam_over_bumps/sw_hump_2d_c0p1_n.py b/proteus/tests/SWEs/dam_over_bumps/sw_hump_2d_c0p1_n.py
index eff1d47a90..c070e5387c 100644
--- a/proteus/tests/SWEs/dam_over_bumps/sw_hump_2d_c0p1_n.py
+++ b/proteus/tests/SWEs/dam_over_bumps/sw_hump_2d_c0p1_n.py
@@ -8,8 +8,7 @@
 reflecting_BCs=ct.opts.reflecting_BCs
 refinement=ct.opts.refinement
 runCFL=0.5
-timeIntegration_sw2d = "SSP33"
-#timeIntegration_sw2d = "FE"
+sspOrder = 3
 
 multilevelNonlinearSolver  = Newton
 if (ct.LUMPED_MASS_MATRIX==1):
@@ -20,12 +19,8 @@
 
 timeIntegration = ct.SW2DCV.RKEV 
 stepController = Min_dt_controller
-if timeIntegration_sw2d == "SSP33": #mwf hack
-    timeOrder = 3
-    nStagesTime = 3
-else:
-    timeOrder = 1
-    nStagesTime = 1
+timeOrder = sspOrder
+nStagesTime = sspOrder
 
 rtol_u[0] = 1.0e-4
 rtol_u[1] = 1.0e-4
diff --git a/proteus/tests/SWEs/dam_over_bumps/sw_hump_2d_p.py b/proteus/tests/SWEs/dam_over_bumps/sw_hump_2d_p.py
index e37a16bf8d..747c7056d2 100644
--- a/proteus/tests/SWEs/dam_over_bumps/sw_hump_2d_p.py
+++ b/proteus/tests/SWEs/dam_over_bumps/sw_hump_2d_p.py
@@ -33,7 +33,6 @@
 #This is relevant just when use_second_order_NonFlatB_with_EV_stabilization=True
 cE=1
 LUMPED_MASS_MATRIX=0
-#reflecting_BCs=1
 
 bt = domain.boundaryTags
 bt['front'] = bt['bottom']
diff --git a/proteus/tests/SWEs/malpasset/sw_hump_2d_c0p1_n.py b/proteus/tests/SWEs/malpasset/sw_hump_2d_c0p1_n.py
index 3ce955f978..4755e109d3 100644
--- a/proteus/tests/SWEs/malpasset/sw_hump_2d_c0p1_n.py
+++ b/proteus/tests/SWEs/malpasset/sw_hump_2d_c0p1_n.py
@@ -9,7 +9,7 @@
 
 reflecting_BCs=ct.opts.reflecting_BCs
 refinement=ct.opts.refinement
-runCFL=0.5
+runCFL=0.33
 #timeIntegration_sw2d = "SSP33"
 timeIntegration_sw2d = "FE"
 
@@ -19,7 +19,6 @@
 else:
     levelNonlinearSolver = ExplicitConsistentMassMatrixShallowWaterEquationsSolver
 
-
 timeIntegration = ct.SW2DCV.RKEV 
 stepController = Min_dt_controller
 if timeIntegration_sw2d == "SSP33": #mwf hack
@@ -81,7 +80,7 @@
 levelLinearSolver = LU
 
 #conservativeFlux = {0:'pwl'}
-tnList = [0]+[float(n)*ct.T/float(ct.nDTout) for n in range(1,ct.nDTout+1)]
+tnList = [0.,1E-6]+[float(n)*ct.T/float(ct.nDTout) for n in range(1,ct.nDTout+1)]
 #archiveFlag = ArchiveFlags.EVERY_SEQUENCE_STEP
 
 p = PointGauges(gauges=(( ('h'), ((5550.0,4400.0, 0),
diff --git a/proteus/tests/SWEs/malpasset/sw_hump_2d_p.py b/proteus/tests/SWEs/malpasset/sw_hump_2d_p.py
index 4587965c4c..e4794da49b 100644
--- a/proteus/tests/SWEs/malpasset/sw_hump_2d_p.py
+++ b/proteus/tests/SWEs/malpasset/sw_hump_2d_p.py
@@ -21,8 +21,8 @@
 g = 9.81
 # PARAMETERS #
 mannings=0.033 #0.025 #0.025#m^{-1/3} after Ying etal 2009
-cE=5
-LUMPED_MASS_MATRIX=1
+cE=1
+LUMPED_MASS_MATRIX=0
 
 domain = None#RectangularDomain(L=L)
 meshfile = "mal_50sec"
diff --git a/proteus/tests/SWEs/parabolid_with_friction/1D/small_scales/L2_batch.py b/proteus/tests/SWEs/parabolid_with_friction/1D/L2_batch.py
similarity index 76%
rename from proteus/tests/SWEs/parabolid_with_friction/1D/small_scales/L2_batch.py
rename to proteus/tests/SWEs/parabolid_with_friction/1D/L2_batch.py
index ec3463cdd7..144ae5b8ca 100644
--- a/proteus/tests/SWEs/parabolid_with_friction/1D/small_scales/L2_batch.py
+++ b/proteus/tests/SWEs/parabolid_with_friction/1D/L2_batch.py
@@ -4,7 +4,6 @@
 comm = Comm.get()
 ctx = Context.get()
 
-
 # simulation flags for error analysis
 #
 # simFlagsList is initialized in proteus.iproteus
@@ -12,9 +11,10 @@
 
 simFlagsList[0]['errorQuantities']=['u']
 simFlagsList[0]['errorTypes']= ['numericalSolution'] #compute error in soln and glob. mass bal
-simFlagsList[0]['errorNorms']= ['L1','L2'] #compute L2 norm in space or H0 or ...
+simFlagsList[0]['errorNorms']= ['L1','L2','LI'] #compute L2 norm in space or H0 or ...
 simFlagsList[0]['errorTimes']= ['Last'] #'All', 'Last'
 simFlagsList[0]['echo']=True
+simFlagsList[0]['echoRelativeErrors']=True
 
 #
 start
diff --git a/proteus/tests/SWEs/parabolid_with_friction/1D/small_scales/clean.sh b/proteus/tests/SWEs/parabolid_with_friction/1D/clean.sh
similarity index 100%
rename from proteus/tests/SWEs/parabolid_with_friction/1D/small_scales/clean.sh
rename to proteus/tests/SWEs/parabolid_with_friction/1D/clean.sh
diff --git a/proteus/tests/SWEs/parabolid_with_friction/1D/small_scales/sw_hump_2d_c0p1_n.py b/proteus/tests/SWEs/parabolid_with_friction/1D/sw_hump_2d_c0p1_n.py
similarity index 85%
rename from proteus/tests/SWEs/parabolid_with_friction/1D/small_scales/sw_hump_2d_c0p1_n.py
rename to proteus/tests/SWEs/parabolid_with_friction/1D/sw_hump_2d_c0p1_n.py
index 743ff40ae8..837e10c009 100644
--- a/proteus/tests/SWEs/parabolid_with_friction/1D/small_scales/sw_hump_2d_c0p1_n.py
+++ b/proteus/tests/SWEs/parabolid_with_friction/1D/sw_hump_2d_c0p1_n.py
@@ -3,9 +3,8 @@
 from sw_hump_2d_p import *
 
 refinement=3
-runCFL=0.5
-timeIntegration_sw2d = "SSP33"
-#timeIntegration_sw2d = "FE"
+runCFL=0.4
+SSPOrder=3 
 
 multilevelNonlinearSolver  = Newton
 if (LUMPED_MASS_MATRIX==1):
@@ -15,12 +14,8 @@
 
 timeIntegration = SW2DCV.RKEV 
 stepController = Min_dt_controller
-if timeIntegration_sw2d == "SSP33": #mwf hack
-    timeOrder = 3
-    nStagesTime = 3
-else:
-    timeOrder = 1
-    nStagesTime = 1
+timeOrder = SSPOrder
+nStagesTime = SSPOrder
 
 rtol_u[0] = 1.0e-4
 rtol_u[1] = 1.0e-4
@@ -74,4 +69,4 @@
 levelLinearSolver = LU
 
 #conservativeFlux = {0:'pwl'}
-#tnList=[0.,1E-6]+[float(n)*T/float(nDTout) for n in range(1,nDTout+1)]
+tnList=[0.,1E-6]+[float(n)*T/float(nDTout) for n in range(1,nDTout+1)]
diff --git a/proteus/tests/SWEs/parabolid_with_friction/1D/small_scales/sw_hump_2d_p.py b/proteus/tests/SWEs/parabolid_with_friction/1D/sw_hump_2d_p.py
similarity index 95%
rename from proteus/tests/SWEs/parabolid_with_friction/1D/small_scales/sw_hump_2d_p.py
rename to proteus/tests/SWEs/parabolid_with_friction/1D/sw_hump_2d_p.py
index 86377a0774..e48136d1ba 100644
--- a/proteus/tests/SWEs/parabolid_with_friction/1D/small_scales/sw_hump_2d_p.py
+++ b/proteus/tests/SWEs/parabolid_with_friction/1D/sw_hump_2d_p.py
@@ -7,24 +7,24 @@
 
 nd=2
 
-T=1
+T=1000.00
 nDTout=10
 
-L=(10.0,1.0)
+L=(8000.0,800.0)
 g = 9.81
 domain = RectangularDomain(L=L,x=[0,0,0])
 
-h0=1
-a=3
-B=1
-k=0.5
+h0=10
+a=3000
+B=2
+k=0.001
 p = np.sqrt(8*g*h0)/a
 s = np.sqrt(p**2 - k**2)/2.
 mannings=k
 
-#This is relevant just when use_second_order_NonFlatB_with_EV_stabilization=True
 cE=1
 LUMPED_MASS_MATRIX=0
+LINEAR_FRICTION=1
 
 bt = domain.boundaryTags
 bt['front'] = bt['bottom']
@@ -161,5 +161,5 @@ def getDFBC_v(x,flag):
 #########################################
 bathymetry={0:bathymetry_function}
 LevelModelType = SW2DCV.LevelModel
-coefficients = SW2DCV.Coefficients(g=g,bathymetry=bathymetry,cE=cE,LUMPED_MASS_MATRIX=LUMPED_MASS_MATRIX,mannings=mannings)
+coefficients = SW2DCV.Coefficients(g=g,bathymetry=bathymetry,cE=cE,LUMPED_MASS_MATRIX=LUMPED_MASS_MATRIX,LINEAR_FRICTION=LINEAR_FRICTION,mannings=mannings)
 
diff --git a/proteus/tests/SWEs/parabolid_with_friction/2D/L2_batch.py b/proteus/tests/SWEs/parabolid_with_friction/2D/L2_batch.py
index 857a1e5183..e469c46ba6 100644
--- a/proteus/tests/SWEs/parabolid_with_friction/2D/L2_batch.py
+++ b/proteus/tests/SWEs/parabolid_with_friction/2D/L2_batch.py
@@ -12,10 +12,10 @@
 
 simFlagsList[0]['errorQuantities']=['u']
 simFlagsList[0]['errorTypes']= ['numericalSolution'] #compute error in soln and glob. mass bal
-simFlagsList[0]['errorNorms']= ['L1'] #compute L2 norm in space or H0 or ...
+simFlagsList[0]['errorNorms']= ['L1','LI'] #compute L2 norm in space or H0 or ...
 simFlagsList[0]['errorTimes']= ['Last'] #'All', 'Last'
 simFlagsList[0]['echo']=True
-
+simFlagsList[0]['echoRelativeErrors']=True
 #
 start
 quit
diff --git a/proteus/tests/SWEs/parabolid_with_friction/2D/sw_hump_2d_c0p1_n.py b/proteus/tests/SWEs/parabolid_with_friction/2D/sw_hump_2d_c0p1_n.py
index 74dd27906a..92489e9884 100644
--- a/proteus/tests/SWEs/parabolid_with_friction/2D/sw_hump_2d_c0p1_n.py
+++ b/proteus/tests/SWEs/parabolid_with_friction/2D/sw_hump_2d_c0p1_n.py
@@ -2,8 +2,8 @@
 from proteus.default_n import *
 from sw_hump_2d_p import *
 
-refinement=6
-runCFL=0.25
+refinement=2
+runCFL=0.5
 timeIntegration_sw2d = "SSP33"
 #timeIntegration_sw2d = "FE"
 
diff --git a/proteus/tests/SWEs/parabolid_with_friction/2D/sw_hump_2d_p.py b/proteus/tests/SWEs/parabolid_with_friction/2D/sw_hump_2d_p.py
index 44afe7b444..13a5a3134f 100644
--- a/proteus/tests/SWEs/parabolid_with_friction/2D/sw_hump_2d_p.py
+++ b/proteus/tests/SWEs/parabolid_with_friction/2D/sw_hump_2d_p.py
@@ -8,16 +8,16 @@
 
 nd=2
 
-T=1. #2*math.pi/omega
+T=1000. #2*math.pi/omega
 nDTout=10
 
-L=(10.0,10.0)
+L=(10000.0,10000.0)
 g = 9.81
 # PARAMETERS #
-h0=1.0
-a=3.
-B=2.
-k=0.5
+h0=10.0
+a=3000.
+B=5.
+k=0.002
 p = np.sqrt(8*g*h0)/a
 s = np.sqrt(p**2 - k**2)/2.
 mannings=k
@@ -27,6 +27,7 @@
 #This is relevant just when use_second_order_NonFlatB_with_EV_stabilization=True
 cE=1
 LUMPED_MASS_MATRIX=0
+LINEAR_FRICTION=1
 
 bt = domain.boundaryTags
 bt['front'] = bt['bottom']
@@ -145,5 +146,5 @@ def getDFBC_v(x,flag):
 #########################################
 bathymetry={0:bathymetry_function}
 LevelModelType = SW2DCV.LevelModel
-coefficients = SW2DCV.Coefficients(g=g,bathymetry=bathymetry,cE=cE,LUMPED_MASS_MATRIX=LUMPED_MASS_MATRIX,mannings=mannings)
+coefficients = SW2DCV.Coefficients(g=g,bathymetry=bathymetry,cE=cE,LUMPED_MASS_MATRIX=LUMPED_MASS_MATRIX,LINEAR_FRICTION=LINEAR_FRICTION,mannings=mannings)
 
diff --git a/setup.py b/setup.py
index 6a396e5abb..7875d0b6f2 100644
--- a/setup.py
+++ b/setup.py
@@ -39,7 +39,7 @@
         break
 
 setup(name='proteus',
-      version='1.4.1.dev0',
+      version='1.4.2',
       classifiers=[
           'Development Status :: 4 - Beta',
           'Environment :: Console',
@@ -416,9 +416,13 @@
                                   libraries=['hdf5','stdc++','m',PROTEUS_DAETK_LIB]+PROTEUS_PETSC_LIBS+PROTEUS_MPI_LIBS+PROTEUS_HDF5_LIBS,
                                   extra_link_args=PROTEUS_EXTRA_LINK_ARGS + PROTEUS_PETSC_EXTRA_LINK_ARGS,
                                   extra_compile_args=PROTEUS_EXTRA_COMPILE_ARGS + PROTEUS_PETSC_EXTRA_COMPILE_ARGS),
-                    Extension("mprans.cNCLS",["proteus/mprans/cNCLS.pyx"],depends=["proteus/mprans/NCLS.h"], language="c++",
+                    Extension("mprans.cNCLS",["proteus/mprans/cNCLS.pyx"],
+                              depends=["proteus/mprans/NCLS.h"] + ["proteus/ModelFactory.h","proteus/CompKernel.h"],
+                              language="c++",
                               include_dirs=[numpy.get_include(), 'proteus']),
-                   Extension("mprans.cMCorr",["proteus/mprans/cMCorr.pyx"],depends=["proteus/mprans/MCorr.h"], define_macros=[('PROTEUS_LAPACK_H',PROTEUS_LAPACK_H),
+                   Extension("mprans.cMCorr",["proteus/mprans/cMCorr.pyx"],
+                             depends=["proteus/mprans/MCorr.h"] + ["proteus/ModelFactory.h","proteus/CompKernel.h"],
+                             define_macros=[('PROTEUS_LAPACK_H',PROTEUS_LAPACK_H),
                                             ('PROTEUS_LAPACK_INTEGER',PROTEUS_LAPACK_INTEGER),
                                             ('PROTEUS_BLAS_H',PROTEUS_BLAS_H)],language="c++",
                              include_dirs=[numpy.get_include(), 'proteus'],
@@ -427,7 +431,8 @@
                              libraries=['m',PROTEUS_LAPACK_LIB,PROTEUS_BLAS_LIB],
                              extra_compile_args=PROTEUS_EXTRA_COMPILE_ARGS,
                              extra_link_args=PROTEUS_EXTRA_LINK_ARGS),
-                   Extension("mprans.cRANS2P",["proteus/mprans/cRANS2P.pyx"], depends=["proteus/mprans/RANS2P.h"] + ["proteus/ModelFactory.h","proteus/CompKernel.h"],
+                   Extension("mprans.cRANS2P",["proteus/mprans/cRANS2P.pyx"],
+                             depends=["proteus/mprans/RANS2P.h"] + ["proteus/ModelFactory.h","proteus/CompKernel.h"],
                              language="c++", include_dirs=[numpy.get_include(), 'proteus']),
                    Extension("mprans.cRANS2P2D",["proteus/mprans/cRANS2P2D.pyx"],
                              depends=["proteus/mprans/RANS2P2D.h"] + ["proteus/ModelFactory.h","proteus/CompKernel.h"],
@@ -497,6 +502,10 @@
                   (os.path.join(proteus_install_path,'tests','griffiths_lane_6'),
                    ['proteus/tests/griffiths_lane_6/richards_expected.h5',
                     'proteus/tests/griffiths_lane_6/elastoplastic_expected.h5']),
+                  (os.path.join(proteus_install_path,'tests','levelset'),
+                   ['proteus/tests/levelset/rotation/rotation_c0p1cg_vbdf_2_level_1_expected.h5',
+                    'proteus/tests/levelset/vortex2D/vortex_c0p1cg_vbdf_2_level_1_expected.h5',
+                    'proteus/tests/levelset/vortex/vortex_c0p1cg_bdf_2_level_1_expected.h5']),
                   (os.path.join(proteus_install_path,'tests','LS_with_edgeBased_EV','VOF','comparison_files'),
                    ['proteus/tests/LS_with_edgeBased_EV/VOF/comparison_files/vof_level_3_SmoothnessBased.h5',
                     'proteus/tests/LS_with_edgeBased_EV/VOF/comparison_files/vof_level_3_EV1.h5',