  We consider the sample to sample fluctuations that occur in the value of a
thermodynamic quantity $P$ in an ensemble of finite systems with quenched
disorder, at equilibrium. The variance of $P$, $V_{P}$, which characterizes
these fluctuations is calculated as a function of the systems' linear size $l$,
focusing on the behavior at the critical point. The specific model considered
is the bond-disordered Ashkin-Teller model on a square lattice. Using Monte
Carlo simulations, several bond-disordered Ashkin-Teller
  models were examined, including the bond-disordered Ising model and the
bond-disordered four-state Potts model. It was found that far from criticality
  the energy, magnetization, specific heat and susceptibility are strongly self
averaging, that is $V_{P}\sim l^{-d}$ (where $d=2$ is the dimension). At
criticality though, the results indicate that the magnetization $M$ and the
susceptibility $\chi$ are non self averaging, i.e. $\frac{V_{\chi}}{\chi^{2}},
  \frac{V_{M}}{M^{2}}\not \rightarrow 0$. The energy $E$ at criticality is
weakly self averaging, that is $V_{E}\sim l^{-y_{v}}$ with $0<y_{v}<d$. Less
conclusively, and possibly only as a transient behavior, the specific heat too
is found to be weakly self averaging. A phenomenological theory of finite size
scaling for disordered systems is developed. Its main prediction is that when
the specific heat exponent $\alpha<0$ ($\alpha$ of the disordered model) then,
for a quantity $P$ which scales as $l^{\rho}$ at criticality, its variance
$V_{P}$ will scale asymptotically as $l^{2\rho+\frac{\alpha}{\nu}}$. we found
very good agreement between the theory and the data for $V_{\chi}$ and $V_{E}$.
