I have made a mistake. It turns out that inmate.f multiplies the coefficients of thermal expansion by the moduli. Thus alp( ) holds the stress coefficients. What you need only do is have
! Membrane forces
norm(1) = (dd(1,1)*eps(1) + dd(1,2)*eps(2) + dd(1,4)*eps(3))*thk
norm(2) = (dd(2,1)*eps(1) + dd(2,2)*eps(2) + dd(2,4)*eps(3))*thk
norm(3) = (dd(4,1)*eps(1) + dd(4,2)*eps(2) + dd(4,4)*eps(3))*thk
! Modify in-plane forces by constant temperature change
norm(1:2) = norm(1:2) - alp(1:2)*d(9)*thk
The rest of the original code is correct.