You will have to program it yourself. The easiest way to do this would be to write a separate element that only implements your boundary conditions:
q . n = h (T^4 - T_o^4)
So if you were doing a 2D problem with 4-node quads you could would be writing a 2-node element that you
would overlay on the edges of your mesh.
See
http://www.ce.berkeley.edu/~sanjay/FEAP/lysmer.html for an example implementation of a boundary condition element, albeit for radiation of elastic waves (and linear for that matter).
