The programmer manual describes how displacements are transmitted to the element program. Basically, the array ul(ndf,nen,*) gives all the information you need.
ul(ndf,nen,1) are values for the current time, which we usually denote as t_n+1 -- for ndf dof's, nen nodes.
ul(ndf,nen,2) are values of the difference between solutions at t_n+1 and t_n, the previous step.
thus the solution at t_n is obtained as ul(ndf,nen,1) - ul(ndf,nen,2)
For flux boundary conditions on your element you will need to perform the integrations and compute the loading. For thermal problems there are some convec1d.f convec2d.f convec3d.f routines in the thermal subdirectory you can consult for how these are implemented in feap.