Mathematica 9.0.1 now provides extended FITS file import and export functionality, including arbitrary header fields. However, there is a bug and a hiccup. I've submitted the bug report to Wolfram support. Meanwhile, here are the issues and workarounds:
To demonstrate the bug, we'll generate a FITS file as follows:
Here's metadata that will become additional Header Unit content. We will generate n additional keys, where depending on the value of n, the Data Unit will be stepped upon by a faultily written Header. Mathematica automatically generates 13 header keys upon Export. Each key/value/comment triple occupies an 80-byte "card". The header must be in multiples of 2880 byte "blocks" (equivalent to 36 cards), so padding (ascii char 20) must conditionally be added after the card entries to make this so.
Evidently Mathematica's padding logic is broken for the case where exactly 36 cards (n=23 user + 13 automatic) go into the header. This can be verified by executing the code below, and examining the resulting 3 files (say with DS9 FITS viewer).
One can vary the value of n below and see what happens for n<23, n=23 (busted) and n>23.