We measure correlation functions of the nonperturbatively renormalized energy-momentum tensor in $N_f=2+1$ full QCD at finite temperature by applying the gradient flow method both to the gauge and quark fields. Our main interest is to study the conservation law of the energy-momentum tensor and to test whether the linear response relation is properly realized for the entropy density. By using the linear response relation we calculate the specific heat from the correlation function. We adopt the nonperturbatively improved Wilson fermion and Iwasaki gauge action at a fine lattice spacing $=0.07$ fm. In this paper the temperature is limited to a single value $T=232$ MeV. The $u$, $d$ quark mass is rather heavy with $m_pi/m_rho=0.63$ while the $s$ quark mass is set to approximately its physical value.