diff --git a/CHANGELOG.markdown b/CHANGELOG.markdown
new file mode 100644
--- /dev/null
+++ b/CHANGELOG.markdown
@@ -0,0 +1,4 @@
+0.2.0.0
+-------
+* Add partial output API
+
diff --git a/Data/NetCDF.hs b/Data/NetCDF.hs
--- a/Data/NetCDF.hs
+++ b/Data/NetCDF.hs
@@ -1,3 +1,4 @@
+{-# LANGUAGE ScopedTypeVariables #-}
 -- | Bindings to the Unidata NetCDF data access library.
 --
 --   As well as conventional low-level FFI bindings to the functions
@@ -16,7 +17,7 @@
 -- > ...
 -- > type SVRet = IO (Either NcError (SV.Vector a))
 -- > ...
--- >   enc <- openFile "tst.nc" ReadMode
+-- >   enc <- openFile "tst.nc"
 -- >   case enc of
 -- >     Right nc -> do
 -- >       eval <- get nc "varname" :: SVRet CDouble
@@ -34,7 +35,7 @@
 -- > ...
 -- > type RepaRet3 a = IO (Either NcError (R.Array F R.DIM3 a))
 -- > ...
--- >   enc <- openFile "tst.nc" ReadMode
+-- >   enc <- openFile "tst.nc"
 -- >   case enc of
 -- >     Right nc -> do
 -- >       eval <- get nc "varname" :: RepaRet3 CDouble
@@ -45,8 +46,11 @@
        , module Data.NetCDF.Metadata
        , NcStorable (..)
        , IOMode (..)
-       , openFile, closeFile, withFile
-       , get1, get, getA, getS ) where
+       , openFile, createFile, closeFile
+       , withReadFile, withCreateFile
+       , get1, get, getA, getS
+       , put1, put, putA, putS
+       , coardsScale ) where
 
 import Data.NetCDF.Raw
 import Data.NetCDF.Types
@@ -56,52 +60,73 @@
 import Data.NetCDF.Store
 import Data.NetCDF.Utils
 
-import Control.Applicative ((<$>))
 import Control.Exception (bracket)
-import Control.Monad (forM, void)
-import Control.Monad.IO.Class
-import Control.Monad.Trans.Reader
+import Control.Monad (forM, forM_, void)
 import Control.Error
-import Data.List
 import qualified Data.Map as M
 import Foreign.C
 import System.IO (IOMode (..))
 
--- | Open a NetCDF file and read all metadata.
-openFile :: FilePath -> IOMode -> NcIO NcInfo
-openFile p mode = runAccess "openFile" p $ do
-  ncid <- chk $ nc_open p (ncIOMode mode)
+-- | Open an existing NetCDF file for read-only access and read all
+-- metadata: the returned 'NcInfo' value contains all the information
+-- about dimensions, variables and attributes in the file.
+openFile :: FilePath -> NcIO (NcInfo NcRead)
+openFile p = runAccess "openFile" p $ do
+  ncid <- chk $ nc_open p (ncIOMode ReadMode)
   (ndims, nvars, nattrs, unlim) <- chk $ nc_inq ncid
   dims <- forM [0..ndims-1] (read1Dim ncid unlim)
   attrs <- forM [0..nattrs-1] (read1Attr ncid ncGlobal)
   vars <- forM [0..nvars-1] (read1Var ncid dims)
   let mkMap nf = foldl (\m v -> M.insert (nf v) v m) M.empty
       dimmap = mkMap ncDimName dims
-      attmap = mkMap ncAttrName attrs
+      attmap = M.fromList attrs
       varmap = mkMap ncVarName vars
       varidmap = M.fromList $ zip (map ncVarName vars) [0..]
   return $ NcInfo p dimmap varmap attmap ncid varidmap
 
+-- | Create a new NetCDF file, ready for write-only access.  The
+-- 'NcInfo' parameter contains all the information about dimensions,
+-- 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)
+  newds <- forM (M.toList ds) (write1Dim ncid . snd)
+  let dimids = M.fromList $ zip (M.keys ds) newds
+  forM_ (M.toList as) (write1Attr ncid ncGlobal)
+  newvs <- forM (M.toList vs) (write1Var ncid dimids . snd)
+  let varids = M.fromList $ zip (M.keys vs) newvs
+  chk $ nc_enddef ncid
+  return $ NcInfo n ds vs as ncid varids
+
 -- | Close a NetCDF file.
-closeFile :: NcInfo -> IO ()
+closeFile :: NcInfo a -> IO ()
 closeFile (NcInfo _ _ _ _ ncid _) = void $ nc_close ncid
 
--- | Bracket file use: a little different from the standard 'withFile'
--- function because of error handling.
-withFile :: FilePath -> IOMode
-         -> (NcInfo -> IO r) -> (NcError -> IO r) -> IO r
-withFile p m ok err = bracket
-                      (openFile p m)
+-- | Bracket read-only file use: a little different from the standard
+-- 'withFile' function because of error handling.
+withReadFile :: FilePath
+             -> (NcInfo NcRead -> IO r) -> (NcError -> IO r) -> IO r
+withReadFile p ok e = bracket
+                      (openFile p)
                       (either (const $ return ()) closeFile)
-                      (either err ok)
+                      (either e ok)
 
+-- | Bracket write-only file use: a little different from the standard
+-- 'withFile' function because of error handling.
+withCreateFile :: NcInfo NcWrite
+               -> (NcInfo NcWrite -> IO r) -> (NcError -> IO r) -> IO r
+withCreateFile nc ok e = bracket
+                         (createFile nc)
+                         (either (const $ return ()) closeFile)
+                         (either e ok)
+
 -- | Read a single value from an open NetCDF file.
-get1 :: NcStorable a => NcInfo -> NcVar -> [Int] -> NcIO a
+get1 :: NcStorable a => NcInfo NcRead -> NcVar -> [Int] -> NcIO a
 get1 nc var idxs = runAccess "get1" (ncName nc) $
   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 -> NcVar -> NcIO (s a)
+get :: (NcStorable a, NcStore s) => NcInfo NcRead -> NcVar -> NcIO (s a)
 get nc var = runAccess "get" (ncName nc) $ do
   let ncid = ncId nc
       varid = (ncVarIds nc) M.! (ncVarName var)
@@ -110,7 +135,7 @@
 
 -- | Read a slice of a variable from an open NetCDF file.
 getA :: (NcStorable a, NcStore s)
-     => NcInfo -> NcVar -> [Int] -> [Int] -> NcIO (s a)
+     => NcInfo NcRead -> NcVar -> [Int] -> [Int] -> NcIO (s a)
 getA nc var start count = runAccess "getA" (ncName nc) $ do
   let ncid = ncId nc
       varid = (ncVarIds nc) M.! (ncVarName var)
@@ -118,49 +143,131 @@
 
 -- | Read a strided slice of a variable from an open NetCDF file.
 getS :: (NcStorable a, NcStore s)
-     => NcInfo -> NcVar -> [Int] -> [Int] -> [Int] -> NcIO (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
       varid = (ncVarIds nc) M.! (ncVarName var)
   chk $ get_vars ncid varid start count stride
 
 
+-- | Write a single value to an open NetCDF file.
+put1 :: NcStorable a => NcInfo NcWrite -> NcVar -> [Int] -> a -> NcIO ()
+put1 nc var idxs val = runAccess "put1" (ncName nc) $
+  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 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)
+     => NcInfo NcWrite -> NcVar -> [Int] -> [Int] -> s a -> NcIO ()
+putA nc var start count val = runAccess "putA" (ncName nc) $ do
+  let ncid = ncId nc
+      varid = (ncVarIds nc) M.! (ncVarName var)
+  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)
+     => 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
+      varid = (ncVarIds nc) M.! (ncVarName var)
+  chk $ put_vars ncid varid start count stride val
+
+
 -- | Helper function to read a single NC dimension.
 read1Dim :: Int -> Int -> Int -> Access NcDim
 read1Dim ncid unlim dimid = do
   (name, len) <- chk $ nc_inq_dim ncid dimid
   return $ NcDim name len (dimid == unlim)
 
+-- | Helper function to write a single NC dimension.
+write1Dim :: Int -> NcDim -> Access Int
+write1Dim ncid (NcDim name len unlim) = do
+  chk $ nc_def_dim ncid name (if unlim then ncUnlimitedLength else len)
+
 -- | Helper function to read a single NC attribute.
-read1Attr :: Int -> Int -> Int -> Access NcAttr
+read1Attr :: Int -> Int -> Int -> Access (Name, NcAttr)
 read1Attr ncid varid attid = do
   n <- chk $ nc_inq_attname ncid varid attid
   (itype, len) <- chk $ nc_inq_att ncid varid n
   a <- readAttr ncid varid n (toEnum itype) len
-  return a
+  return (n, a)
 
+-- | Helper function to write a single NC attribute.
+write1Attr :: Int -> Int -> (Name, NcAttr) -> Access ()
+write1Attr ncid varid (n, a) = writeAttr ncid varid n a
+
 -- | Helper function to read metadata for a single NC variable.
 read1Var :: Int -> [NcDim] -> Int -> Access NcVar
 read1Var ncid dims varid = do
   (n, itype, nvdims, vdimids, nvatts) <- chk $ nc_inq_var ncid varid
   let vdims = map (dims !!) $ take nvdims vdimids
   vattrs <- forM [0..nvatts-1] (read1Attr ncid varid)
-  let vattmap = foldl (\m a -> M.insert (ncAttrName a) a m) M.empty vattrs
+  let vattmap = foldl (\m (nm, a) -> M.insert nm a m) M.empty vattrs
   return $ NcVar n (toEnum itype) vdims vattmap
 
+-- | Helper function to write metadata for a single NC variable.
+write1Var :: Int -> M.Map Name Int -> NcVar -> Access Int
+write1Var ncid dimidmap (NcVar n t dims as) = do
+  let dimids = map ((dimidmap M.!) . ncDimName) dims
+  varid <- chk $ nc_def_var ncid n (fromEnum t) (length dims) dimids
+  forM_ (M.toList as) (write1Attr ncid varid)
+  return varid
+
 -- | Read an attribute from a NetCDF variable with error handling.
 readAttr :: Int -> Int -> String -> NcType -> Int -> Access NcAttr
-readAttr nc var n NcChar l = readAttr' nc var n l nc_get_att_text
-readAttr nc var n NcShort l = readAttr' nc var n l nc_get_att_short
-readAttr nc var n NcInt l = readAttr' nc var n l nc_get_att_int
-readAttr nc var n NcFloat l = readAttr' nc var n l nc_get_att_float
-readAttr nc var n NcDouble l = readAttr' nc var n l nc_get_att_double
-readAttr nc var n NcUInt l = readAttr' nc var n l nc_get_att_uint
-readAttr _ _ n _ _ = return $ NcAttr n ([0] :: [CInt])
+readAttr nc var n NcByte l =
+  readAttr' nc var n l (NcAttrByte . map fromIntegral) nc_get_att_uchar
+readAttr nc var n NcChar l = readAttr' nc var n l NcAttrChar nc_get_att_text
+readAttr nc var n NcShort l = readAttr' nc var n l NcAttrShort nc_get_att_short
+readAttr nc var n NcInt l = readAttr' nc var n l NcAttrInt nc_get_att_int
+readAttr nc var n NcFloat l = readAttr' nc var n l NcAttrFloat nc_get_att_float
+readAttr nc var n NcDouble l =
+  readAttr' nc var n l NcAttrDouble nc_get_att_double
 
 -- | Helper function for attribute reading.
-readAttr' :: Show a => Int -> Int -> String -> Int
+readAttr' :: Show a => Int -> Int -> String -> Int -> ([a] -> NcAttr)
           -> (Int -> Int -> String -> Int -> IO (Int, [a])) -> Access NcAttr
-readAttr' nc var n l rf = NcAttr n <$> (chk $ rf nc var n l)
+readAttr' nc var n l w rf = chk $ do
+  tmp <- rf nc var n l
+  return $ (fst tmp, w $ snd tmp)
 
+-- | Write an attribute to a NetCDF variable with error handling.
+writeAttr :: Int -> Int -> String -> NcAttr -> Access ()
+writeAttr nc var n (NcAttrByte v) = writeAttr' nc var n id v nc_put_att_uchar
+writeAttr nc var n (NcAttrChar v) = writeAttr' nc var n id v nc_put_att_text
+writeAttr nc var n (NcAttrShort v) =
+  writeAttr' nc var n fromIntegral v nc_put_att_short
+writeAttr nc var n (NcAttrInt v) =
+  writeAttr' nc var n fromIntegral v nc_put_att_int
+writeAttr nc var n (NcAttrFloat v) =
+  writeAttr' nc var n realToFrac v nc_put_att_float
+writeAttr nc var n (NcAttrDouble v) =
+  writeAttr' nc var n realToFrac v nc_put_att_double
 
+-- | Helper function for attribute writeing.
+writeAttr' :: Int -> Int -> String -> (a -> b) -> [a]
+          -> (Int -> Int -> String -> Int -> [b] -> IO Int) -> Access ()
+writeAttr' nc var n conv vs wf = chk $ wf nc var n (length vs) (map conv vs)
+
+
+-- | Apply COARDS value scaling.
+coardsScale :: forall a b s. (NcStorable a, NcStorable b, FromNcAttr a,
+                              NcStore s, Real a, Fractional b)
+             => NcVar -> s a -> s b
+coardsScale v din = smap xform din
+  where offset = fromMaybe 0.0 $
+                 ncVarAttr v "add_offset" >>= fromAttr :: CDouble
+        scale = fromMaybe 1.0 $
+                ncVarAttr v "scale_factor" >>= fromAttr :: CDouble
+        fill = ncVarAttr v "_FillValue" >>= fromAttr :: Maybe a
+        xform x = case fill of
+          Nothing -> realToFrac $ realToFrac x * scale + offset
+          Just f -> if x == f
+                    then realToFrac f
+                    else realToFrac $ realToFrac x * scale + offset
diff --git a/Data/NetCDF/HMatrix.hs b/Data/NetCDF/HMatrix.hs
new file mode 100644
--- /dev/null
+++ b/Data/NetCDF/HMatrix.hs
@@ -0,0 +1,40 @@
+{-# LANGUAGE FlexibleInstances #-}
+{-# OPTIONS_GHC -fno-warn-orphans #-}
+-- | NetCDF store instance for HMatrix vectors and matrices.
+module Data.NetCDF.HMatrix
+       ( HVector (..)
+       , HRowMajorMatrix (..)
+       , HColumnMajorMatrix (..)
+       ) where
+
+import Data.NetCDF.Store
+
+import qualified Numeric.Container as C
+import Data.Packed.Foreign
+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
+  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 = HRowMajorMatrix $
+                       matrixFromVector RowMajor (last s) $
+                       unsafeFromForeignPtr p 0 (Prelude.product s)
+  smap f (HRowMajorMatrix m) = HRowMajorMatrix $ C.mapMatrix f m
+
+newtype HColumnMajorMatrix a = HColumnMajorMatrix (C.Matrix a)
+
+instance NcStore HColumnMajorMatrix where
+  toForeignPtr (HColumnMajorMatrix m) = fst $ unsafeMatrixToForeignPtr m
+  fromForeignPtr p s = HColumnMajorMatrix $
+                       matrixFromVector ColumnMajor (last s) $
+                       unsafeFromForeignPtr p 0 (Prelude.product s)
+  smap f (HColumnMajorMatrix m) = HColumnMajorMatrix $ C.mapMatrix f m
diff --git a/Data/NetCDF/Metadata.hs b/Data/NetCDF/Metadata.hs
--- a/Data/NetCDF/Metadata.hs
+++ b/Data/NetCDF/Metadata.hs
@@ -1,69 +1,177 @@
-{-# LANGUAGE ExistentialQuantification, StandaloneDeriving #-}
+{-# LANGUAGE FlexibleInstances, UndecidableInstances, EmptyDataDecls #-}
 -- | NetCDF file metadata handling: when a NetCDF file is opened,
 -- metadata defining the dimensions, variables and attributes in the
 -- file are read all at once to create a value of type `NcInfo`.
 
 module Data.NetCDF.Metadata
-       ( NcDim (..)
-       , NcAttr (..)
+       ( Name
+       , NcDim (..)
+       , NcAttr (..), ToNcAttr (..), FromNcAttr (..)
        , NcVar (..)
-       , NcInfo (..)
+       , NcInfo (..), NcRead, NcWrite
        , ncDim, ncAttr, ncVar, ncVarAttr
+       , emptyNcInfo
+       , addNcDim, addNcVar, addNcAttr, addNcVarAttr, (#)
        ) where
 
 import Data.NetCDF.Types
 import qualified Data.Map as M
+import Data.Word
+import Foreign.C
 
+-- | Names for dimensions, variables, attributes.
+type Name = String
 
 -- | Information about a dimension: name, number of entries and
 -- whether unlimited.
-data NcDim = NcDim { ncDimName      :: String
+data NcDim = NcDim { ncDimName      :: Name
                    , ncDimLength    :: Int
                    , ncDimUnlimited :: Bool
                    } deriving Show
 
 -- | Attribute value.
-data NcAttr = forall a. Show a => NcAttr { ncAttrName :: String
-                                         , ncAttrVals :: [a] }
+data NcAttr = NcAttrByte [Word8]
+            | NcAttrChar [Char]
+            | NcAttrShort [CShort]
+            | NcAttrInt [CInt]
+            | NcAttrFloat [CFloat]
+            | NcAttrDouble [CDouble]
+            deriving (Eq, Show)
 
-deriving instance Show NcAttr
+-- | Conversion from Haskell types to attribute values.
+class ToNcAttrHelp a where
+  toAttrHelp :: [a] -> NcAttr
 
+instance ToNcAttrHelp Word8 where toAttrHelp = NcAttrByte
+instance ToNcAttrHelp Char where toAttrHelp = NcAttrChar
+instance ToNcAttrHelp CShort where toAttrHelp = NcAttrShort
+instance ToNcAttrHelp CInt where toAttrHelp = NcAttrInt
+instance ToNcAttrHelp CFloat where toAttrHelp = NcAttrFloat
+instance ToNcAttrHelp CDouble where toAttrHelp = NcAttrDouble
+
+class ToNcAttr a where
+  toAttr :: a -> NcAttr
+
+instance ToNcAttrHelp a => ToNcAttr a where
+  toAttr x = toAttrHelp [x]
+instance ToNcAttrHelp a => ToNcAttr [a] where
+  toAttr xs = toAttrHelp xs
+
+-- | Conversion from attribute values to Haskell types.
+class FromNcAttr a where
+  fromAttr :: NcAttr -> Maybe a
+
+instance FromNcAttr Word8 where
+  fromAttr (NcAttrByte [x]) = Just x
+  fromAttr _ = Nothing
+instance FromNcAttr [Word8] where
+  fromAttr (NcAttrByte xs) = Just xs
+  fromAttr _ = Nothing
+instance FromNcAttr Char where
+  fromAttr (NcAttrChar [x]) = Just x
+  fromAttr _ = Nothing
+instance FromNcAttr [Char] where
+  fromAttr (NcAttrChar xs) = Just xs
+  fromAttr _ = Nothing
+instance FromNcAttr CShort where
+  fromAttr (NcAttrShort [x]) = Just x
+  fromAttr _ = Nothing
+instance FromNcAttr [CShort] where
+  fromAttr (NcAttrShort xs) = Just xs
+  fromAttr _ = Nothing
+instance FromNcAttr CInt where
+  fromAttr (NcAttrInt [x]) = Just x
+  fromAttr _ = Nothing
+instance FromNcAttr [CInt] where
+  fromAttr (NcAttrInt xs) = Just xs
+  fromAttr _ = Nothing
+instance FromNcAttr CFloat where
+  fromAttr (NcAttrFloat [x]) = Just x
+  fromAttr _ = Nothing
+instance FromNcAttr [CFloat] where
+  fromAttr (NcAttrFloat xs) = Just xs
+  fromAttr _ = Nothing
+instance FromNcAttr CDouble where
+  fromAttr (NcAttrDouble [x]) = Just x
+  fromAttr _ = Nothing
+instance FromNcAttr [CDouble] where
+  fromAttr (NcAttrDouble xs) = Just xs
+  fromAttr _ = Nothing
+
+
 -- | Information about a variable: name, type, dimensions and
 -- attributes.
-data NcVar = NcVar { ncVarName  :: String
+data NcVar = NcVar { ncVarName  :: Name
                    , ncVarType  :: NcType
                    , ncVarDims  :: [NcDim]
-                   , ncVarAttrs :: M.Map String NcAttr
+                   , ncVarAttrs :: M.Map Name NcAttr
                    } deriving Show
 
+-- | Type tags for NcInfo values.
+data NcRead
+data NcWrite
+
 -- | Metadata information for a whole NetCDF file.
-data NcInfo = NcInfo { ncName :: FilePath
-                       -- ^ File name.
-                     , ncDims :: M.Map String NcDim
-                       -- ^ Dimensions defined in file.
-                     , ncVars :: M.Map String NcVar
-                       -- ^ Variables defined in file.
-                     , ncAttrs :: M.Map String NcAttr
-                       -- ^ Global attributes defined in file.
-                     , ncId :: NcId
-                       -- ^ Low-level file access ID.
-                     , ncVarIds :: M.Map String NcId
-                       -- ^ Low-level IDs for variables.
-                     } deriving Show
+data NcInfo a = NcInfo { ncName :: FilePath
+                         -- ^ File name.
+                       , ncDims :: M.Map Name NcDim
+                         -- ^ Dimensions defined in file.
+                       , ncVars :: M.Map Name NcVar
+                         -- ^ Variables defined in file.
+                       , ncAttrs :: M.Map Name NcAttr
+                         -- ^ Global attributes defined in file.
+                       , ncId :: NcId
+                         -- ^ Low-level file access ID.
+                       , ncVarIds :: M.Map Name NcId
+                         -- ^ Low-level IDs for variables.
+                       } deriving Show
 
 
 -- | Extract dimension metadata by name.
-ncDim :: NcInfo -> String -> Maybe NcDim
+ncDim :: NcInfo a -> Name -> Maybe NcDim
 ncDim nc n = M.lookup n $ ncDims nc
 
 -- | Extract a global attribute by name.
-ncAttr :: NcInfo -> String -> Maybe NcAttr
+ncAttr :: NcInfo a -> Name -> Maybe NcAttr
 ncAttr nc n = M.lookup n $ ncAttrs nc
 
 -- | Extract variable metadata by name.
-ncVar :: NcInfo -> String -> Maybe NcVar
+ncVar :: NcInfo a -> Name -> Maybe NcVar
 ncVar nc n = M.lookup n $ ncVars nc
 
 -- | Extract an attribute for a given variable by name.
-ncVarAttr :: NcVar -> String -> Maybe NcAttr
+ncVarAttr :: NcVar -> Name -> Maybe NcAttr
 ncVarAttr v n = M.lookup n $ ncVarAttrs v
+
+
+-- | Empty NcInfo value to build on.
+emptyNcInfo :: FilePath -> NcInfo NcWrite
+emptyNcInfo n = NcInfo n M.empty M.empty M.empty ncInvalidId M.empty
+
+-- | Add a new dimension to an NcInfo value.
+addNcDim :: NcDim -> NcInfo NcWrite -> NcInfo NcWrite
+addNcDim dim@(NcDim name _ _) (NcInfo n ds vs as fid vids) =
+  NcInfo n (M.insert name chkdim ds) vs as fid vids
+  where chkdim = dim { ncDimUnlimited =
+                          ncDimUnlimited dim &&
+                          (not . or . (map ncDimUnlimited) . M.elems $ ds) }
+
+-- | Add a new variable to an NcInfo value.
+addNcVar :: NcVar -> NcInfo NcWrite -> NcInfo NcWrite
+addNcVar var@(NcVar name _ _ _) (NcInfo n ds vs as fid vids) =
+  NcInfo n ds (M.insert name var vs) as fid vids
+
+-- | Add a new global attribute to an NcInfo value.
+addNcAttr :: Name -> NcAttr -> NcInfo NcWrite -> NcInfo NcWrite
+addNcAttr name att (NcInfo n ds vs as fid vids) =
+  NcInfo n ds vs (M.insert name att as) fid vids
+
+-- | Add a new attribute to an NcVar value.
+addNcVarAttr :: Name -> NcAttr -> NcVar -> NcVar
+addNcVarAttr name att (NcVar n t ds as) = NcVar n t ds (M.insert name att as)
+
+infixl 8 #
+
+-- | Handy postfix function application.
+(#) :: a -> (a -> b) -> b
+(#) = flip ($)
diff --git a/Data/NetCDF/PutGet.hs b/Data/NetCDF/PutGet.hs
--- a/Data/NetCDF/PutGet.hs
+++ b/Data/NetCDF/PutGet.hs
@@ -14,12 +14,9 @@
 import Foreign.Marshal.Alloc
 import Control.Applicative ((<$>))
 import Control.Monad (liftM)
-import Data.Word
 
-import Data.NetCDF.Types
 import Data.NetCDF.Storable
 import Data.NetCDF.Store
-import Data.NetCDF.Raw.PutGet
 
 put_var1 :: NcStorable a => Int -> Int -> [Int] -> a -> IO Int
 put_var1 nc var idxs v = do
diff --git a/Data/NetCDF/Raw.hs b/Data/NetCDF/Raw.hs
--- a/Data/NetCDF/Raw.hs
+++ b/Data/NetCDF/Raw.hs
@@ -4,10 +4,10 @@
        ( module Data.NetCDF.Raw.Base
        , module Data.NetCDF.Raw.PutGet
        , module Data.NetCDF.Raw.Attributes
-       , module Data.NetCDF.Raw.NetCDF4
+--       , module Data.NetCDF.Raw.NetCDF4
        ) where
 
 import Data.NetCDF.Raw.Base
 import Data.NetCDF.Raw.PutGet
 import Data.NetCDF.Raw.Attributes
-import Data.NetCDF.Raw.NetCDF4
+--import Data.NetCDF.Raw.NetCDF4
diff --git a/Data/NetCDF/Raw/NetCDF4.chs b/Data/NetCDF/Raw/NetCDF4.chs
deleted file mode 100644
--- a/Data/NetCDF/Raw/NetCDF4.chs
+++ /dev/null
@@ -1,151 +0,0 @@
-{-# LANGUAGE ForeignFunctionInterface #-}
-
--- | Raw bindings for NetCDF4 features (not implemented).
-
-module Data.NetCDF.Raw.NetCDF4 where
-
-#include <netcdf.h>
-
---------------------------------------------------------------------------------
---
---  UNSUPPORTED FEATURES (NetCDF-4)
---
---------------------------------------------------------------------------------
-
--- USER DEFINED TYPES
---
--- int nc_def_compound(int ncid, size_t size,
---                     const char *name, nc_type *typeidp);
--- int nc_insert_compound(int ncid, nc_type xtype, const char *name,
---                        size_t offset, nc_type field_typeid);
--- int nc_insert_array_compound(int ncid, nc_type xtype, const char *name,
---                              size_t offset, nc_type field_typeid,
---                              int ndims, const int *dim_sizes);
--- int nc_inq_type(int ncid, nc_type xtype, char *name, size_t *size);
--- int nc_inq_compound(int ncid, nc_type xtype, char *name, size_t *sizep,
---                     size_t *nfieldsp);
--- int nc_inq_compound_name(int ncid, nc_type xtype, char *name);
--- int nc_inq_compound_size(int ncid, nc_type xtype, size_t *sizep);
--- int nc_inq_compound_nfields(int ncid, nc_type xtype, size_t *nfieldsp);
--- int nc_inq_compound_fieldname(int ncid, nc_type xtype, int fieldid,
---                               char *name);
--- int nc_inq_compound_fieldindex(int ncid, nc_type xtype, const char *name,
---                                int *fieldidp);
--- int nc_inq_compound_fieldoffset(int ncid, nc_type xtype, int fieldid,
---                                 size_t *offsetp);
--- int nc_inq_compound_fieldtype(int ncid, nc_type xtype, int fieldid,
---                               nc_type *field_typeidp);
--- int nc_inq_compound_fieldndims(int ncid, nc_type xtype, int fieldid,
---                                int *ndimsp);
--- int nc_inq_compound_fielddim_sizes(int ncid, nc_type xtype, int fieldid,
---                                    int *dim_sizes);
--- typedef struct {
---     size_t len; /**< Length of VL data (in base type units) */
---     void *p;    /**< Pointer to VL data */
--- } nc_vlen_t;
--- int nc_def_vlen(int ncid, const char *name,
---                 nc_type base_typeid, nc_type *xtypep);
--- int nc_inq_vlen(int ncid, nc_type xtype, char *name, size_t *datum_sizep,
---                 nc_type *base_nc_typep);
--- int nc_free_vlen(nc_vlen_t *vl);
--- int nc_put_vlen_element(int ncid, int typeid1, void *vlen_element,
---                         size_t len, const void *data);
--- int nc_get_vlen_element(int ncid, int typeid1, const void *vlen_element,
---                         size_t *len, void *data);
--- int nc_free_string(size_t len, char **data);
--- int nc_inq_user_type(int ncid, nc_type xtype, char *name, size_t *size,
---                      nc_type *base_nc_typep, size_t *nfieldsp, int *classp);
--- int nc_def_enum(int ncid, nc_type base_typeid, const char *name,
---                 nc_type *typeidp);
--- int nc_insert_enum(int ncid, nc_type xtype, const char *name,
---                    const void *value);
--- int nc_inq_enum_member(int ncid, nc_type xtype, int idx, char *name,
---                        void *value);
--- int nc_inq_enum_ident(int ncid, nc_type xtype, long long value,
---                       char *identifier);
--- int nc_def_opaque(int ncid, size_t size, const char *name, nc_type *xtypep);
--- int nc_inq_opaque(int ncid, nc_type xtype, char *name, size_t *sizep);
--- int nc_inq_typeid(int ncid, const char *name, nc_type *typeidp);
--- int nc_inq_compound_field(int ncid, nc_type xtype, int fieldid, char *name,
---                           size_t *offsetp, nc_type *field_typeidp,
---                           int *ndimsp, int *dim_sizesp);
--- int nc_free_vlens(size_t len, nc_vlen_t vlens[]);
--- int nc_inq_enum(int ncid, nc_type xtype, char *name, nc_type *base_nc_typep,
---                 size_t *base_sizep, size_t *num_membersp);
-
--- GROUPS
---
--- int nc_inq_ncid(int ncid, const char *name, int *grp_ncid);
--- int nc_inq_grps(int ncid, int *numgrps, int *ncids);
--- int nc_inq_grpname(int ncid, char *name);
--- int nc_inq_grpname_full(int ncid, size_t *lenp, char *full_name);
--- int nc_inq_grpname_len(int ncid, size_t *lenp);
--- int nc_inq_grp_parent(int ncid, int *parent_ncid);
--- int nc_inq_grp_ncid(int ncid, const char *grp_name, int *grp_ncid);
--- int nc_inq_grp_full_ncid(int ncid, const char *full_name, int *grp_ncid);
--- int nc_inq_varids(int ncid, int *nvars, int *varids);
--- int nc_inq_dimids(int ncid, int *ndims, int *dimids, int include_parents);
--- int nc_inq_typeids(int ncid, int *ntypes, int *typeids);
--- int nc_def_grp(int parent_ncid, const char *name, int *new_ncid);
-
--- NETCDF-4 VARIABLES
---
--- int nc_def_var_deflate(int ncid, int varid, int shuffle, int deflate,
---                        int deflate_level);
--- int nc_inq_var_deflate(int ncid, int varid, int *shufflep,
---                        int *deflatep, int *deflate_levelp);
--- int nc_def_var_fletcher32(int ncid, int varid, int fletcher32);
--- int nc_inq_var_fletcher32(int ncid, int varid, int *fletcher32p);
--- int nc_def_var_chunking(int ncid, int varid, int storage,
---                         const size_t *chunksizesp);
--- int nc_inq_var_chunking(int ncid, int varid, int *storagep,
---                         size_t *chunksizesp);
--- int nc_def_var_fill(int ncid, int varid, int no_fill,
---                     const void *fill_value);
--- int nc_inq_var_fill(int ncid, int varid, int *no_fill, void *fill_valuep);
--- int nc_def_var_endian(int ncid, int varid, int endian);
--- int nc_inq_var_endian(int ncid, int varid, int *endianp);
--- int nc_set_chunk_cache(size_t size, size_t nelems, float preemption);
--- int nc_get_chunk_cache(size_t *sizep, size_t *nelemsp, float *preemptionp);
--- int nc_set_var_chunk_cache(int ncid, int varid, size_t size, size_t nelems,
---                            float preemption);
--- int nc_get_var_chunk_cache(int ncid, int varid, size_t *sizep,
---                            size_t *nelemsp, float *preemptionp);
--- int nc_inq_var_szip(int ncid, int varid, int *options_maskp,
---                     int *pixels_per_blockp);
-
--- Variable length strings are only supported in NetCDF-4 files, which
--- we don't handle yet.
---
--- int nc_put_var_string(int ncid, int varid, const char **op);
--- int nc_get_var_string(int ncid, int varid, char **ip);
--- int nc_put_var1_string(int ncid, int varid, const size_t *indexp,
---                        const char **op);
--- int nc_get_var1_string(int ncid, int varid, const size_t *indexp,
---                        char **ip);
--- int nc_put_vara_string(int ncid, int varid, const size_t *startp,
---                        const size_t *countp, const char **op);
--- int nc_get_vara_string(int ncid, int varid, const size_t *startp,
---                        const size_t *countp, char **ip);
--- int nc_put_vars_string(int ncid, int varid, const size_t *startp,
---                        const size_t *countp, const ptrdiff_t *stridep,
---                        const char **op);
--- int nc_get_vars_string(int ncid, int varid, const size_t *startp,
---                        const size_t *countp, const ptrdiff_t *stridep,
---                        char **ip);
--- int nc_put_varm_string(int ncid, int varid, const size_t *startp,
---                        const size_t *countp, const ptrdiff_t *stridep,
---                        const ptrdiff_t * imapp, const char **op);
--- int nc_get_varm_string(int ncid, int varid, const size_t *startp,
---                        const size_t *countp, const ptrdiff_t *stridep,
---                        const ptrdiff_t * imapp, char **ip);
--- int nc_put_att_string(int ncid, int varid, const char *name,
---                       size_t len, const char **op);
--- int nc_get_att_string(int ncid, int varid, const char *name, char **ip);
-
--- MISCELLANEOUS
---
--- int nc_var_par_access(int ncid, int varid, int par_access);
--- int nc_inq_type_equal(int ncid1, nc_type typeid1, int ncid2,
---                       nc_type typeid2, int *equal);
--- int nc_set_default_format(int format, int *old_formatp);
diff --git a/Data/NetCDF/Repa.hs b/Data/NetCDF/Repa.hs
new file mode 100644
--- /dev/null
+++ b/Data/NetCDF/Repa.hs
@@ -0,0 +1,16 @@
+{-# LANGUAGE FlexibleInstances #-}
+{-# OPTIONS_GHC -fno-warn-orphans #-}
+-- | NetCDF store instance for Repa arrays.
+module Data.NetCDF.Repa where
+
+import Data.NetCDF.Store
+
+import Data.Array.Repa
+import qualified Data.Array.Repa as R
+import qualified Data.Array.Repa.Repr.ForeignPtr as RF
+import Data.Array.Repa.Repr.ForeignPtr (F)
+
+instance Shape sh => NcStore (Array F sh) where
+  toForeignPtr = RF.toForeignPtr
+  fromForeignPtr p s = RF.fromForeignPtr (shapeOfList $ reverse s) p
+  smap f s = computeS $ R.map f s
diff --git a/Data/NetCDF/Storable.hs b/Data/NetCDF/Storable.hs
--- a/Data/NetCDF/Storable.hs
+++ b/Data/NetCDF/Storable.hs
@@ -11,7 +11,6 @@
 import Foreign.Storable
 
 import Data.NetCDF.Types
-import Data.NetCDF.Store
 import Data.NetCDF.Raw.PutGet
 
 
@@ -102,52 +101,4 @@
   ffi_get_vara = nc_get_vara_double'_
   ffi_put_vars = nc_put_vars_double'_
   ffi_get_vars = nc_get_vars_double'_
-
-instance NcStorable CUChar where
-  -- #define NC_UBYTE 7 /**< unsigned 1 byte int */
-  ncType _ = NcUByte
-  ffi_put_var1 = nc_put_var1_uchar'_
-  ffi_get_var1 = nc_get_var1_uchar'_
-  ffi_put_var = nc_put_var_uchar'_
-  ffi_get_var = nc_get_var_uchar'_
-  ffi_put_vara = nc_put_vara_uchar'_
-  ffi_get_vara = nc_get_vara_uchar'_
-  ffi_put_vars = nc_put_vars_uchar'_
-  ffi_get_vars = nc_get_vars_uchar'_
-
-instance NcStorable CUShort where
-  -- #define NC_USHORT 8 /**< unsigned 2-byte int */
-  ncType _ = NcUShort
-  ffi_put_var1 = nc_put_var1_ushort'_
-  ffi_get_var1 = nc_get_var1_ushort'_
-  ffi_put_var = nc_put_var_ushort'_
-  ffi_get_var = nc_get_var_ushort'_
-  ffi_put_vara = nc_put_vara_ushort'_
-  ffi_get_vara = nc_get_vara_ushort'_
-  ffi_put_vars = nc_put_vars_ushort'_
-  ffi_get_vars = nc_get_vars_ushort'_
-
-instance NcStorable CUInt where
-  -- #define NC_UINT 9 /**< unsigned 4-byte int */
-  ncType _ = NcUInt
-  ffi_put_var1 = nc_put_var1_uint'_
-  ffi_get_var1 = nc_get_var1_uint'_
-  ffi_put_var = nc_put_var_uint'_
-  ffi_get_var = nc_get_var_uint'_
-  ffi_put_vara = nc_put_vara_uint'_
-  ffi_get_vara = nc_get_vara_uint'_
-  ffi_put_vars = nc_put_vars_uint'_
-  ffi_get_vars = nc_get_vars_uint'_
-
-instance NcStorable CLong where
-  -- #define NC_INT64 10 /**< signed 8-byte int */
-  ncType _ = NcInt64
-  ffi_put_var1 = nc_put_var1_long'_
-  ffi_get_var1 = nc_get_var1_long'_
-  ffi_put_var = nc_put_var_long'_
-  ffi_get_var = nc_get_var_long'_
-  ffi_put_vara = nc_put_vara_long'_
-  ffi_get_vara = nc_get_vara_long'_
-  ffi_put_vars = nc_put_vars_long'_
-  ffi_get_vars = nc_get_vars_long'_
 
diff --git a/Data/NetCDF/Store.hs b/Data/NetCDF/Store.hs
--- a/Data/NetCDF/Store.hs
+++ b/Data/NetCDF/Store.hs
@@ -5,19 +5,10 @@
 
 module Data.NetCDF.Store where
 
-import Data.List (reverse)
 import Foreign.Storable
 import Foreign.ForeignPtr
 
-import Data.Vector.Storable hiding (reverse)
 
-import Data.Array.Repa
-import qualified Data.Array.Repa as R
-import qualified Data.Array.Repa.Repr.ForeignPtr as RF
-import Data.Array.Repa.Repr.ForeignPtr (F)
-
-import Data.NetCDF.Types
-
 -- | 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
@@ -26,13 +17,3 @@
   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
-
-instance NcStore Vector where
-  toForeignPtr = fst . unsafeToForeignPtr0
-  fromForeignPtr p s = unsafeFromForeignPtr0 p (Prelude.product s)
-  smap = Data.Vector.Storable.map
-
-instance Shape sh => NcStore (Array F sh) where
-  toForeignPtr = RF.toForeignPtr
-  fromForeignPtr p s = RF.fromForeignPtr (shapeOfList $ reverse s) p
-  smap f s = computeS $ R.map f s
diff --git a/Data/NetCDF/Types.hs b/Data/NetCDF/Types.hs
--- a/Data/NetCDF/Types.hs
+++ b/Data/NetCDF/Types.hs
@@ -5,7 +5,6 @@
 
 import Control.Exception
 import Data.Typeable
-import System.FilePath
 import System.IO (IOMode (..))
 
 -- | Representation for NetCDF external data types.
@@ -16,12 +15,6 @@
             | NcInt     -- ^ Signed 4 byte integer
             | NcFloat   -- ^ Single precision floating point
             | NcDouble  -- ^ Double precision floating point
-            | NcUByte   -- ^ Unsigned 1 byte int
-            | NcUShort  -- ^ Unsigned 2-byte int
-            | NcUInt    -- ^ Unsigned 4-byte int
-            | NcInt64   -- ^ Signed 8-byte int
-            | NcUInt64  -- ^ Unsigned 8-byte int
-            | NcString  -- ^ String
             deriving (Eq, Show)
 
 instance Enum NcType where
@@ -31,12 +24,6 @@
   fromEnum NcInt = 4
   fromEnum NcFloat = 5
   fromEnum NcDouble = 6
-  fromEnum NcUByte = 7
-  fromEnum NcUShort = 8
-  fromEnum NcUInt = 9
-  fromEnum NcInt64 = 10
-  fromEnum NcUInt64 = 11
-  fromEnum NcString = 12
   toEnum n = case n of
     1  -> NcByte
     2  -> NcChar
@@ -44,12 +31,6 @@
     4  -> NcInt
     5  -> NcFloat
     6  -> NcDouble
-    7  -> NcUByte
-    8  -> NcUShort
-    9  -> NcUInt
-    10 -> NcInt64
-    11 -> NcUInt64
-    12 -> NcString
     _ -> throw (NcInvalidType n)
 
 -- | Internal representation of NetCDF IDs.
@@ -73,3 +54,11 @@
 -- | Fake variable identifier for global attributes.
 ncGlobal :: Int
 ncGlobal = -1
+
+-- | Fake length value for defining unlimited dimensions.
+ncUnlimitedLength :: Int
+ncUnlimitedLength = 0
+
+-- | Fake file identifier for unopened files.
+ncInvalidId :: NcId
+ncInvalidId = -1
diff --git a/Data/NetCDF/Utils.hs b/Data/NetCDF/Utils.hs
--- a/Data/NetCDF/Utils.hs
+++ b/Data/NetCDF/Utils.hs
@@ -26,6 +26,11 @@
   status :: a -> Int
   proj :: a -> OutType a
 
+instance Checkable Int where
+  type OutType Int = ()
+  status s = s
+  proj _ = ()
+
 instance Checkable (Int, a) where
   type OutType (Int, a) = a
   status (s, _) = s
diff --git a/Data/NetCDF/Vector.hs b/Data/NetCDF/Vector.hs
new file mode 100644
--- /dev/null
+++ b/Data/NetCDF/Vector.hs
@@ -0,0 +1,13 @@
+{-# LANGUAGE FlexibleInstances #-}
+{-# OPTIONS_GHC -fno-warn-orphans #-}
+-- | NetCDF store instance for Storable Vectors.
+module Data.NetCDF.Vector where
+
+import Data.NetCDF.Store
+
+import qualified Data.Vector.Storable as SV
+
+instance NcStore SV.Vector where
+  toForeignPtr = fst . SV.unsafeToForeignPtr0
+  fromForeignPtr p s = SV.unsafeFromForeignPtr0 p (Prelude.product s)
+  smap = SV.map
diff --git a/README.markdown b/README.markdown
new file mode 100644
--- /dev/null
+++ b/README.markdown
@@ -0,0 +1,46 @@
+hnetcdf
+=======
+
+Haskell NetCDF library: as well as conventional low-level FFI bindings
+to the functions in the NetCDF library (in the `Data.NetCDF.Raw`
+modules), `hnetcdf` provides a higher-level Haskell interface
+(currently only for reading data).  This higher-level interface aims
+to provide a "container polymorphic" view of NetCDF data allowing
+NetCDF variables to be read into `Storable` `Vectors` and Repa arrays
+easily.
+
+For example:
+
+``` haskell
+import Data.NetCDF
+import Foreign.C
+import qualified Data.Vector.Storable as SV
+...
+type SVRet = IO (Either NcError (SV.Vector a))
+...
+  enc <- openFile "tst.nc" ReadMode
+  case enc of
+    Right nc -> do
+      eval <- get nc "varname" :: SVRet CDouble
+      ...
+```
+
+gets the full contents of a NetCDF variable as a `Storable` `Vector`,
+while the following code reads the same variable (assumed to be
+three-dimensional) into a Repa array:
+
+``` haskell
+import Data.NetCDF
+import Foreign.C
+import qualified Data.Array.Repa as R
+import qualified Data.Array.Repa.Eval as RE
+import Data.Array.Repa.Repr.ForeignPtr (F)
+...
+type RepaRet3 a = IO (Either NcError (R.Array F R.DIM3 a))
+...
+  enc <- openFile "tst.nc" ReadMode
+  case enc of
+    Right nc -> do
+      eval <- get nc "varname" :: RepaRet3 CDouble
+      ...
+```
diff --git a/hnetcdf.cabal b/hnetcdf.cabal
--- a/hnetcdf.cabal
+++ b/hnetcdf.cabal
@@ -1,5 +1,5 @@
 name:                hnetcdf
-version:             0.1.0.0
+version:             0.2.0.0
 synopsis:            Haskell NetCDF library
 description:
   Bindings to the Unidata NetCDF library, along with a higher-level
@@ -12,6 +12,7 @@
 copyright:           Copyright (2013) Ian Ross
 category:            Data
 build-type:          Simple
+extra-source-files:  README.markdown CHANGELOG.markdown
 cabal-version:       >=1.8
 homepage:            https://github.com/ian-ross/hnetcdf
 
@@ -21,6 +22,9 @@
 
 Library
   exposed-modules:     Data.NetCDF
+                       Data.NetCDF.Vector
+                       Data.NetCDF.Repa
+                       Data.NetCDF.HMatrix
                        Data.NetCDF.Metadata
                        Data.NetCDF.PutGet
                        Data.NetCDF.Storable
@@ -31,18 +35,20 @@
                        Data.NetCDF.Raw.Base
                        Data.NetCDF.Raw.PutGet
                        Data.NetCDF.Raw.Attributes
-                       Data.NetCDF.Raw.NetCDF4
                        Data.NetCDF.Raw.Utils
+--                       Data.NetCDF.Raw.NetCDF4
   build-depends:       base                        >= 4.3      && < 5,
                        containers                  >= 0.4      && < 0.6,
-                       filepath                    >= 1.3      && < 1.4,
-                       transformers                >= 0.2      && < 0.4,
                        either                      >= 3.4.1    && < 3.5,
                        errors                      >= 1.4.2    && < 1.5,
-                       vector                      >= 0.9      && < 0.11,
-                       repa                        >= 3.2      && < 3.3
+                       filepath                    >= 1.3      && < 1.4,
+                       hmatrix                     >= 0.15.2.0 && < 0.16,
+                       repa,
+                       transformers                >= 0.2      && < 0.5,
+                       vector                      >= 0.9      && < 0.11
   build-tools:         c2hs
   extra-libraries:     netcdf
+  ghc-options:         -Wall
 
 Test-Suite test-raw-metadata
   type:                exitcode-stdio-1.0
@@ -52,7 +58,7 @@
                        base                        >= 4.3      && < 5,
                        containers                  >= 0.4      && < 0.6,
                        vector                      >= 0.9      && < 0.11,
-                       repa                        >= 3.2      && < 3.3,
+                       repa,
                        test-framework              >= 0.3,
                        test-framework-hunit,
                        test-framework-quickcheck2  >= 0.3,
@@ -68,7 +74,7 @@
                        base                        >= 4.3      && < 5,
                        containers                  >= 0.4      && < 0.6,
                        vector                      >= 0.9      && < 0.11,
-                       repa                        >= 3.2      && < 3.3,
+                       repa,
                        directory                   >= 1.1      && < 1.3,
                        test-framework              >= 0.3,
                        test-framework-hunit,
@@ -85,7 +91,7 @@
                        base                        >= 4.3      && < 5,
                        containers                  >= 0.4      && < 0.6,
                        vector                      >= 0.9      && < 0.11,
-                       repa                        >= 3.2      && < 3.3,
+                       repa,
                        directory                   >= 1.1      && < 1.3,
                        test-framework              >= 0.3,
                        test-framework-hunit,
@@ -102,8 +108,9 @@
                        base                        >= 4.3      && < 5,
                        containers                  >= 0.4      && < 0.6,
                        errors                      >= 1.4.2    && < 1.5,
+                       hmatrix                     >= 0.15.2.0 && < 0.16,
+                       repa,
                        vector                      >= 0.9      && < 0.11,
-                       repa                        >= 3.2      && < 3.3,
                        directory                   >= 1.1      && < 1.3,
                        test-framework              >= 0.3,
                        test-framework-hunit,
@@ -112,3 +119,14 @@
                        QuickCheck
   extra-libraries:     netcdf
 
+Test-Suite test-put
+  type:                exitcode-stdio-1.0
+  hs-source-dirs:      test
+  main-is:             test-put.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
+  extra-libraries:     netcdf
diff --git a/test/test-get.hs b/test/test-get.hs
--- a/test/test-get.hs
+++ b/test/test-get.hs
@@ -8,9 +8,13 @@
 import Data.Array.Repa.Repr.ForeignPtr (F)
 import Control.Error
 import Control.Monad
-
+import Numeric.Container as C
 import Foreign.C
+
 import Data.NetCDF
+import qualified Data.NetCDF.Vector as V
+import qualified Data.NetCDF.Repa as R
+import qualified Data.NetCDF.HMatrix as H
 
 main :: IO ()
 main = defaultMain tests
@@ -41,6 +45,9 @@
     , testCase "Read (float, Repa)" (getVarRepa infile "vf" fwit)
     , testCase "Read (int, Repa)" (getVarRepa infile "vi" iwit)
     , testCase "Read (short, Repa)" (getVarRepa infile "vs" swit)
+    , testCase "Read (float, hmatrix)" (getVarHMV infile "vf" fwit)
+    , testCase "Read (int, hmatrix)" (getVarHMV infile "vi" iwit)
+    , testCase "Read (short, hmatrix)" (getVarHMV infile "vs" swit)
     ],
     testGroup "Variable slice access functions"
     [ testCase "Read (float, SV)" (getVarASV infile "vf" fwit)
@@ -49,6 +56,9 @@
     , testCase "Read (float, Repa)" (getVarARepa infile "vf" fwit)
     , testCase "Read (int, Repa)" (getVarARepa infile "vi" iwit)
     , testCase "Read (short, Repa)" (getVarARepa infile "vs" swit)
+    , testCase "Read (float, hmatrix)" (getVarAHMV infile "vf" fwit)
+    , testCase "Read (int, hmatrix)" (getVarAHMV infile "vi" iwit)
+    , testCase "Read (short, hmatrix)" (getVarAHMV infile "vs" swit)
     ],
     testGroup "Strided slice access functions"
     [ testCase "Read (float, SV)" (getVarSSV infile "vf" fwit)
@@ -57,6 +67,9 @@
     , testCase "Read (float, Repa)" (getVarSRepa infile "vf" fwit)
     , testCase "Read (int, Repa)" (getVarSRepa infile "vi" iwit)
     , testCase "Read (short, Repa)" (getVarSRepa infile "vs" swit)
+    , testCase "Read (float, hmatrix)" (getVarSHMV infile "vf" fwit)
+    , testCase "Read (int, hmatrix)" (getVarSHMV infile "vi" iwit)
+    , testCase "Read (short, hmatrix)" (getVarSHMV infile "vs" swit)
     ]
   ]
 
@@ -64,6 +77,8 @@
 type SVRet a = IO (Either NcError (SV.Vector a))
 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))
 
 
 --------------------------------------------------------------------------------
@@ -75,7 +90,7 @@
 getVar1 :: forall a. (Eq a, Num a, Show a, NcStorable a)
         => FilePath -> String -> a -> Assertion
 getVar1 f v _ = do
-  enc <- openFile f ReadMode
+  enc <- openFile f
   assertBool "failed to open file" $ isRight enc
   let Right nc = enc
   case ncVar nc v of
@@ -95,31 +110,7 @@
                   (val == trueval)
   void $ closeFile nc
 
--- getVar1 :: FilePath -> String -> Assertion
--- getVar1 f v = do
---   enc <- openFile f ReadMode
---   assertBool "failed to open file" $ isRight enc
---   let Right nc = enc
---   case ncVar nc v of
---     Nothing -> assertBool "missing variable" False
---     Just var -> do
---       let [nz, ny, nx] = map ncDimLength $ ncVarDims var
---           typ = ncVarType var
---       forM_ [1..nx] $ \ix -> do
---         forM_ [1..ny] $ \iy -> do
---           forM_ [1..nz] $ \iz -> do
---             eval <- get1 nc var [iz-1, iy-1, ix-1] :: IO (Either NcError a)
---             case eval of
---               Left e -> assertBool ("read error: " ++ show e) False
---               Right val -> do
---                 let trueval = fromIntegral $ ix + 10 * iy + 100 * iz
---                 assertBool ("value error: " ++ show val ++
---                             " instead of " ++ show trueval)
---                   (val == trueval)
---   void $ closeFile nc
 
-
-
 --------------------------------------------------------------------------------
 --
 --  WHOLE VARIABLE READ
@@ -129,7 +120,7 @@
 getVarSV :: forall a. (Eq a, Num a, Show a, NcStorable a)
          => FilePath -> String -> a -> Assertion
 getVarSV f v _ = do
-  enc <- openFile f ReadMode
+  enc <- openFile f
   assertBool "failed to open file" $ isRight enc
   let Right nc = enc
   case ncVar nc v of
@@ -153,7 +144,7 @@
 getVarRepa :: forall a. (Eq a, Num a, Show a, NcStorable a)
            => FilePath -> String -> a -> Assertion
 getVarRepa f v _ = do
-  enc <- openFile f ReadMode
+  enc <- openFile f
   assertBool "failed to open file" $ isRight enc
   let Right nc = enc
   case ncVar nc v of
@@ -174,8 +165,32 @@
             (and $ zipWith (==) tstvals truevals)
   void $ closeFile nc
 
+getVarHMV :: forall a. (Eq a, Num a, Show a, NcStorable a)
+          => FilePath -> String -> a -> Assertion
+getVarHMV f v _ = do
+  enc <- openFile f
+  assertBool "failed to open file" $ isRight enc
+  let Right nc = enc
+  case ncVar nc v of
+    Nothing -> assertBool "missing variable" False
+    Just var -> do
+      eval <- get nc var :: HMVRet a
+      case eval of
+        Left e -> assertBool ("read error: " ++ show e) False
+        Right (H.HVector val) -> do
+          let [nz, ny, nx] = map ncDimLength $ ncVarDims var
+              truevals = [fromIntegral $ ix + 10 * iy + 100 * iz |
+                          ix <- [1..nx], iy <- [1..ny], iz <- [1..nz]]
+              idxs = [(iz - 1) * ny * nx + (iy - 1) * nx + (ix - 1) |
+                      ix <- [1..nx], iy <- [1..ny], iz <- [1..nz]]
+              tstvals = [val C.@> i | i <- idxs]
+          assertBool ("value error: " ++ show tstvals ++
+                      " instead of " ++ show truevals)
+            (and $ zipWith (==) tstvals truevals)
+  void $ closeFile nc
 
 
+
 --------------------------------------------------------------------------------
 --
 --  VARIABLE SLICE READ
@@ -185,7 +200,7 @@
 getVarASV :: forall a. (Num a, Show a, Eq a, NcStorable a)
           => FilePath -> String -> a -> Assertion
 getVarASV f v _ = do
-  enc <- openFile f ReadMode
+  enc <- openFile f
   assertBool "failed to open file" $ isRight enc
   let Right nc = enc
   case ncVar nc v of
@@ -232,7 +247,7 @@
 getVarARepa :: forall a. (Num a, Show a, Eq a, NcStorable a)
             => FilePath -> String -> a -> Assertion
 getVarARepa f v _ = do
-  enc <- openFile f ReadMode
+  enc <- openFile f
   assertBool "failed to open file" $ isRight enc
   let Right nc = enc
   case ncVar nc v of
@@ -276,7 +291,54 @@
               assertBool ("3: value error: " ++ show tstvals ++
                           " instead of " ++ show truevals) (tstvals == truevals)
 
+getVarAHMV :: forall a. (Num a, Show a, Eq a, NcStorable a)
+           => FilePath -> String -> a -> Assertion
+getVarAHMV f v _ = do
+  enc <- openFile f
+  assertBool "failed to open file" $ isRight enc
+  let Right nc = enc
+  case ncVar nc v of
+    Nothing -> assertBool "missing variable" False
+    Just var -> do
+      let [nz, ny, nx] = map ncDimLength $ ncVarDims var
+      forM_ [1..nx] $ \ix -> do
+        forM_ [1..ny] $ \iy -> do
+          let start = [0, iy-1, ix-1]
+              count = [nz, 1, 1]
+          eval <- getA nc var start count :: HMVRet a
+          case eval of
+            Left e -> assertBool ("read error: " ++ show e) False
+            Right vals -> do
+              let truevals = H.HVector $ C.buildVector (fromIntegral nz)
+                             (\iz -> fromIntegral $ ix + 10 * iy + 100 * (iz+1))
+              assertBool ("1: value error: " ++ show vals ++
+                          " instead of " ++ show truevals) (vals == truevals)
+      forM_ [1..nx] $ \ix -> do
+        forM_ [1..nz] $ \iz -> do
+          let start = [iz-1, 0, ix-1]
+              count = [1, ny, 1]
+          eval <- getA nc var start count :: HMVRet a
+          case eval of
+            Left e -> assertBool ("read error: " ++ show e) False
+            Right vals -> do
+              let truevals = H.HVector $ C.buildVector (fromIntegral ny)
+                             (\iy -> fromIntegral $ ix + 10 * (iy+1) + 100 * iz)
+              assertBool ("2: value error: " ++ show vals ++
+                          " instead of " ++ show truevals) (vals == truevals)
+      forM_ [1..ny] $ \iy -> do
+        forM_ [1..nz] $ \iz -> do
+          let start = [iz-1, iy-1, 0]
+              count = [1, 1, nx]
+          eval <- getA nc var start count :: HMVRet a
+          case eval of
+            Left e -> assertBool ("read error: " ++ show e) False
+            Right vals -> do
+              let truevals = H.HVector $ C.buildVector (fromIntegral nx)
+                             (\ix -> fromIntegral $ (ix+1) + 10 * iy + 100 * iz)
+              assertBool ("3: value error: " ++ show vals ++
+                          " instead of " ++ show truevals) (vals == truevals)
 
+
 --------------------------------------------------------------------------------
 --
 --  STRIDED SLICE READ
@@ -286,7 +348,7 @@
 getVarSSV :: forall a. (Num a, Show a, Eq a, NcStorable a)
           => FilePath -> String -> a -> Assertion
 getVarSSV f v _ = do
-  enc <- openFile f ReadMode
+  enc <- openFile f
   assertBool "failed to open file" $ isRight enc
   let Right nc = enc
   case ncVar nc v of
@@ -342,7 +404,7 @@
 getVarSRepa :: forall a. (Num a, Show a, Eq a, NcStorable a)
             => FilePath -> String -> a -> Assertion
 getVarSRepa f v _ = do
-  enc <- openFile f ReadMode
+  enc <- openFile f
   assertBool "failed to open file" $ isRight enc
   let Right nc = enc
   case ncVar nc v of
@@ -403,3 +465,62 @@
                 assertBool ("3: value error: " ++ show tstvals ++
                             " instead of " ++ show truevals)
                   (tstvals == truevals)
+
+getVarSHMV :: forall a. (Num a, Show a, Eq a, NcStorable a)
+          => FilePath -> String -> a -> Assertion
+getVarSHMV f v _ = do
+  enc <- openFile f
+  assertBool "failed to open file" $ isRight enc
+  let Right nc = enc
+  case ncVar nc v of
+    Nothing -> assertBool "missing variable" False
+    Just var -> do
+      let [nz, ny, nx] = map ncDimLength $ ncVarDims var
+      forM_ [1..nx] $ \ix -> do
+        forM_ [1..ny] $ \iy -> do
+          forM_ [1..nz] $ \s -> do
+            let start = [0, iy-1, ix-1]
+                count = [(nz + s - 1) `div` s, 1, 1]
+                stride = [s, 1, 1]
+            eval <- getS nc var start count stride :: HMVRet a
+            case eval of
+              Left e -> assertBool ("read error: " ++ show e) False
+              Right vals -> do
+                let truevals = H.HVector $ C.buildVector
+                               (fromIntegral $ (nz + s - 1) `div` s)
+                               (\iz -> fromIntegral $
+                                       ix + 10 * iy + 100 * (iz * s + 1))
+                assertBool ("1: value error: " ++ show vals ++
+                            " instead of " ++ show truevals) (vals == truevals)
+      forM_ [1..nx] $ \ix -> do
+        forM_ [1..nz] $ \iz -> do
+          forM_ [1..ny] $ \s -> do
+            let start = [iz-1, 0, ix-1]
+                count = [1, (ny + s - 1) `div` s, 1]
+                stride = [1, s, 1]
+            eval <- getS nc var start count stride :: HMVRet a
+            case eval of
+              Left e -> assertBool ("read error: " ++ show e) False
+              Right vals -> do
+                let truevals = H.HVector $ C.buildVector
+                               (fromIntegral $ (ny + s - 1) `div` s)
+                               (\iy -> fromIntegral $
+                                       ix + 10 * (iy * s + 1) + 100 * iz)
+                assertBool ("2: value error: " ++ show vals ++
+                            " instead of " ++ show truevals) (vals == truevals)
+      forM_ [1..ny] $ \iy -> do
+        forM_ [1..nz] $ \iz -> do
+          forM_ [1..nx] $ \s -> do
+            let start = [iz-1, iy-1, 0]
+                count = [1, 1, (nx + s - 1) `div` s]
+                stride = [1, 1, s]
+            eval <- getS nc var start count stride :: HMVRet a
+            case eval of
+              Left e -> assertBool ("read error: " ++ show e) False
+              Right vals -> do
+                let truevals = H.HVector $ C.buildVector
+                               (fromIntegral $ (nx + s - 1) `div` s)
+                               (\ix -> fromIntegral $
+                                       (ix * s + 1) + 10 * iy + 100 * iz)
+                assertBool ("3: value error: " ++ show vals ++
+                            " instead of " ++ show truevals) (vals == truevals)
diff --git a/test/test-put.hs b/test/test-put.hs
new file mode 100644
--- /dev/null
+++ b/test/test-put.hs
@@ -0,0 +1,33 @@
+{-# LANGUAGE ScopedTypeVariables #-}
+
+import qualified Data.Map as M
+import Data.Maybe
+import qualified Data.Vector.Storable as SV
+import Control.Error
+import Control.Monad
+import Foreign.C
+
+import Data.NetCDF
+import Data.NetCDF.Vector
+
+main :: IO ()
+main = do
+  let xdim = NcDim "x" 10 False
+      ydim = NcDim "y" 10 False
+      nc = emptyNcInfo "tst-put.nc" #
+           addNcDim xdim #
+           addNcDim ydim #
+           addNcVar (NcVar "x" NcDouble [xdim] M.empty) #
+           addNcVar (NcVar "y" NcDouble [ydim] M.empty) #
+           addNcVar (NcVar "z" NcDouble [xdim, ydim] M.empty)
+  let write nc = do
+        let xvar = fromJust $ ncVar nc "x"
+            yvar = fromJust $ ncVar nc "y"
+            zvar = fromJust $ ncVar nc "z"
+        put nc xvar $ SV.fromList [1..10 :: CDouble]
+        put nc yvar $ SV.fromList [1..10 :: CDouble]
+        put nc zvar $ SV.fromList [ 100 * x + y |
+                                    x <- [1..10 :: CDouble],
+                                    y <- [1..10 :: CDouble] ]
+        return ()
+  withCreateFile nc write (putStrLn . ("ERROR: " ++) . show)