diff --git a/.travis.yml b/.travis.yml
index e03daae4..649df162 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -1,6 +1,7 @@
 language: python
 python:
   - 2.7
+  - 3.5
 before_install:
   - sudo apt-get update -qq
   - sudo apt-get install -qq gfortran liblapack-pic
diff --git a/rptest/1d/test_rp1_meqn2.py b/rptest/1d/test_rp1_meqn2.py
index e121ac79..2a100dea 100644
--- a/rptest/1d/test_rp1_meqn2.py
+++ b/rptest/1d/test_rp1_meqn2.py
@@ -1,4 +1,6 @@
 
+from __future__ import absolute_import
+from __future__ import print_function
 import test_rp1_meqn2
 import numpy
 
@@ -13,6 +15,6 @@
 qr = [1.,0.]
 wave,s,amdq,apdq = test_rp1_meqn2.rp1_driver.call_rp1(2,ql,qr)
 
-print "s = ",s
-print "wave = "
-print wave
+print("s = ",s)
+print("wave = ")
+print(wave)
diff --git a/src/__init__.py b/src/__init__.py
index 8507663c..e3280710 100755
--- a/src/__init__.py
+++ b/src/__init__.py
@@ -4,91 +4,93 @@
 Wave propagation Riemann solvers implemented in Python and Fortran.
 """
 
+from __future__ import absolute_import
+from __future__ import print_function
 rp_solver_list_1d = []
 rp_solver_list_2d = []
 rp_solver_list_3d = []
 
 # Import 1d Riemann solvers
-from advection_1D_py import advection_1D
-from vc_advection_1D_py import vc_advection_1D
-from acoustics_1D_py import acoustics_1D
-from burgers_1D_py import burgers_1D
-from shallow_1D_py import shallow_roe_1D, shallow_hll_1D, shallow_exact_1D
-from euler_1D_py import euler_roe_1D, euler_hll_1D, euler_hllc_1D, euler_exact_1D
-from nonlinear_elasticity_1D_py import nonlinear_elasticity_1D
-import static
+from . import advection_1D_py
+from . import vc_advection_1D_py
+from . import acoustics_1D_py
+from . import burgers_1D_py
+from . import shallow_1D_py
+from . import euler_1D_py
+from . import nonlinear_elasticity_1D_py
+from . import static
 
-import acoustics_1D_constants
-import acoustics_variable_1D_constants
-import advection_1D_constants
-import burgers_1D_constants
-import euler_with_efix_1D_constants
-import nonlinear_elasticity_fwave_1D_constants
-import reactive_euler_with_efix_1D_constants
-import shallow_roe_with_efix_1D_constants
-import traffic_1D_constants
-import traffic_vc_1D_constants
-import acoustics_2D_constants
-import acoustics_mapped_2D_constants
-import advection_2D_constants
-import burgers_2D_constants
-import euler_mapgrid_2D_constants
-import euler_5wave_2D_constants
-import euler_4wave_2D_constants
-import kpp_2D_constants
-import psystem_2D_constants
-import shallow_roe_with_efix_2D_constants
-import shallow_sphere_2D_constants
-import vc_acoustics_2D_constants
-import vc_advection_2D_constants
-import vc_elasticity_2D_constants
-import vc_acoustics_3D_constants
-import euler_3D_constants
-import burgers_3D_constants
-import vc_advection_3D_constants
+from . import acoustics_1D_constants
+from . import acoustics_variable_1D_constants
+from . import advection_1D_constants
+from . import burgers_1D_constants
+from . import euler_with_efix_1D_constants
+from . import nonlinear_elasticity_fwave_1D_constants
+from . import reactive_euler_with_efix_1D_constants
+from . import shallow_roe_with_efix_1D_constants
+from . import traffic_1D_constants
+from . import traffic_vc_1D_constants
+from . import acoustics_2D_constants
+from . import acoustics_mapped_2D_constants
+from . import advection_2D_constants
+from . import burgers_2D_constants
+from . import euler_mapgrid_2D_constants
+from . import euler_5wave_2D_constants
+from . import euler_4wave_2D_constants
+from . import kpp_2D_constants
+from . import psystem_2D_constants
+from . import shallow_roe_with_efix_2D_constants
+from . import shallow_sphere_2D_constants
+from . import vc_acoustics_2D_constants
+from . import vc_advection_2D_constants
+from . import vc_elasticity_2D_constants
+from . import vc_acoustics_3D_constants
+from . import euler_3D_constants
+from . import burgers_3D_constants
+from . import vc_advection_3D_constants
 
 
 try:
-    import acoustics_1D
-    import acoustics_variable_1D
-    import acoustics_1D_ptwise
-    import advection_1D
-    import advection_1D_ptwise
-    import burgers_1D
-    import euler_with_efix_1D
-    import nonlinear_elasticity_fwave_1D
-    import reactive_euler_with_efix_1D
-    import shallow_roe_with_efix_1D
-    import traffic_1D
-    import traffic_vc_1D
-    import acoustics_2D
-    import acoustics_mapped_2D
-    import acoustics_2D_ptwise
-    import advection_2D
-    import burgers_2D
-    import euler_mapgrid_2D
-    import euler_5wave_2D
-    import euler_4wave_2D
-    import kpp_2D
-    import psystem_2D
-    import shallow_roe_with_efix_2D
-    import shallow_sphere_2D
-    import vc_acoustics_2D
-    import vc_advection_2D
-    import vc_elasticity_2D
-    import vc_acoustics_3D
-    import euler_3D
-    import burgers_3D
-    import vc_advection_3D
+    from . import acoustics_1D
+    from . import acoustics_variable_1D
+    from . import acoustics_1D_ptwise
+    from . import advection_1D
+    from . import advection_1D_ptwise
+    from . import burgers_1D
+    from . import euler_with_efix_1D
+    from . import nonlinear_elasticity_fwave_1D
+    from . import reactive_euler_with_efix_1D
+    from . import shallow_roe_with_efix_1D
+    from . import traffic_1D
+    from . import traffic_vc_1D
+    from . import acoustics_2D
+    from . import acoustics_mapped_2D
+    from . import acoustics_2D_ptwise
+    from . import advection_2D
+    from . import burgers_2D
+    from . import euler_mapgrid_2D
+    from . import euler_5wave_2D
+    from . import euler_4wave_2D
+    from . import kpp_2D
+    from . import psystem_2D
+    from . import shallow_roe_with_efix_2D
+    from . import shallow_sphere_2D
+    from . import vc_acoustics_2D
+    from . import vc_advection_2D
+    from . import vc_elasticity_2D
+    from . import vc_acoustics_3D
+    from . import euler_3D
+    from . import burgers_3D
+    from . import vc_advection_3D
 except ImportError as e:
     import traceback
-    print "********************************************************************"
-    print 'Warning: Some Riemannn solvers were not able to be imported.'
-    print ' Did you run "pip install" in your clawpack directory?'
+    print("********************************************************************")
+    print('Warning: Some Riemannn solvers were not able to be imported.')
+    print(' Did you run "pip install" in your clawpack directory?')
     traceback.print_exc()
-    print "********************************************************************"
+    print("********************************************************************")
 
 import os
 if os.path.exists('./layered_shallow_water_1D.so'):
     import layered_shallow_water_1D
-    import layered_shallow_water_1D_constants
+    from . import layered_shallow_water_1D_constants
diff --git a/src/acoustics_1D_py.py b/src/acoustics_1D_py.py
index 6b6c561f..5be120c3 100644
--- a/src/acoustics_1D_py.py
+++ b/src/acoustics_1D_py.py
@@ -33,6 +33,8 @@
 #                     http://www.opensource.org/licenses/
 # ============================================================================
 
+from __future__ import absolute_import
+from six.moves import range
 num_eqn = 2
 num_waves = 2
 
@@ -78,7 +80,7 @@ def acoustics_1D(q_l,q_r,aux_l,aux_r,problem_data):
     s[1,:] = problem_data['cc']
     
     # Compute the left going and right going fluctuations
-    for m in xrange(num_eqn):
+    for m in range(num_eqn):
         amdq[m,:] = s[0,:] * wave[m,0,:]
         apdq[m,:] = s[1,:] * wave[m,1,:]
     
diff --git a/src/advection_1D_py.py b/src/advection_1D_py.py
index e4054670..930cadcf 100644
--- a/src/advection_1D_py.py
+++ b/src/advection_1D_py.py
@@ -20,6 +20,7 @@
 # ============================================================================
 
 # Riemann solver constants
+from __future__ import absolute_import
 num_eqn = 1
 num_waves = 1
 
diff --git a/src/burgers_1D_py.py b/src/burgers_1D_py.py
index 6682e64a..8e1cb868 100644
--- a/src/burgers_1D_py.py
+++ b/src/burgers_1D_py.py
@@ -17,6 +17,7 @@
 #                     http://www.opensource.org/licenses/
 # ============================================================================
 
+from __future__ import absolute_import
 num_eqn = 1
 num_waves = 1
 
diff --git a/src/euler_1D_py.py b/src/euler_1D_py.py
index d3bc511c..f1c1da9b 100644
--- a/src/euler_1D_py.py
+++ b/src/euler_1D_py.py
@@ -39,7 +39,9 @@
 #                     http://www.opensource.org/licenses/
 # ============================================================================
 
+from __future__ import absolute_import
 import numpy as np
+from six.moves import range
 
 num_eqn = 3
 
@@ -101,8 +103,8 @@ def euler_roe_1D(q_l,q_r,aux_l,aux_r,problem_data):
     else:
         # Godunov update
         s_index = np.zeros((2,num_rp))
-        for m in xrange(num_eqn):
-            for mw in xrange(num_waves):
+        for m in range(num_eqn):
+            for mw in range(num_waves):
                 s_index[0,:] = s[mw,:]
                 amdq[m,:] += np.min(s_index,axis=0) * wave[m,mw,:]
                 apdq[m,:] += np.max(s_index,axis=0) * wave[m,mw,:]
@@ -175,8 +177,8 @@ def euler_hll_1D(q_l,q_r,aux_l,aux_r,problem_data):
     
     # Compute variations
     s_index = np.zeros((2,num_rp))
-    for m in xrange(num_eqn):
-        for mw in xrange(num_waves):
+    for m in range(num_eqn):
+        for mw in range(num_waves):
             s_index[0,:] = s[mw,:]
             amdq[m,:] += np.min(s_index,axis=0) * wave[m,mw,:]
             apdq[m,:] += np.max(s_index,axis=0) * wave[m,mw,:]
@@ -281,8 +283,8 @@ def euler_hllc_1D(q_l,q_r,aux_l,aux_r,problem_data):
     
     # Compute variations
     s_index = np.zeros((2,num_rp))
-    for m in xrange(num_eqn):
-        for mw in xrange(num_waves):
+    for m in range(num_eqn):
+        for mw in range(num_waves):
             s_index[0,:] = s[mw,:]
             amdq[m,:] += np.min(s_index,axis=0) * wave[m,mw,:]
             apdq[m,:] += np.max(s_index,axis=0) * wave[m,mw,:]
diff --git a/src/nonlinear_elasticity_1D_py.py b/src/nonlinear_elasticity_1D_py.py
index 5ab9902f..32eeda08 100644
--- a/src/nonlinear_elasticity_1D_py.py
+++ b/src/nonlinear_elasticity_1D_py.py
@@ -28,7 +28,9 @@
 #                     http://www.opensource.org/licenses/
 # ============================================================================
 
+from __future__ import absolute_import
 import numpy as np
+from six.moves import range
 
 def nonlinear_elasticity_1D(q_l,q_r,aux_l,aux_r,problem_data):
     r"""
@@ -81,7 +83,7 @@ def nonlinear_elasticity_1D(q_l,q_r,aux_l,aux_r,problem_data):
     s[1,:] = cr
     
     # Compute the left going and right going fluctuations
-    for m in xrange(num_eqn):
+    for m in range(num_eqn):
         amdq[m,:] = fwave[m,0,:]
         apdq[m,:] = fwave[m,1,:]
     
diff --git a/src/riemann_interactive.py b/src/riemann_interactive.py
index 9d385577..a09648c9 100644
--- a/src/riemann_interactive.py
+++ b/src/riemann_interactive.py
@@ -11,6 +11,7 @@
 """
 
 
+from __future__ import absolute_import
 import numpy as np
 import jinja2
 import json
diff --git a/src/riemann_tools.py b/src/riemann_tools.py
index c1be7fba..0c1e754e 100644
--- a/src/riemann_tools.py
+++ b/src/riemann_tools.py
@@ -1,3 +1,5 @@
+from __future__ import absolute_import
+from __future__ import print_function
 from matplotlib import animation
 from clawpack.visclaw.JSAnimation import IPython_display
 from IPython.display import display
@@ -6,6 +8,7 @@
 import numpy as np
 import matplotlib.pyplot as plt
 import matplotlib
+from six.moves import range
 
 
 sympy.init_printing(use_latex='mathjax')
@@ -55,19 +58,19 @@ def riemann_solution(solver,q_l,q_r,aux_l=None,aux_r=None,t=0.2,problem_data=Non
     num_states = num_waves + 1
     
     if verbose:
-        print 'States in Riemann solution:'
+        print('States in Riemann solution:')
         states_sym = sympy.Matrix(states)
         display([states_sym[:,k] for k in range(num_states)])
     
-        print 'Waves (jumps between states):'
+        print('Waves (jumps between states):')
         wave_sym = sympy.Matrix(wave[:,:,0])
         display([wave_sym[:,k] for k in range(num_waves)])
     
-        print "Speeds: "
+        print("Speeds: ")
         s_sym = sympy.Matrix(s)
         display(s_sym.T)
     
-        print "amdq, apdq: "
+        print("amdq, apdq: ")
         amdq_sym = sympy.Matrix(amdq).T
         apdq_sym = sympy.Matrix(apdq).T
         display([amdq_sym, apdq_sym])
diff --git a/src/setup.py b/src/setup.py
index 3b4e3441..718f7825 100644
--- a/src/setup.py
+++ b/src/setup.py
@@ -3,6 +3,7 @@
 #
 # f2py -c ../../rp1_layered_shallow_water.f90 -m layered_shallow_water_1D
 
+from __future__ import absolute_import
 one_d_ptwise_riemann = ['acoustics',
                         'advection']
 
diff --git a/src/shallow_1D_py.py b/src/shallow_1D_py.py
index 6a6b9fd9..bc010cdd 100644
--- a/src/shallow_1D_py.py
+++ b/src/shallow_1D_py.py
@@ -36,7 +36,9 @@
 #                     http://www.opensource.org/licenses/
 # ============================================================================
 
+from __future__ import absolute_import
 import numpy as np
+from six.moves import range
 
 num_eqn = 2
 num_waves = 2
@@ -97,8 +99,8 @@ def shallow_roe_1D(q_l,q_r,aux_l,aux_r,problem_data):
         raise NotImplementedError("Entropy fix has not been implemented.")
     else:
         s_index = np.zeros((2,num_rp))
-        for m in xrange(num_eqn):
-            for mw in xrange(num_waves):
+        for m in range(num_eqn):
+            for mw in range(num_waves):
                 s_index[0,:] = s[mw,:]
                 amdq[m,:] += np.min(s_index,axis=0) * wave[m,mw,:]
                 apdq[m,:] += np.max(s_index,axis=0) * wave[m,mw,:]
@@ -165,8 +167,8 @@ def shallow_hll_1D(q_l,q_r,aux_l,aux_r,problem_data):
     
     # Compute variations
     s_index = np.zeros((2,num_rp))
-    for m in xrange(num_eqn):
-        for mw in xrange(num_waves):
+    for m in range(num_eqn):
+        for mw in range(num_waves):
             s_index[0,:] = s[mw,:]
             amdq[m,:] += np.min(s_index,axis=0) * wave[m,mw,:]
             apdq[m,:] += np.max(s_index,axis=0) * wave[m,mw,:]
@@ -222,8 +224,8 @@ def shallow_fwave_1d(q_l, q_r, aux_l, aux_r, problem_data):
     fwave[0,1,:] = beta2
     fwave[1,1,:] = beta2 * s[1,:]
 
-    for m in xrange(num_eqn):
-        for mw in xrange(num_waves):
+    for m in range(num_eqn):
+        for mw in range(num_waves):
             amdq[m,:] += (s[mw,:] < 0.0) * fwave[m,mw,:]
             apdq[m,:] += (s[mw,:] >= 0.0) * fwave[m,mw,:]
 
diff --git a/src/vc_advection_1D_py.py b/src/vc_advection_1D_py.py
index c0a66a21..dc88d23a 100644
--- a/src/vc_advection_1D_py.py
+++ b/src/vc_advection_1D_py.py
@@ -9,6 +9,7 @@
 Note that this is the color equation, not the conservative advection equation.
 """
 
+from __future__ import absolute_import
 import numpy as np
 
 # Riemann solver constants
diff --git a/src/write_constants.py b/src/write_constants.py
index 1c1cbd54..db6463b5 100644
--- a/src/write_constants.py
+++ b/src/write_constants.py
@@ -8,7 +8,10 @@
 - names/indices of q fields
 - names/indices of aux fields
 """
+from __future__ import absolute_import
+from __future__ import print_function
 import os
+from six.moves import zip
 
 for filename in os.listdir('.'):
     # Only process normal Riemann solver files
@@ -53,7 +56,7 @@
         elif filename.startswith('rpn3'):
             out_name = filename.split('.')[0][5:]+'_3D_constants.py'
         out_name = './'+out_name
-        print 'writing ' + out_name
+        print('writing ' + out_name)
         f = open(out_name,'w')
         f.write('num_waves = '+str(num_waves)+'\n')
         f.write('num_eqn   = '+str(num_eqn)+'\n')