Maybe so:
j1 = 9/10;
j3 = 9/10;
j4 = 1/10;
jz = 1/2;
d1 = 1/10;
d2 = 1/10;
z = 4;
r1 = 1/2*(Cos[x] + Cos[y]);
r2 = Cos[x]*Cos[y];
a = 2*(j1 - j2 + j2*r2 - j1*r1)*z*sa + d1*(2*sa - 1) + jz*sb;
b = 2*(j3 - j4 + j4*r2 - j3*r1)*z*sb + d2*(2*sb - 1) + jz*sa;
c = jz*Sqrt[sa*sb];
ecl = (j2 - j1)*n*z*sa^2 - d1*n*sa^2 + (j4 - j3)*n*z*sb^2 -
d2*n*sb^2 - jz*n*sa*sb;
w1 = ((a - b)*((a + b)^2 - 4*c^2) + ((a + b)^2 + 4*c^2)*
Sqrt[(a + b)^2 - 4*c^2])/(2*((a + b)^2 - 4*c^2));
w2 = ((-a + b)*((a + b)^2 - 4*c^2) + ((a + b)^2 + 4*c^2)*
Sqrt[(a + b)^2 - 4*c^2])/(2*((a + b)^2 - 4*c^2));
q1 = w1/(Exp[w1/t] - 1);
q2 = w2/(Exp[w2/t] - 1);
sa = 1/2;
sb = 1/2;
j2 = 3/10;
n = 1;(* I assume *)
FUNC[T_?NumericQ] :=(1/(4 \[Pi]^2))*NIntegrate[Evaluate[D[q1 + q2, t] /. t -> T], {x, -Pi, Pi}, {y, -Pi, Pi}, Method -> "LocalAdaptive"]
FUNC[2] (* for t = 2 *)
(* 1.59329 *)
Plot[FUNC[T], {T, 1/100, 1}]
