A disorder-dependent Gaussian variational approach is applied to the d-dimensional ferromagnetic XY model in a random field. The randomness yields a non extensive contribution to the variational free energy, implying a random mass term in correlation functions. The Imry-Ma low temperature result, concerning the existence (d > 4) or absence (d < 4) of long-range order is obtained in a transparent way. The physical picture which emerges below d = 4 is that of a marginally stable mixture of domains. We also calculate within this variational scheme, disorder dependent correlation functions, as well as the probability distribution of the Imry-Ma domain size.