diff --git a/fsl_glm.cc b/fsl_glm.cc
index 9a5b54ecf47d002ec2df341173a7bfa37d11aa67..f6801ccb6b889ffe4265a728e86f21719348af01 100644
--- a/fsl_glm.cc
+++ b/fsl_glm.cc
@@ -155,12 +155,16 @@ int setup(){
 		
 		// create mask
 		if(fnmask.value()>""){
+			if(debug.value())
+				cout << "Reading mask file " << fnmask.value() << endl;
 			read_volume(mask,fnmask.value());
 			if(!samesize(tmpdata[0],mask)){
 				cerr << "ERROR: Mask image does not match input image" << endl;
 				return 1;
 			};
 		}else{
+			if(debug.value())
+				cout << "Creating mask image" << endl;
 			mask=tmpdata[0]*0.0+1.0;
 			data=tmpdata.matrix(mask);
 			Melodic::update_mask(mask,data);
@@ -169,6 +173,8 @@ int setup(){
 		voxels = data.Ncols();
 			
 		if(perfvn.value()){
+			if(debug.value())
+				cout << "Perform MELODIC variance normalisation (and demeaning)" << endl;			
 			data = remmean(data,1);
 			vnscales = Melodic::varnorm(data);		
 		}
@@ -177,28 +183,44 @@ int setup(){
 		data = read_ascii_matrix(fnin.value());	
 
 	if(fsl_imageexists(fndesign.value())){//read design
+		if(debug.value())
+			cout << "Reading design file "<< fndesign.value()<< endl;
 		volume4D<float> tmpdata;
 		read_volume4D(tmpdata,fndesign.value());
 		if(!samesize(tmpdata[0],mask)){
 			cerr << "ERROR: GLM design does not match input image in size" << endl;
 			return 1;
 		}
+		if(debug.value())
+			cout << "Transposing data" << endl;
 		design = tmpdata.matrix(mask).t();
 		data = data.t();
 	}else{
 		design = read_ascii_matrix(fndesign.value());
 	}
 
-	if(perf_demean.value())
+	if(perf_demean.value()){
+		if(debug.value())
+			cout << "De-meaning the data matrix" << endl;
 		data = remmean(data,1);
-	if(normdat.value())
+	}
+	if(normdat.value()){
+		if(debug.value())
+			cout << "Normalising data matrix to unit std-deviation" << endl;
 		data =  SP(data,ones(data.Nrows(),1)*pow(stdev(data,1),-1));
+	}
 
 	meanR=mean(data,1);
-	if(perf_demean.value())	
+	if(perf_demean.value()){
+		if(debug.value())
+			cout << "De-meaning design matrix" << endl;	
 		design = remmean(design,1);
-	if(normdes.value())
+	}
+	if(normdes.value()){
+		if(debug.value())
+			cout << "Normalising design matrix to unit std-deviation" << endl;
 		design =  SP(design,ones(design.Nrows(),1)*pow(stdev(design,1),-1));
+	}
 	if(fncontrasts.value()>""){//read contrast		
 		contrasts = read_ascii_matrix(fncontrasts.value());
 		if(!(contrasts.Ncols()==design.Ncols())){
@@ -268,6 +290,7 @@ int main(int argc,char *argv[]){
 			options.add(perfvn);
 			options.add(perf_demean);
 			options.add(help);
+			options.add(debug);
 			options.add(outcope);
 			options.add(outz);
 			options.add(outt);