As g is the same on all boundaries, you do not need subdomains and could simply use the variational formulation:
F = (u*v*dx - T_old*v*dx + dt*inner(a1*grad(u), grad(v))*dx - dt*f*v*dx - dt*g*v*ds)
Otherwise I have no specific comments as what you are doing is quite straight forward.
On the visualization side, you could use Paraview to visualize the data (with eiter XDMF or VTK/PVD files).