diff --git a/CHANGELOG.markdown b/CHANGELOG.markdown
--- a/CHANGELOG.markdown
+++ b/CHANGELOG.markdown
@@ -1,3 +1,9 @@
+0.3.0.0
+-------
+* Breaking change: fix problem with hmatrix API -- replace
+  `HRowMajorMatrix` and `HColumnMajorMatrix` with single `HMatrix`
+  type and manage data ordering internally.
+
 0.2.1.0
 -------
 * Update to work with hmatrix 0.16
diff --git a/Data/NetCDF.hs b/Data/NetCDF.hs
--- a/Data/NetCDF.hs
+++ b/Data/NetCDF.hs
@@ -1,4 +1,4 @@
-{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE ScopedTypeVariables, ConstraintKinds #-}
 -- | Bindings to the Unidata NetCDF data access library.
 --
 --   As well as conventional low-level FFI bindings to the functions
@@ -89,7 +89,7 @@
 -- variables and attributes in the file.
 createFile :: NcInfo NcWrite -> NcIO (NcInfo NcWrite)
 createFile (NcInfo n ds vs as _ _) = runAccess "createFile" n $ do
-  ncid <- chk $ nc_create n (ncIOMode WriteMode)
+  ncid <- chk $ nc_create n ncClobber
   newds <- forM (M.toList ds) (write1Dim ncid . snd)
   let dimids = M.fromList $ zip (M.keys ds) newds
   forM_ (M.toList as) (write1Attr ncid ncGlobal)
@@ -130,7 +130,8 @@
   chk $ get_var1 (ncId nc) ((ncVarIds nc) M.! (ncVarName var)) idxs
 
 -- | Read a whole variable from an open NetCDF file.
-get :: (NcStorable a, NcStore s) => NcInfo NcRead -> NcVar -> NcIO (s a)
+get :: (NcStorable a, NcStore s, NcStoreExtraCon s a) =>
+       NcInfo NcRead -> NcVar -> NcIO (s a)
 get nc var = runAccess "get" (ncName nc) $ do
   let ncid = ncId nc
       varid = (ncVarIds nc) M.! (ncVarName var)
@@ -138,7 +139,7 @@
   chk $ get_var ncid varid sz
 
 -- | Read a slice of a variable from an open NetCDF file.
-getA :: (NcStorable a, NcStore s)
+getA :: (NcStorable a, NcStore s, NcStoreExtraCon s a)
      => NcInfo NcRead -> NcVar -> [Int] -> [Int] -> NcIO (s a)
 getA nc var start count = runAccess "getA" (ncName nc) $ do
   let ncid = ncId nc
@@ -146,7 +147,7 @@
   chk $ get_vara ncid varid start count
 
 -- | Read a strided slice of a variable from an open NetCDF file.
-getS :: (NcStorable a, NcStore s)
+getS :: (NcStorable a, NcStore s, NcStoreExtraCon s a)
      => NcInfo NcRead -> NcVar -> [Int] -> [Int] -> [Int] -> NcIO (s a)
 getS nc var start count stride = runAccess "getS" (ncName nc) $ do
   let ncid = ncId nc
@@ -160,14 +161,15 @@
   chk $ put_var1 (ncId nc) ((ncVarIds nc) M.! (ncVarName var)) idxs val
 
 -- | Write a whole variable to an open NetCDF file.
-put :: (NcStorable a, NcStore s) => NcInfo NcWrite -> NcVar -> s a -> NcIO ()
+put :: (NcStorable a, NcStore s, NcStoreExtraCon s a) =>
+       NcInfo NcWrite -> NcVar -> s a -> NcIO ()
 put nc var val = runAccess "put" (ncName nc) $ do
   let ncid = ncId nc
       varid = (ncVarIds nc) M.! (ncVarName var)
   chk $ put_var ncid varid val
 
 -- | Write a slice of a variable to an open NetCDF file.
-putA :: (NcStorable a, NcStore s)
+putA :: (NcStorable a, NcStore s, NcStoreExtraCon s a)
      => NcInfo NcWrite -> NcVar -> [Int] -> [Int] -> s a -> NcIO ()
 putA nc var start count val = runAccess "putA" (ncName nc) $ do
   let ncid = ncId nc
@@ -175,7 +177,7 @@
   chk $ put_vara ncid varid start count val
 
 -- | Write a strided slice of a variable to an open NetCDF file.
-putS :: (NcStorable a, NcStore s)
+putS :: (NcStorable a, NcStore s, NcStoreExtraCon s a)
      => NcInfo NcWrite -> NcVar -> [Int] -> [Int] -> [Int] -> s a -> NcIO ()
 putS nc var start count stride val = runAccess "putS" (ncName nc) $ do
   let ncid = ncId nc
@@ -264,7 +266,8 @@
 
 -- | Apply COARDS value scaling.
 coardsScale :: forall a b s. (NcStorable a, NcStorable b, FromNcAttr a,
-                              NcStore s, Real a, Fractional b)
+                              NcStore s, Real a, Fractional b,
+                              NcStoreExtraCon s a, NcStoreExtraCon s b)
              => NcVar -> s a -> s b
 coardsScale v din = smap xform din
   where offset = fromMaybe 0.0 $
diff --git a/Data/NetCDF/HMatrix.hs b/Data/NetCDF/HMatrix.hs
--- a/Data/NetCDF/HMatrix.hs
+++ b/Data/NetCDF/HMatrix.hs
@@ -1,48 +1,40 @@
-{-# LANGUAGE FlexibleInstances, MultiParamTypeClasses #-}
-{-# OPTIONS_GHC -fno-warn-orphans #-}
+{-# LANGUAGE FlexibleInstances, ConstraintKinds, TypeFamilies #-}
+{-# LANGUAGE MultiParamTypeClasses, ScopedTypeVariables, GADTs #-}
+{-# OPTIONS_GHC -fno-warn-orphans -fwarn-missing-methods #-}
 -- | NetCDF store instance for HMatrix vectors and matrices.
 module Data.NetCDF.HMatrix
        ( HVector (..)
-       , HRowMajorMatrix (..)
-       , HColumnMajorMatrix (..)
+       , HMatrix (..)
        ) where
 
 import Data.NetCDF.Store
 
 import qualified Numeric.Container as C
+import qualified Numeric.LinearAlgebra.Util as C
 import Foreign.C
-import Data.Packed.Foreign
+import Data.Packed.Vector
 import Data.Packed.Development
 
 newtype HVector a = HVector (C.Vector a)
                   deriving (Eq, Show)
 
 instance NcStore HVector where
-  toForeignPtr (HVector v) = (\(x, _, _) -> x) $ unsafeToForeignPtr v
+  type NcStoreExtraCon HVector a = C.Element a
+  toForeignPtr (HVector v) = fst3 $ unsafeToForeignPtr v
   fromForeignPtr p s = HVector $ unsafeFromForeignPtr p 0 (Prelude.product s)
   smap f (HVector v) = HVector $ C.mapVector f v
 
-newtype HRowMajorMatrix a = HRowMajorMatrix (C.Matrix a)
-
-instance NcStore HRowMajorMatrix where
-  toForeignPtr (HRowMajorMatrix m) = fst $ unsafeMatrixToForeignPtr m
-  fromForeignPtr p s =
-    let c = last s
-        d = product s
-    in HRowMajorMatrix $ matrixFromVector RowMajor (d `div` c) (last s) $
-       unsafeFromForeignPtr p 0 (Prelude.product s)
-  smap f (HRowMajorMatrix m) = HRowMajorMatrix $ C.mapMatrix f m
-
-newtype HColumnMajorMatrix a = HColumnMajorMatrix (C.Matrix a)
+data HMatrix a = HMatrix (C.Matrix a)
 
-instance NcStore HColumnMajorMatrix where
-  toForeignPtr (HColumnMajorMatrix m) = fst $ unsafeMatrixToForeignPtr m
+instance NcStore HMatrix where
+  type NcStoreExtraCon HMatrix a = C.Element a
+  toForeignPtr (HMatrix m) = fst3 $ unsafeToForeignPtr $ C.flatten m
   fromForeignPtr p s =
     let c = last s
         d = product s
-    in HColumnMajorMatrix $ matrixFromVector ColumnMajor (d `div` c) (last s) $
+    in HMatrix $ matrixFromVector RowMajor (d `div` c) (last s) $
        unsafeFromForeignPtr p 0 (Prelude.product s)
-  smap f (HColumnMajorMatrix m) = HColumnMajorMatrix $ C.mapMatrix f m
+  smap f (HMatrix m) = HMatrix $ C.mapMatrix f m
 
 instance C.Element CShort
 instance C.Element CInt
@@ -51,3 +43,9 @@
 
 instance C.Container C.Vector CFloat
 instance C.Container C.Vector CDouble
+
+instance C.Indexable (C.Vector CFloat) CFloat where (!) = (@>)
+instance C.Indexable (C.Vector CDouble) CDouble where (!) = (@>)
+
+fst3 :: (a, b, c) -> a
+fst3 (x, _, _) = x
diff --git a/Data/NetCDF/PutGet.hs b/Data/NetCDF/PutGet.hs
--- a/Data/NetCDF/PutGet.hs
+++ b/Data/NetCDF/PutGet.hs
@@ -1,3 +1,4 @@
+{-# LANGUAGE ConstraintKinds #-}
 -- | Mid-level interface for reading and writing NetCDF data.  The
 -- functions in "Data.NetCDF" provide a more convenient interface.
 
@@ -38,14 +39,16 @@
       v <- peek iv
       return $ (fromIntegral res, v)
 
-put_var :: (NcStorable a, NcStore s) => Int -> Int -> s a -> IO Int
+put_var :: (NcStorable a, NcStore s, NcStoreExtraCon s a) =>
+           Int -> Int -> s a -> IO Int
 put_var nc var v = do
   let ncid = fromIntegral nc
       varid = fromIntegral var
       is = toForeignPtr v
   withForeignPtr is $ \isp -> fromIntegral <$> ffi_put_var ncid varid isp
 
-get_var :: (NcStorable a, NcStore s) => Int -> Int -> [Int] -> IO (Int, s a)
+get_var :: (NcStorable a, NcStore s, NcStoreExtraCon s a) =>
+           Int -> Int -> [Int] -> IO (Int, s a)
 get_var nc var sz = do
   let ncid = fromIntegral nc
       varid = fromIntegral var
@@ -54,7 +57,7 @@
     res <- ffi_get_var ncid varid isp
     return (fromIntegral res, fromForeignPtr is sz)
 
-put_vara :: (NcStorable a, NcStore s) =>
+put_vara :: (NcStorable a, NcStore s, NcStoreExtraCon s a) =>
             Int -> Int -> [Int] -> [Int] -> s a -> IO Int
 put_vara nc var start count v = do
   let ncid = fromIntegral nc
@@ -66,7 +69,7 @@
         res <- ffi_put_vara ncid varid startp countp isp
         return $ fromIntegral res
 
-get_vara :: (NcStorable a, NcStore s)
+get_vara :: (NcStorable a, NcStore s, NcStoreExtraCon s a)
          => Int -> Int -> [Int] -> [Int] -> IO (Int, s a)
 get_vara nc var start count = do
   let ncid = fromIntegral nc
@@ -78,7 +81,7 @@
         res <- ffi_get_vara ncid varid s c isp
         return (fromIntegral res, fromForeignPtr is (filter (>1) count))
 
-put_vars :: (NcStorable a, NcStore s) =>
+put_vars :: (NcStorable a, NcStore s, NcStoreExtraCon s a) =>
             Int -> Int -> [Int] -> [Int] -> [Int] -> s a -> IO Int
 put_vars nc var start count stride v = do
   let ncid = fromIntegral nc
@@ -91,7 +94,7 @@
           res <- ffi_put_vars ncid varid startp countp stridep isp
           return $ fromIntegral res
 
-get_vars :: (NcStorable a, NcStore s)
+get_vars :: (NcStorable a, NcStore s, NcStoreExtraCon s a)
          => Int -> Int -> [Int] -> [Int] -> [Int] -> IO (Int, s a)
 get_vars nc var start count stride = do
   let ncid = fromIntegral nc
diff --git a/Data/NetCDF/Repa.hs b/Data/NetCDF/Repa.hs
--- a/Data/NetCDF/Repa.hs
+++ b/Data/NetCDF/Repa.hs
@@ -1,4 +1,4 @@
-{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE FlexibleInstances, ConstraintKinds, TypeFamilies #-}
 {-# OPTIONS_GHC -fno-warn-orphans #-}
 -- | NetCDF store instance for Repa arrays.
 module Data.NetCDF.Repa where
diff --git a/Data/NetCDF/Store.hs b/Data/NetCDF/Store.hs
--- a/Data/NetCDF/Store.hs
+++ b/Data/NetCDF/Store.hs
@@ -1,4 +1,4 @@
-{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE FlexibleInstances, ConstraintKinds, TypeFamilies #-}
 -- | The /store polymorphism/ for the functions to get the values of
 -- NetCDF variables relies on a simple `NcStore` typeclass for
 -- converting between /store/ values and @ForeignPtr@s.
@@ -7,13 +7,20 @@
 
 import Foreign.Storable
 import Foreign.ForeignPtr
-
+import GHC.Exts
 
 -- | Class representing containers suitable for storing values read
 -- from NetCDF variables.  Just has methods to convert back and forth
 -- between the store and a foreign pointer, and to perform simple
--- mapping over the store.
+-- mapping over the store.  The NcStoreExtraCon associated
+-- constraint-kinded type is used to track extra class constraints on
+-- elements needed for some store types.
 class NcStore s where
-  toForeignPtr :: Storable e => s e -> ForeignPtr e
-  fromForeignPtr :: Storable e => ForeignPtr e -> [Int] -> s e
-  smap :: (Storable a, Storable b) => (a -> b) -> s a -> s b
+  type NcStoreExtraCon s a :: Constraint
+  type NcStoreExtraCon s a = ()
+  toForeignPtr :: (Storable e, NcStoreExtraCon s e) =>
+                  s e -> ForeignPtr e
+  fromForeignPtr :: (Storable e, NcStoreExtraCon s e) =>
+                    ForeignPtr e -> [Int] -> s e
+  smap :: (Storable a, Storable b, NcStoreExtraCon s a, NcStoreExtraCon s b) =>
+          (a -> b) -> s a -> s b
diff --git a/Data/NetCDF/Types.hs b/Data/NetCDF/Types.hs
--- a/Data/NetCDF/Types.hs
+++ b/Data/NetCDF/Types.hs
@@ -62,3 +62,12 @@
 -- | Fake file identifier for unopened files.
 ncInvalidId :: NcId
 ncInvalidId = -1
+
+-- | 'mode' flags for nc_create (see netcdf.h)
+ncClobber, ncNoClobber, nc64bitOffset, ncNetCDF4, ncClassicModel :: Int
+ncClobber      = 0
+ncNoClobber    = 0x0004
+nc64bitOffset  = 0x0200
+ncNetCDF4      = 0x1000
+ncClassicModel = 0x0100
+
diff --git a/Data/NetCDF/Vector.hs b/Data/NetCDF/Vector.hs
--- a/Data/NetCDF/Vector.hs
+++ b/Data/NetCDF/Vector.hs
@@ -1,4 +1,4 @@
-{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE FlexibleInstances, ConstraintKinds, TypeFamilies #-}
 {-# OPTIONS_GHC -fno-warn-orphans #-}
 -- | NetCDF store instance for Storable Vectors.
 module Data.NetCDF.Vector where
diff --git a/examples/example3.hs b/examples/example3.hs
--- a/examples/example3.hs
+++ b/examples/example3.hs
@@ -10,7 +10,7 @@
 import Numeric.Container
 
 type VRet a = IO (Either NcError (HVector a))
-type MRet a = IO (Either NcError (HRowMajorMatrix a))
+type MRet a = IO (Either NcError (HMatrix a))
 
 main :: IO ()
 main = do
@@ -53,9 +53,9 @@
   -- scaling, convert units.
   let (Just zvar) = ncVar nc "z500"
   putStrLn $ "z500 dims: " ++ show (map ncDimName $ ncVarDims zvar)
-  Right slice1tmp <- getA nc zvar [0, 0, 0] [1, nlat, nlon] :: MRet CShort
-  let (HRowMajorMatrix slice1tmp2) =
-        coardsScale zvar slice1tmp :: HRowMajorMatrix CDouble
+  Right slice1tmp <-
+    getA nc zvar [0, 0, 0] [1, nlat, nlon] :: MRet CShort
+  let (HMatrix slice1tmp2) = coardsScale zvar slice1tmp :: HMatrix CDouble
       slice1 = cmap ((/ 9.8) . realToFrac) slice1tmp2 :: Matrix Double
   putStrLn $ "size slice1 = " ++
     show (rows slice1) ++ " x " ++ show (cols slice1)
diff --git a/hnetcdf.cabal b/hnetcdf.cabal
--- a/hnetcdf.cabal
+++ b/hnetcdf.cabal
@@ -1,5 +1,5 @@
 name:                hnetcdf
-version:             0.2.2.0
+version:             0.3.0.0
 synopsis:            Haskell NetCDF library
 description:
   Bindings to the Unidata NetCDF library, along with a higher-level
@@ -131,6 +131,24 @@
                        errors                      >= 1.4.2    && < 1.5,
                        vector                      >= 0.9      && < 0.11,
                        directory                   >= 1.1      && < 1.3
+  extra-libraries:     netcdf
+
+Test-Suite test-hmatrix
+  type:                exitcode-stdio-1.0
+  hs-source-dirs:      test
+  main-is:             test-hmatrix.hs
+  build-depends:       hnetcdf,
+                       base                        >= 4.3      && < 5,
+                       containers                  >= 0.4      && < 0.6,
+                       errors                      >= 1.4.2    && < 1.5,
+                       vector                      >= 0.9      && < 0.11,
+                       directory                   >= 1.1      && < 1.3,
+                       hmatrix                     >= 0.16.0.2 && < 0.17,
+                       test-framework              >= 0.3,
+                       test-framework-hunit,
+                       test-framework-quickcheck2  >= 0.3,
+                       HUnit,
+                       QuickCheck
   extra-libraries:     netcdf
 
 Executable example1
diff --git a/test/test-get.hs b/test/test-get.hs
--- a/test/test-get.hs
+++ b/test/test-get.hs
@@ -78,7 +78,7 @@
 type RepaRet3 a = IO (Either NcError (R.Array F R.DIM3 a))
 type RepaRet1 a = IO (Either NcError (R.Array F R.DIM1 a))
 type HMVRet a = IO (Either NcError (H.HVector a))
-type HMMRet a = IO (Either NcError (H.HRowMajorMatrix a))
+type HMMRet a = IO (Either NcError (H.HMatrix a))
 
 
 --------------------------------------------------------------------------------
@@ -165,7 +165,7 @@
             (and $ zipWith (==) tstvals truevals)
   void $ closeFile nc
 
-getVarHMV :: forall a. (Eq a, Num a, Show a, NcStorable a)
+getVarHMV :: forall a. (Eq a, Num a, Show a, NcStorable a, Element a)
           => FilePath -> String -> a -> Assertion
 getVarHMV f v _ = do
   enc <- openFile f
@@ -291,7 +291,7 @@
               assertBool ("3: value error: " ++ show tstvals ++
                           " instead of " ++ show truevals) (tstvals == truevals)
 
-getVarAHMV :: forall a. (Num a, Show a, Eq a, NcStorable a)
+getVarAHMV :: forall a. (Num a, Show a, Eq a, NcStorable a, Element a)
            => FilePath -> String -> a -> Assertion
 getVarAHMV f v _ = do
   enc <- openFile f
@@ -466,7 +466,7 @@
                             " instead of " ++ show truevals)
                   (tstvals == truevals)
 
-getVarSHMV :: forall a. (Num a, Show a, Eq a, NcStorable a)
+getVarSHMV :: forall a. (Num a, Show a, Eq a, NcStorable a, Element a)
           => FilePath -> String -> a -> Assertion
 getVarSHMV f v _ = do
   enc <- openFile f
diff --git a/test/test-hmatrix.hs b/test/test-hmatrix.hs
new file mode 100644
--- /dev/null
+++ b/test/test-hmatrix.hs
@@ -0,0 +1,67 @@
+{-# LANGUAGE ScopedTypeVariables #-}
+
+import Test.Framework (defaultMain, testGroup, Test)
+import Test.Framework.Providers.HUnit
+import Test.Framework.Providers.QuickCheck2
+import Test.HUnit hiding (Test, assert)
+import Test.QuickCheck
+import Test.QuickCheck.Monadic
+
+import qualified Data.Map as M
+import Data.Maybe
+import qualified Data.Vector.Storable as SV
+import Numeric.LinearAlgebra.HMatrix
+import Control.Error
+import Control.Monad
+import Foreign.C
+
+import Data.NetCDF
+import Data.NetCDF.Vector
+import Data.NetCDF.HMatrix
+
+type HMatrixRet a = IO (Either NcError (HMatrix a))
+
+main :: IO ()
+main = defaultMain tests
+
+tests :: [Test]
+tests =
+  [ testGroup "HMatrix handling"
+    [ testCase "Data ordering" hmatrixDataOrdering
+    ]
+  ]
+
+hmatrixDataOrdering :: Assertion
+hmatrixDataOrdering = do
+  let xdim = NcDim "x" 5 False
+      ydim = NcDim "y" 10 False
+      ncout = emptyNcInfo "tst-hmatrix.nc" #
+           addNcDim xdim #
+           addNcDim ydim #
+           addNcVar (NcVar "x" NcDouble [xdim] M.empty) #
+           addNcVar (NcVar "y" NcDouble [ydim] M.empty) #
+           addNcVar (NcVar "z1" NcDouble [ydim, xdim] M.empty) #
+           addNcVar (NcVar "z2" NcDouble [ydim, xdim] M.empty)
+      z1 = cmap realToFrac $ matrix 5 [1..50] :: Matrix CDouble
+      z2 = fromColumns $ toColumns z1 :: Matrix CDouble
+  let write nc = do
+        let xvar = fromJust $ ncVar nc "x"
+            yvar = fromJust $ ncVar nc "y"
+            z1var = fromJust $ ncVar nc "z1"
+            z2var = fromJust $ ncVar nc "z2"
+        put nc xvar $ SV.fromList [1..5 :: CDouble]
+        put nc yvar $ SV.fromList [1..10 :: CDouble]
+        put nc z1var $ HMatrix z1
+        put nc z2var $ HMatrix z2
+        return ()
+  withCreateFile ncout write (putStrLn . ("ERROR: " ++) . show)
+  let read nc = do
+        let z1var = fromJust $ ncVar nc "z1"
+            z2var = fromJust $ ncVar nc "z2"
+        Right (HMatrix z1) <- get nc z1var :: HMatrixRet CDouble
+        Right (HMatrix z2) <- get nc z2var :: HMatrixRet CDouble
+        return (z1, z2)
+  (z1tst, z2tst) <-
+    withReadFile "tst-hmatrix.nc" read (error . ("ERROR: " ++) . show)
+  assertBool "Straightforward write/read" (flatten z1tst == flatten z1)
+  assertBool "Columnar write/read" (flatten z2tst == flatten z2)