diff --git a/src/plot.F90 b/src/plot.F90
index 075061bebf..d2e81308f5 100644
--- a/src/plot.F90
+++ b/src/plot.F90
@@ -71,108 +71,6 @@ contains
     err = 0
   end function openmc_plot_geometry
 
-!===============================================================================
-! DRAW_MESH_LINES draws mesh line boundaries on an image
-!===============================================================================
-
-  subroutine draw_mesh_lines(pl, data)
-    type(ObjectPlot), intent(in)    :: pl
-    integer,          intent(inout) :: data(:,:,:)
-
-    logical :: in_mesh
-    integer :: out_, in_  ! pixel location
-    integer :: rgb(3)     ! RGB color for meshlines pixels
-    integer :: outrange(2), inrange(2) ! range of pixel locations
-    integer :: i, j       ! loop indices
-    integer :: plus
-    integer :: ijk_ll(3)  ! mesh bin ijk indicies of plot lower left
-    integer :: ijk_ur(3)  ! mesh bin ijk indicies of plot upper right
-    integer :: outer, inner
-    real(8) :: frac
-    real(8) :: width(3)   ! real widths of the plot
-    real(8) :: xyz_ll_plot(3)  ! lower left xyz of plot image
-    real(8) :: xyz_ur_plot(3)  ! upper right xyz of plot image
-    real(8) :: xyz_ll(3)  ! lower left xyz
-    real(8) :: xyz_ur(3)  ! upper right xyz
-    type(RegularMesh) :: m
-
-    rgb(:) = pl % meshlines_color % rgb
-
-    select case (pl % basis)
-      case(PLOT_BASIS_XY)
-        outer = 1
-        inner = 2
-      case(PLOT_BASIS_XZ)
-        outer = 1
-        inner = 3
-      case(PLOT_BASIS_YZ)
-        outer = 2
-        inner = 3
-    end select
-
-    xyz_ll_plot = pl % origin
-    xyz_ur_plot = pl % origin
-
-    xyz_ll_plot(outer) = pl % origin(outer) - pl % width(1) / TWO
-    xyz_ll_plot(inner) = pl % origin(inner) - pl % width(2) / TWO
-    xyz_ur_plot(outer) = pl % origin(outer) + pl % width(1) / TWO
-    xyz_ur_plot(inner) = pl % origin(inner) + pl % width(2) / TWO
-
-    width = xyz_ur_plot - xyz_ll_plot
-
-    m = meshes(pl % index_meshlines_mesh)
-    call m % get_indices(xyz_ll_plot, ijk_ll(:m % n_dimension()), in_mesh)
-    call m % get_indices(xyz_ur_plot, ijk_ur(:m % n_dimension()), in_mesh)
-
-    ! sweep through all meshbins on this plane and draw borders
-    do i = ijk_ll(outer), ijk_ur(outer)
-      do j = ijk_ll(inner), ijk_ur(inner)
-        ! check if we're in the mesh for this ijk
-        if (i > 0 .and. i <= m % dimension(outer) .and. &
-              j > 0 .and. j <= m % dimension(inner)) then
-
-          ! get xyz's of lower left and upper right of this mesh cell
-          xyz_ll(outer) = m % lower_left(outer) + m % width(outer) * (i - 1)
-          xyz_ll(inner) = m % lower_left(inner) + m % width(inner) * (j - 1)
-          xyz_ur(outer) = m % lower_left(outer) + m % width(outer) * i
-          xyz_ur(inner) = m % lower_left(inner) + m % width(inner) * j
-
-          ! map the xyz ranges to pixel ranges
-
-          frac = (xyz_ll(outer) - xyz_ll_plot(outer)) / width(outer)
-          outrange(1) = int(frac * real(pl % pixels(1), 8))
-          frac = (xyz_ur(outer) - xyz_ll_plot(outer)) / width(outer)
-          outrange(2) = int(frac * real(pl % pixels(1), 8))
-
-          frac = (xyz_ur(inner) - xyz_ll_plot(inner)) / width(inner)
-          inrange(1) = int((ONE - frac) * real(pl % pixels(2), 8))
-          frac = (xyz_ll(inner) - xyz_ll_plot(inner)) / width(inner)
-          inrange(2) = int((ONE - frac) * real(pl % pixels(2), 8))
-
-          ! draw lines
-          do out_ = outrange(1), outrange(2)
-            do plus = 0, pl % meshlines_width
-              data(:, out_ + 1, inrange(1) + plus + 1) = rgb
-              data(:, out_ + 1, inrange(2) + plus + 1) = rgb
-              data(:, out_ + 1, inrange(1) - plus + 1) = rgb
-              data(:, out_ + 1, inrange(2) - plus + 1) = rgb
-            end do
-          end do
-          do in_ = inrange(1), inrange(2)
-            do plus = 0, pl % meshlines_width
-              data(:, outrange(1) + plus + 1, in_ + 1) = rgb
-              data(:, outrange(2) + plus + 1, in_ + 1) = rgb
-              data(:, outrange(1) - plus + 1, in_ + 1) = rgb
-              data(:, outrange(2) - plus + 1, in_ + 1) = rgb
-            end do
-          end do
-
-        end if
-      end do
-    end do
-
-  end subroutine draw_mesh_lines
-
 !===============================================================================
 ! CREATE_VOXEL outputs a binary file that can be input into silomesh for 3D
 ! geometry visualization.  It works the same way as create_ppm by dragging a