diff --git a/proteus/m_comp_co2/m_comp_co2.py b/proteus/m_comp_co2/m_comp_co2.py
index 180fd88998..4b8862ddfb 100644
--- a/proteus/m_comp_co2/m_comp_co2.py
+++ b/proteus/m_comp_co2/m_comp_co2.py
@@ -1,6 +1,5 @@
 from __future__ import division
 from builtins import range
-#from past.utils import old_div
 import os
 import proteus
 from .cm_comp_co2 import *
diff --git a/proteus/richards/ADR.py b/proteus/richards/ADR.py
index 94d582c2ea..24163e5231 100644
--- a/proteus/richards/ADR.py
+++ b/proteus/richards/ADR.py
@@ -1,6 +1,5 @@
 from __future__ import division
 from builtins import range
-from past.utils import old_div
 import proteus
 from .cADR import *
 import numpy as np
@@ -58,7 +57,7 @@ def choose_dt(self):
         comm = Comm.get()
         maxCFL = 1.0e-6
         maxCFL = max(maxCFL, comm.globalMax(self.cfl.max()))
-        self.dt = old_div(self.runCFL, maxCFL)
+        self.dt = self.runCFL / maxCFL
         if self.dtLast is None:
             self.dtLast = self.dt
         self.t = self.tLast + self.dt
@@ -98,7 +97,7 @@ def updateStage(self):
                 logEvent("Second stage of SSP33 method", level=4)
                 for ci in range(self.nc):
                     self.u_dof_stage[ci][self.lstage][:] = self.transport.u[ci].dof
-                    self.u_dof_stage[ci][self.lstage] *= old_div(1., 4.)
+                    self.u_dof_stage[ci][self.lstage] *= 1. / 4.
                     self.u_dof_stage[ci][self.lstage] += 3. / 4. * self.u_dof_last[ci]
                     # Update u_dof_old
                     self.transport.u_dof_old[:] = self.u_dof_stage[ci][self.lstage]
@@ -106,7 +105,7 @@ def updateStage(self):
                 logEvent("Third stage of SSP33 method", level=4)
                 for ci in range(self.nc):
                     self.u_dof_stage[ci][self.lstage][:] = self.transport.u[ci].dof
-                    self.u_dof_stage[ci][self.lstage] *= old_div(2.0, 3.0)
+                    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]
                     # update u_dof_old
                     self.transport.u_dof_old[:] = self.u_dof_last[ci]
@@ -123,7 +122,7 @@ def updateStage(self):
                 logEvent("Second stage of SSP22 method", level=4)
                 for ci in range(self.nc):
                     self.u_dof_stage[ci][self.lstage][:] = self.transport.u[ci].dof
-                    self.u_dof_stage[ci][self.lstage][:] *= old_div(1., 2.)
+                    self.u_dof_stage[ci][self.lstage][:] *= 1. / 2.
                     self.u_dof_stage[ci][self.lstage][:] += 1. / 2. * self.u_dof_last[ci]
                     # update u_dof_old
                     self.transport.u_dof_old[:] = self.u_dof_last[ci]
@@ -600,7 +599,7 @@ def __init__(self,
         self.edge_based_cfl = np.zeros(self.u[0].dof.shape)+100
         #mesh
         #self.q['x'] = np.zeros((self.mesh.nElements_global,self.nQuadraturePoints_element,3),'d')
-        self.q[('dV_u', 0)] = (old_div(1.0, self.mesh.nElements_global)) * np.ones((self.mesh.nElements_global, self.nQuadraturePoints_element), 'd')
+        self.q[('dV_u', 0)] = (1.0 / self.mesh.nElements_global) * np.ones((self.mesh.nElements_global, self.nQuadraturePoints_element), 'd')
         self.ebqe['x'] = np.zeros((self.mesh.nExteriorElementBoundaries_global,self.nElementBoundaryQuadraturePoints_elementBoundary,3),'d')
         self.q[('u',0)] = np.zeros((self.mesh.nElements_global,self.nQuadraturePoints_element),'d')
         self.q[('grad(u)',0)] = np.zeros((self.mesh.nElements_global,self.nQuadraturePoints_element,self.nSpace_global),'d')
@@ -723,14 +722,14 @@ def __init__(self,
         if 'penalty' in self.ebq_global:
             for ebN in range(self.mesh.nElementBoundaries_global):
                 for k in range(self.nElementBoundaryQuadraturePoints_elementBoundary):
-                    self.ebq_global['penalty'][ebN,k] = old_div(self.numericalFlux.penalty_constant,(self.mesh.elementBoundaryDiametersArray[ebN]**self.numericalFlux.penalty_power))
+                    self.ebq_global['penalty'][ebN,k] = self.numericalFlux.penalty_constant / (self.mesh.elementBoundaryDiametersArray[ebN]**self.numericalFlux.penalty_power)
         #penalty term
         #cek move  to Numerical flux initialization
         if 'penalty' in self.ebqe:
             for ebNE in range(self.mesh.nExteriorElementBoundaries_global):
                 ebN = self.mesh.exteriorElementBoundariesArray[ebNE]
                 for k in range(self.nElementBoundaryQuadraturePoints_elementBoundary):
-                    self.ebqe['penalty'][ebNE,k] = old_div(self.numericalFlux.penalty_constant,self.mesh.elementBoundaryDiametersArray[ebN]**self.numericalFlux.penalty_power)
+                    self.ebqe['penalty'][ebNE,k] = self.numericalFlux.penalty_constant / self.mesh.elementBoundaryDiametersArray[ebN]**self.numericalFlux.penalty_power
         logEvent(memory("numericalFlux","OneLevelTransport"),level=4)
         self.elementEffectiveDiametersArray  = self.mesh.elementInnerDiametersArray
         #use post processing tools to get conservative fluxes, None by default
diff --git a/proteus/richards/Richards.py b/proteus/richards/Richards.py
index ad01ddcc2b..a199260d29 100644
--- a/proteus/richards/Richards.py
+++ b/proteus/richards/Richards.py
@@ -1,6 +1,5 @@
 from __future__ import division
 from builtins import range
-#from past.utils import old_div
 import proteus
 from .cRichards import *
 import numpy as np
diff --git a/scripts/marinExtractHeight.py b/scripts/marinExtractHeight.py
index 418ae391fe..cf214f9642 100755
--- a/scripts/marinExtractHeight.py
+++ b/scripts/marinExtractHeight.py
@@ -68,7 +68,7 @@
             phi = pdata.GetArray("phid").GetTuple1(i)
 
             if (phi > 0.0) and (phi_old < 0.0):
-                height = old_div((float(i-1) + (old_div(phi_old,(phi_old-phi)))),float(opts.resolution))
+                height = (float(i-1) + phi_old / (phi_old - phi)) / float(opts.resolution)
             phi_old=phi
 
         if (height > 0.0):