diff --git a/src/boxclaw/__init__.py b/src/boxclaw/__init__.py
index 401dc5173..3e346d3c2 100644
--- a/src/boxclaw/__init__.py
+++ b/src/boxclaw/__init__.py
@@ -3,6 +3,9 @@
 import os
 import logging, logging.config
 
+import fboxlib
+fboxlib.open()
+
 # Default logging configuration file
 _DEFAULT_LOG_CONFIG_PATH = os.path.join(os.path.dirname(__file__),'log.config')
 del os
diff --git a/src/boxclaw/cfl.py b/src/boxclaw/cfl.py
index 9142492f6..5a0878eb5 100644
--- a/src/boxclaw/cfl.py
+++ b/src/boxclaw/cfl.py
@@ -1,4 +1,6 @@
 
+import fboxlib
+
 class CFL(object):
     """Parallel CFL object, responsible for computing the
     Courant-Friedrichs-Lewy condition across all processes.
@@ -10,14 +12,14 @@ def __init__(self, global_max):
 
     def get_global_max(self):
         r"""
-        Compute the maximum CFL number over all processes for the current step.
+        Compute the maximum CFL number over all processes for the
+        current step.
 
         This is used to determine whether the CFL condition was
         violated and adjust the timestep.
         """
-        import boxlib
         self._reduce_vec.array = self._local_max
-        self._global_max = boxlib.bl[0].ReduceRealMax(self._local_max)
+        self._global_max = fboxlib.reduce_max(self._local_max)
         return self._global_max
 
     def get_cached_max(self):
@@ -27,6 +29,4 @@ def set_local_max(self,new_local_max):
         self._local_max = new_local_max
 
     def update_global_max(self,new_local_max):
-        import boxlib
-        self._global_max = boxlib.bl[0].ReduceRealMax(new_local_max)
-
+        self._global_max = fboxlib.reduce_max(new_local_max)
diff --git a/src/boxclaw/controller.py b/src/boxclaw/controller.py
index 546099e50..ae0565c8b 100644
--- a/src/boxclaw/controller.py
+++ b/src/boxclaw/controller.py
@@ -2,6 +2,7 @@
 Module for BoxClaw controller class.
 """
 
+import fboxlib
 from clawpack import pyclaw
 
 class Controller(pyclaw.controller.Controller):
@@ -14,8 +15,7 @@ def __init__(self):
         self.output_format = 'multifab'
 
     def is_proc_0(self):
-        import boxlib
-        return boxlib.rank() == 0
+        return fboxlib.mpi_rank() == 0
 
     def log_info(self, str):
         import logging
diff --git a/src/boxclaw/geometry.py b/src/boxclaw/geometry.py
index 26658c6a9..73978385f 100644
--- a/src/boxclaw/geometry.py
+++ b/src/boxclaw/geometry.py
@@ -4,6 +4,9 @@
 Module containing BoxClaw geometry.
 """
 
+import fboxlib
+import math
+
 from clawpack import pyclaw
 from clawpack.pyclaw import geometry as pyclaw_geometry
 
@@ -16,35 +19,17 @@ class Patch(pyclaw_geometry.Patch):
     __doc__ += pyclaw.util.add_parent_doc(pyclaw_geometry.Patch)
 
     def __deepcopy__(self,memo):
-        # don't recreate boxarrays
         return self
 
     def __init__(self,dimensions):
 
-        import boxlib
-        import numpy as np
-
         super(Patch,self).__init__(dimensions)
 
-        self._ba, self._gbox = self._create_boxarray()
-
-        is_per = np.asarray(self.num_dim * [ 1 ], np.int32)
-        rb     = boxlib.RealBox(boxlib.lo(self._gbox), boxlib.hi(self._gbox))
-        self._geom = boxlib.bl[self.num_dim].Geometry(self._gbox, rb, 0, is_per)
+        self._la = self._create_layout()
 
-        # XXX: create a multifab from the boxarray to get geometry information
-        tmp = boxlib.MultiFab(self._ba)
+        bxs = self._la.local_boxes
+        lo, hi = self._la.get_box(bxs[0])
 
-        fab = None
-        for i in range(tmp.size()):
-            fab = tmp[i]
-            if fab is not None:
-                break
-
-        assert(fab is not None)
-
-        lo = boxlib.lo(fab.box())
-        hi = boxlib.hi(fab.box())
         grid_dimensions = []
         for i in range(self.num_dim):
             lower = (lo[i]+0) * self.delta[i]
@@ -52,7 +37,7 @@ def __init__(self,dimensions):
             num_cells = hi[i]-lo[i]+1
 
             grid_dimensions.append(pyclaw_geometry.Dimension(lower,upper,
-                                        num_cells,name=dimensions[i].name))
+                                                             num_cells,name=dimensions[i].name))
 
             if lower == self.lower_global[i]:
                 grid_dimensions[-1].on_lower_boundary = True
@@ -66,17 +51,9 @@ def __init__(self,dimensions):
 
         self.grid = pyclaw_geometry.Grid(grid_dimensions)
 
-        del tmp
-
 
-    def _create_boxarray(self):
-        """Returns a BoxLib BoxArray."""
-
-        # note that boxlib supports more than one box per processor
-
-        import boxlib
-        import math
-        import numpy as np
+    def _create_layout(self):
+        """Returns a FBoxLib layout."""
 
         dx = self.delta
         lg = self.lower_global
@@ -85,12 +62,8 @@ def _create_boxarray(self):
         lo = [ int(lg[d]/dx[d]) for d in range(self.num_dim) ]
         hi = [ lo[d]+nc[d]-1    for d in range(self.num_dim) ]
 
-        box = boxlib.Box(lo, hi)
-        ba  = boxlib.BoxArray([box])
-
         max_sizes = self.num_dim * [ 1 ]
-
-        nprocs = boxlib.size()
+        nprocs = fboxlib.mpi_size()
         if (self.num_dim > 1) and (nprocs % 2**self.num_dim == 0):
             # divide domain into cubes
             nproc_per_dim = nprocs / 2**self.num_dim + 1
@@ -102,14 +75,10 @@ def _create_boxarray(self):
             for d in range(1, self.num_dim):
                 max_sizes[d] = hi[d] - lo[d] + 1
 
-        max_sizes = boxlib.bl[self.num_dim].IntVect(*max_sizes)
-        ba.maxSize(max_sizes)
-
-        if boxlib.rank() == 0:
-            logger.info("max_sizes: " + str(max_sizes))
-            logger.info("boxarray:  " + str(ba))
-
-        return ba, box
+        ba = fboxlib.boxarray(boxes=[[lo, hi]])
+        ba.maxsize(max_sizes)
+        la = fboxlib.layout(ba)
+        return la
 
 
 class Domain(pyclaw_geometry.Domain):
diff --git a/src/boxclaw/io/multifab.py b/src/boxclaw/io/multifab.py
index 327a3c1ea..ed588afcf 100644
--- a/src/boxclaw/io/multifab.py
+++ b/src/boxclaw/io/multifab.py
@@ -8,7 +8,7 @@
 frame data.
 """
 
-import boxlib
+import fboxlib
 import pickle
 import numpy as np
 
@@ -56,7 +56,7 @@ def write(solution,frame,path='./',file_prefix='claw',write_aux=False,
             if os.path.exists(f):
                 raise IOError('Cowardly refusing to clobber %s!' % f)
 
-    rank = boxlib.rank()
+    rank = fboxlib.mpi_rank()
     if rank==0:
         pickle_file = open(pickle_filename,'wb')
         # explicitly dumping a dictionary here to help out anybody trying to read the pickle file
@@ -119,27 +119,25 @@ def write(solution,frame,path='./',file_prefix='claw',write_aux=False,
 
                 if finest_level == 0:
                     write('')
-                write(solution.state.patch._gbox)       # grid domain
+                lo, hi = solution.domain.patch._la.get_box(1)
+                write(str(lo) + " " + str(hi))          # grid domain XXX
+                write('')
                 write(frame)                            # time step number
                 write(' '.join(map(str, dx)))           # dx
                 write('0')                              # cartesian coords
                 write('0')                              # boundary data? (0=no)
 
                 # now write info for each amr level (only one in this case)
-                write('0 %d %f' % (mf.size(), solution.state.t))
+                write('0 %d %f' % (mf.nboxes, solution.state.t))
                 write(frame)
-                for i in range(mf.size()):
-                    bx = solution.state.patch._ba.get(i)
-                    lo = bx.smallEnd()
-                    hi = bx.bigEnd()
+                for i in range(1, mf.nboxes+1):
+                    lo, hi = solution.state.patch._la.get_box(i)
                     for d in range(solution.state.num_dim):
                         write("%lf %lf" % (lo[d]*dx[d], (hi[d]+1)*dx[d]))
                 write("Level_0/Cell")
 
         # write cell data
-        mf.writeOut(plt_filename + '/Level_0/Cell')
+        mf.write(plt_filename, 'Level_0/Cell')
 
     if rank==0:
         pickle_file.close()
-
-
diff --git a/src/boxclaw/state.py b/src/boxclaw/state.py
index 02c47d0d5..8cb1faeaa 100644
--- a/src/boxclaw/state.py
+++ b/src/boxclaw/state.py
@@ -1,4 +1,5 @@
 
+import fboxlib
 import clawpack.pyclaw
 import numpy as np
 
@@ -19,20 +20,15 @@ def get_auxbc_from_aux(self,num_ghost,auxbc):
     def _fill_boundary(self, ncomp, nghost, q, qbc):
         mf = self._create_multifab(ncomp, nghost)
         self._copy_into_multifab(mf, nghost, q)
-        mf.FillBoundary()
-        self.patch._geom.FillPeriodicBoundary(mf)
+        mf.fill_boundary()
         self._copy_outof_multifab(mf, qbc)
         del mf
 
     def _create_multifab(self, ncomp, nghost):
-        import boxlib
-        return boxlib.MultiFab(self.patch._ba, ncomp, nghost)
+        return fboxlib.multifab(self.patch._la, ncomp, nghost)
 
     def _get_array(self, mf):
-        for i in range(mf.size()):
-            if mf[i] is not None:
-                return mf[i].get_array()
-        return None
+        return mf.fab(1).array
 
     def _copy_into_multifab(self, mf, ng, q):
 
@@ -51,12 +47,5 @@ def _copy_into_multifab(self, mf, ng, q):
             raise Exception("Assumption (1 <= num_dim <= 3) violated.")
 
     def _copy_outof_multifab(self, mf, q):
-        num_dim = self.patch.num_dim
-        fab     = self._get_array(mf)
+        fab = self._get_array(mf)
         q[...] = np.rollaxis(fab, self.q.ndim-1)
-
-
-
-
-
-
diff --git a/src/boxclaw/util.py b/src/boxclaw/util.py
index d8317dfc0..d9e278d83 100644
--- a/src/boxclaw/util.py
+++ b/src/boxclaw/util.py
@@ -5,9 +5,9 @@
 def boxlib_build_variant_arg_dicts(kernel_languages=('Fortran',)):
     import itertools
 
-    # test petsc4py only if it is available
+    # test fboxlib only if it is available
     try:
-        import boxlib
+        import fboxlib
     except ImportError:
         return []