There are at least a dozen online tools that predict phosphorylation sites, using a variety of methods. Most of them simply scan a single input sequence for known motifs of specific kinase families, built using some fairly standard method (e.g. HMMs).
Surprisingly, there does not seem to be a tool to look for conserved sites in multiple sequences. However, people have used this approach: see for example Comparative Analysis Reveals Conserved Protein Phosphorylation Networks Implicated in Multiple Diseases.
I guess it would not be too difficult to build such a tool, using the alignment module of your favourite Bio* project to extract and score the appropriate columns from an alignment.
If you are interested in the conservation then take a look at Claudia Chica's server (http://conscore.embl.de/html/index.html). It is designed to give a score for protein motifs.
More information is available in the phospho.ELM server paper which has the reference to the actual implementation and some information on the usage of the conservation score in phospho.ELM.
Webservice here: conscore.embl.de/CS.wsdl