How to add a common color-bar to a set of 2d maps produced from evaluating a function along three projections in Mathematica 8.0?

Question

In our simulation, we are assigning neutral hydrogen gas to dark halos as a result of which a distribution of neutral hydrogen is established across each dark halo represented in the form of a distribution of neutral hydrogen column density ($\frac{particles}{cm^{2}}$) across the projected area of the halo. We know that gas absorbers would have different column densities spanning many orders of magnitudes. The final goal of this project is to measure the contribution of each interval of gas system to the total projected area of the dark halo. In other words, we would like to measure the fraction of the total area covered by any interval of contour spanning a range of z values, namely (i) larger than $10^{20.3}$, (ii) between $10^{20.3}$ and $10^{17}$, (iii) between $10^{17}$ and $10^{14}$ and (iv) smaller than $10^{14}$.

I have a text file that contains five columns of data in which columns 1 and 2 correspond to projected distance along dimension 1 and 2 from the center of the dark halo and then columns 3, 4 and 5 are values of the neutral hydrogen column density in projections 0 "z", 1 "x" and 2 "y" respectively. I have transformed this .txt table into actual 2d maps in which the color strength of each pixel would correspond to the magnitude of the function in that pixel (z component in three planes xy, yz, and xz).

My values for columns 3,4 and 5 are mostly zero with few of them on the order of $10^{14}$ to $10^{23}$. I am interested in having three contour levels (less than $10^{14}$, between $10^{14}$ and $10^{17}$, between $10^{17}$ and $10^{20.3}$ and larger than $10^{20.3}). In other words, my 2d map is in fact a set of 3 contour plots that are produced out of .txt table. (Here is the link of a portion of my data: My Data)

Here is my initial code which only yields the 2d maps without any color-bar, axes labels, and fractional coverages:

datam10q = 
  Import["/Users/U of A/Desktop/covering_fractions/m10q.txt", 
   "Table"];
xy10q = datam10q[[All, {1, 2, 3}]];
yz10q = datam10q[[All, {1, 2, 4}]];
xz10q = datam10q[[All, {1, 2, 5}]];
{ListContourPlot[xy10q, 
  Contours -> {{10^14, Yellow}, {10^15.5, Orange}, {10^17, 
     Red}, {10^20.3, Brown}}, 
  ContourShading -> {Yellow, Orange, Red, Brown, Black}, 
  ContourStyle -> None, InterpolationOrder -> 5, 
  BoundaryStyle -> Black, PerformanceGoal -> "Quality"], 
 ListContourPlot[yz10q, 
  Contours -> {{10^14, Yellow}, {10^15.5, Orange}, {10^17, 
     Red}, {10^20.3, Brown}}, 
  ContourShading -> {Yellow, Orange, Red, Brown, Black}, 
  ContourStyle -> None, InterpolationOrder -> 5, 
  BoundaryStyle -> Black, PerformanceGoal -> "Quality"], 
 ListContourPlot[xz10q, 
  Contours -> {{10^14, Yellow}, {10^15.5, Orange}, {10^17, 
     Red}, {10^20.3, Brown}}, 
  ContourShading -> {Yellow, Orange, Red, Brown, Black}, 
  ContourStyle -> None, InterpolationOrder -> 5, 
  BoundaryStyle -> Black, PerformanceGoal -> "Quality"]}

However, the important information that I would like to add to the code and I cannot do it are:

(1) PlotLegend or BarLegend (which seems my Mathematica 8 cannot produce it.)

(2) Axes Label for each one of three maps (namely x-y, y-z, and x-z)

(3) The covering fractions (What fraction of the total area of the map is covered by pixels that have a z value in the range of interval levels between adjacent contours?) I need to print these out somewhere inside the maps and I know that for the two highest density regions, these fractions are on the order of 0.001 and for the third highest region on the order of 0.01.