diff --git a/pals/R/1DAnnualCycle.R b/pals/R/1DAnnualCycle.R
index 5363454..c311a3d 100644
--- a/pals/R/1DAnnualCycle.R
+++ b/pals/R/1DAnnualCycle.R
@@ -8,7 +8,7 @@
 # Gab Abramowitz CCRC, UNSW 2014 (palshelp at gmail dot com)
 #
 AnnualCycle = function(obslabel,acdata,varname,ytext,legendtext,
-	timestepsize,whole,modlabel='no'){
+	timestepsize,whole,plotcolours,modlabel='no'){
 	######
 	errtext = 'ok'
 	metrics = list()
@@ -48,7 +48,6 @@ AnnualCycle = function(obslabel,acdata,varname,ytext,legendtext,
 		}
 	}
 	xloc=c(1:12) # set location of x-coords
-	plotcolours=LineColours()
 	# Plot model output result:
 	yaxmin=min(data_monthly) # y axis minimum in plot
 	yaxmax=max(data_monthly)+0.18*(max(data_monthly)-yaxmin) # y axis maximum in plot
diff --git a/pals/R/1DDiurnalCycle.R b/pals/R/1DDiurnalCycle.R
index 6f75c3e..be2f87d 100644
--- a/pals/R/1DDiurnalCycle.R
+++ b/pals/R/1DDiurnalCycle.R
@@ -7,7 +7,7 @@
 # Gab Abramowitz CCRC, UNSW 2014 (palshelp at gmail dot com)
 #
 DiurnalCycle = function(obslabel,dcdata,varname,ytext,legendtext,
-	timestepsize,whole,modlabel='no',vqcdata=matrix(-1,nrow=1,ncol=1)){
+	timestepsize,whole,plotcolours,modlabel='no',vqcdata=matrix(-1,nrow=1,ncol=1)){
 	errtext = 'ok'
 	metrics = list()
 	if(!whole){ # we need a whole number of years for this to run
@@ -144,7 +144,6 @@ DiurnalCycle = function(obslabel,dcdata,varname,ytext,legendtext,
 	xloc=c(0:(tstepinday-1)) # set location of x-coords in plot
 	yaxmin=min(avday) # y axis minimum in plot
 	yaxmax=max(avday)+(max(avday)-yaxmin)*0.15 # y axis maximum in plot
-	plotcolours=LineColours() 
 	# Now plot each panel:
 	for(k in 1:4){# for each season (DJF, MAM etc)
 		# Plot obs data result:
diff --git a/pals/R/1DPdf.R b/pals/R/1DPdf.R
index b5d1fd5..badfcbd 100644
--- a/pals/R/1DPdf.R
+++ b/pals/R/1DPdf.R
@@ -5,7 +5,7 @@
 # Gab Abramowitz UNSW 2014 (palshelp at gmail dot com)
 
 PALSPdf = function(obslabel,pdfdata,varname,xtext,legendtext,timing,
-	nbins=500,modlabel='no',vqcdata=matrix(-1,nrow=1,ncol=1)){
+	nbins=500,plotcolours,modlabel='no',vqcdata=matrix(-1,nrow=1,ncol=1)){
 	errtext = 'ok'
 	metrics = list()
 	xcut = 1/100
@@ -13,7 +13,6 @@ PALSPdf = function(obslabel,pdfdata,varname,xtext,legendtext,timing,
 	ntsteps = length(pdfdata[,1]) # Number of timesteps in data:
 	tstepinday=86400/timing$tstepsize # number of time steps in a day
 	ndays = ntsteps/tstepinday # number of days in data set
-	plotcolours=LineColours() 
 	xmin = min(pdfdata)
 	xmax = max(pdfdata)
 	allden = list()
@@ -72,14 +71,14 @@ PALSPdf = function(obslabel,pdfdata,varname,xtext,legendtext,timing,
 		scoretext = paste(overl,collapse=', ')
 		text(x=(xlow+(xhigh-xlow)*0.75),y=ymax*0.6,labels=paste('Overlap: ',scoretext,'%',sep=''),pos=4)
 		if(ncurves==2){ # model only
-			metrics[[1]] = list(name='PDFoverlap%',model_value=overl[1])	
+			metrics[[1]] = list(name='%Overlap',model_value=overl[1])	
 		}else if(ncurves==3){
-			metrics[[1]] = list(name='PDFoverlap%',model_value=overl[1],bench_value=list(bench1=overl[2]))	
+			metrics[[1]] = list(name='%Overlap',model_value=overl[1],bench_value=list(bench1=overl[2]))	
 		}else if(ncurves==4){
-			metrics[[1]] = list(name='PDFoverlap%',model_value=overl[1],
+			metrics[[1]] = list(name='%Overlap',model_value=overl[1],
 				bench_value=list(bench1=overl[2],bench2=overl[3]))
 		}else if(ncurves==5){
-			metrics[[1]] = list(name='PDFoverlap%',model_value=overl[1],
+			metrics[[1]] = list(name='%Overlap',model_value=overl[1],
 				bench_value=list(bench1=overl[2],bench2=overl[3],bench3=overl[4]))
 		}
 	}
diff --git a/pals/R/1DScatter.R b/pals/R/1DScatter.R
index c297834..d4cf034 100644
--- a/pals/R/1DScatter.R
+++ b/pals/R/1DScatter.R
@@ -4,65 +4,152 @@
 #
 # Gab Abramowitz UNSW 2014 (palshelp at gmail dot com)
 #
-PALSScatter = function(obslabel,y_data,x_data,varname,vtext,
-	xytext,timestepsize,whole,ebal=FALSE,modlabel='no',
-	vqcdata=matrix(-1,nrow=1,ncol=1)){
+PALSScatter = function(data,varinfo,ebal=FALSE){
 	#
 	errtext = 'ok'
 	metrics = list()
-	ntsteps = length(x_data) # Number of timesteps in data
-	tstepinday=86400/timestepsize # number of time steps in a day
+	ntsteps = length(data$obs$data) # Number of timesteps in data
+	tstepinday=86400/data$obs$timing$tstepsize # number of time steps in a day
 	ndays = ntsteps/tstepinday # number of days in data set
 	# Plot layout:
 	par(mfcol=c(1,2),mar=c(4,4,3,0.5),oma=c(0,0,0,1),
 		mgp=c(2.5,0.7,0),ps=14,tcl=-0.4)
-	if(modlabel=='no'){
-		xtext = paste(xytext[1],':',obslabel)
-		ytext = paste(xytext[2],':',obslabel)
-	}else if(ebal & (modlabel != 'no')){
-		xtext = paste(xytext[1],':',modlabel)
-		ytext = paste(xytext[2],':',modlabel)
+	# Get plot colours, adjuted for any removed benchmarks:
+	plotcols = BenchmarkColours(data$bench,plotobs=FALSE)
+	# Get legend text:
+	legendtext = LegendText(data,plotobs=FALSE)
+	# Define legend line types:
+	linetypes = c(1,3,3,3)
+	# Define x and y-axis text
+	if(ebal){
+		xtext = paste(xytext[1],':',data$model$name)
+		ytext = paste(xytext[2],':',data$model$name)
 	}else{
-		xtext = paste(xytext[1],':',obslabel)
-		ytext = paste(xytext[2],':',modlabel)
+		vtext = bquote(.(tolower(varinfo$PlotName)) ~ ' (' ~.(varinfo$UnitsText) ~ ')')
+		xtext = paste('Observed :',data$obs$name)
+		ytext = paste('Modelled :',data$model$name)
 	}
-	if((vqcdata[1,1] != -1) & (!ebal)){
-		x_mod = x_data[as.logical(vqcdata[,1])]
-		y_mod = y_data[as.logical(vqcdata[,1])]
+	
+	# Prescribe data in scatterplots / regressions, removing qc flagged data if available:
+	if((data$obs$qcexists) & (!ebal)){
+		x_mod = data$obs$data[as.logical(data$obs$qc)]
+		y_mod = data$model$data[as.logical(data$obs$qc)]
+		if(data$bench$exist){
+			b1_mod = data$bench[[ data$bench$index[1] ]]$data[as.logical(data$obs$qc)]
+			if(data$bench$howmany == 2){
+				b2_mod = data$bench[[ data$bench$index[2] ]]$data[as.logical(data$obs$qc)]
+			}else if(data$bench$howmany == 3){
+				b2_mod = data$bench[[ data$bench$index[2] ]]$data[as.logical(data$obs$qc)]
+				b3_mod = data$bench[[ data$bench$index[3] ]]$data[as.logical(data$obs$qc)]
+			}
+		}
 	}else{
-		x_mod = x_data
-		y_mod = y_data
+		x_mod = data$obs$data
+		y_mod = data$model$data
+		if(data$bench$exist){
+			b1_mod = data$bench[[ data$bench$index[1] ]]$data
+			if(data$bench$howmany == 2){
+				b2_mod = data$bench[[ data$bench$index[2] ]]$data
+			}else if(data$bench$howmany == 3){
+				b2_mod = data$bench[[ data$bench$index[2] ]]$data
+				b3_mod = data$bench[[ data$bench$index[3] ]]$data
+			}
+		}
 	}
+	# Define plot boundaries:
 	ymax=max(y_mod)
 	ymin=min(y_mod)
 	xmax=max(x_mod)
 	xmin=min(x_mod)
 	# First plot scatter of per timestep values:
 	plot(x=x_mod,y=y_mod,main=bquote('Per time step' ~ .(vtext)),
-		col='darkblue',xlab=xtext,ylab=ytext,
+		col=plotcols[1],xlab=xtext,ylab=ytext,
 		type='p',pch='.',cex=3,ylim=c(min(ymin,xmin),max(ymax,xmax)),
 		xlim=c(min(ymin,xmin),max(ymax,xmax)))
-	sline = lm(y_mod~x_mod)
-	# Add 1:1 line:
+	# Overplot with better looking points:
+	points(x=x_mod,y=y_mod,pch=20,col=plotcols[1],cex=0.35)
+	# Define per time step regression coefficients:
+	sline = lm(y_mod~x_mod,na.action=na.omit)
+	if(data$bench$exist){
+		b1line = lm(b1_mod~x_mod,na.action=na.omit)
+		if(data$bench$howmany == 2){
+			b2line = lm(b2_mod~x_mod,na.action=na.omit)
+		}else if(data$bench$howmany == 3){
+			b2line = lm(b2_mod~x_mod,na.action=na.omit)
+			b3line = lm(b3_mod~x_mod,na.action=na.omit)
+		}
+	}
+	# Add 1:1 line to plot:
 	abline(a=0,b=1,col='black',lwd=1)
+	# Add benchmark regresion lines to plot:
+	if(data$bench$exist){
+		abline(a=b1line$coefficients[1],b=b1line$coefficients[2],col=plotcols[2],lwd=4,lty=3)
+		if(data$bench$howmany == 2){
+			abline(a=b2line$coefficients[1],b=b2line$coefficients[2],col=plotcols[3],lwd=4,lty=3)
+		}else if(data$bench$howmany == 3){
+			abline(a=b2line$coefficients[1],b=b2line$coefficients[2],col=plotcols[3],lwd=4,lty=3)
+			abline(a=b3line$coefficients[1],b=b3line$coefficients[2],col=plotcols[4],lwd=4,lty=3)
+		}
+	}
 	# Add regresion line to plot
-	abline(a=sline$coefficients[1],b=sline$coefficients[2],col='darkblue',lwd=4)
+	abline(a=sline$coefficients[1],b=sline$coefficients[2],col=plotcols[1],lwd=5)
 	# Add regression parameter text to plot:
-	intdetail = paste('Intercept:',signif(sline$coefficients[1],3))
-	graddetail = paste('Gradient:',signif(sline$coefficients[2],3))
+	if(data$bench$exist){
+		if(data$bench$howmany == 1){
+			intdetail = paste('Intercept:',signif(sline$coefficients[1],2),',',signif(b1line$coefficients[1],2))
+			graddetail = paste('Gradient:',signif(sline$coefficients[2],2),',',signif(b1line$coefficients[2],2))
+		}else if(data$bench$howmany == 2){
+			intdetail = paste('Intercept:',signif(sline$coefficients[1],2),',',
+				signif(b1line$coefficients[1],2),',',signif(b2line$coefficients[1],2))
+			graddetail = paste('Gradient:',signif(sline$coefficients[2],2),',',
+				signif(b1line$coefficients[2],2),',',signif(b2line$coefficients[2],2))
+		}else if(data$bench$howmany == 3){	
+			intdetail = paste('Intercept:',signif(sline$coefficients[1],2),',',
+				signif(b1line$coefficients[1],2),',',signif(b2line$coefficients[1],2),',',signif(b3line$coefficients[1],2))
+			graddetail = paste('Gradient:',signif(sline$coefficients[2],2),',',
+				signif(b1line$coefficients[2],2),',',signif(b2line$coefficients[2],2),',',signif(b3line$coefficients[2],2))
+		}
+	}else{
+		intdetail = paste('Intercept:',signif(sline$coefficients[1],2))
+		graddetail = paste('Gradient:',signif(sline$coefficients[2],2))
+	}
 	yrange = max(ymax,xmax) - min(ymin,xmin)
 	text(x=min(ymin,xmin),y=max(ymax,xmax),labels=intdetail,pos=4)
 	text(x=min(ymin,xmin),y=(min(ymin,xmin)+0.955*yrange),labels=graddetail,pos=4)
-	if((vqcdata[1,1] != -1) & (!ebal)){
-		text(x=(max(xmax,ymax)-yrange*0.5),y=max(ymax,xmax),labels='Gap-filled data removed',pos=4)
+	if((data$obs$qcexists) & (!ebal)){
+		text(x=(max(xmax,ymax)-yrange*0.45),y=max(ymax,xmax),labels='Gap-filled data removed',pos=4)
+	}
+	# Add legend to plot:
+	legend(x=(max(xmax,ymax)-yrange*0.5),y=(max(ymax,xmax)-yrange*0.85),legendtext,
+		lty=linetypes,col=plotcols,lwd=3,bty="n",yjust=0.8)
+	# Define per time step metrics:
+	if(data$bench$exist){
+		if(data$bench$howmany == 1){
+			metrics[[1]] = list(name='Grad',model_value=sline$coefficients[2],
+				bench_value=list(bench1=b1line$coefficients[2]))	
+			metrics[[2]] = list(name='Int',model_value=sline$coefficients[1],
+				bench_value=list(bench1=b1line$coefficients[1]))
+		}else if(data$bench$howmany == 2){
+			metrics[[1]] = list(name='Grad',model_value=sline$coefficients[2],
+				bench_value=list(bench1=b1line$coefficients[2],bench2=b2line$coefficients[2]))	
+			metrics[[2]] = list(name='Int',model_value=sline$coefficients[1],
+				bench_value=list(bench1=b1line$coefficients[1],bench2=b2line$coefficients[1]))
+		}else if(data$bench$howmany == 3){
+			metrics[[1]] = list(name='Grad',model_value=sline$coefficients[2],
+				bench_value=list(bench1=b1line$coefficients[2],bench2=b2line$coefficients[2],
+				bench3=b3line$coefficients[2]))	
+			metrics[[2]] = list(name='Int',model_value=sline$coefficients[1],
+				bench_value=list(bench1=b1line$coefficients[1],bench2=b2line$coefficients[1],
+				bench3=b3line$coefficients[1]))
+		}
+	}else{
+		metrics[[1]] = list(name='Grad',model_value=sline$coefficients[2])	
+		metrics[[2]] = list(name='Int',model_value=sline$coefficients[1])
 	}
-	metrics[[1]] = list(name='Grad',model_value=sline$coefficients[2])	
-	metrics[[2]] = list(name='Int',model_value=sline$coefficients[1])
-	
 	# If this an energy balance plot, add cumulative total:
 	if(ebal){
-		allebal = mean(x_data-y_data)
-		avebal = signif(mean(abs(x_data-y_data)) , 3)
+		allebal = mean(data$obs$data-data$model$data)
+		avebal = signif(mean(abs(data$obs$data-data$model$data)) , 3)
 		alltext = paste('Total imbalance: ',signif(allebal,3),
 			' W/m2 per timestep',sep='')
 		avtext = paste('Mean deviation: ',avebal,' W/m2 per timestep',sep='')
@@ -71,50 +158,142 @@ PALSScatter = function(obslabel,y_data,x_data,varname,vtext,
 		text(x=min(ymin,xmin),y=(min(ymin,xmin)+0.85*yrange),
 			labels=avtext,pos=4)
 	}
-	# Then plot scatter of daily averages:
+	######### Then plot scatter of daily averages #####################
 	# Reshape data into clolumn time-steps-in-day, day rows:
-	x_days=matrix(x_data,ncol=tstepinday,byrow=TRUE)
-	y_days=matrix(y_data,ncol=tstepinday,byrow=TRUE)
+	x_days=matrix(data$obs$data,ncol=tstepinday,byrow=TRUE)
+	y_days=matrix(data$model$data,ncol=tstepinday,byrow=TRUE)
 	xday = c()
 	yday = c()
-	if((vqcdata[1,1] != -1) & (!ebal)){
-		qc_days=matrix(as.logical(vqcdata[,1]),ncol=tstepinday,byrow=TRUE)
+	if(data$bench$exist){
+		b1_days=matrix(data$bench[[ data$bench$index[1] ]]$data,ncol=tstepinday,byrow=TRUE)
+		b1day = c()
+		if(data$bench$howmany == 2){
+			b2_days=matrix(data$bench[[ data$bench$index[2] ]]$data,ncol=tstepinday,byrow=TRUE)
+			b2day = c()
+		}else if(data$bench$howmany == 3){
+			b2_days=matrix(data$bench[[ data$bench$index[2] ]]$data,ncol=tstepinday,byrow=TRUE)
+			b3_days=matrix(data$bench[[ data$bench$index[3] ]]$data,ncol=tstepinday,byrow=TRUE)
+			b2day = c()
+			b3day = c()
+		}	
+	}
+	# Define daily mean time series, potentially removing qc time steps
+	if((data$obs$qcexists) & (!ebal)){
+		qc_days=matrix(as.logical(data$obs$qc),ncol=tstepinday,byrow=TRUE)
 		for(i in 1:ndays){
 			xday[i] = mean(x_days[i,qc_days[i,]])
 			yday[i] = mean(y_days[i,qc_days[i,]])
 		}
+		if(data$bench$exist){
+			if(data$bench$howmany == 1){
+				for(i in 1:ndays){
+					b1day[i] = mean(b1_days[i,qc_days[i,]])
+				}
+			}else if(data$bench$howmany == 2){
+				for(i in 1:ndays){
+					b1day[i] = mean(b1_days[i,qc_days[i,]])
+					b2day[i] = mean(b2_days[i,qc_days[i,]])
+				}
+			}else if(data$bench$howmany == 3){
+				for(i in 1:ndays){
+					b1day[i] = mean(b1_days[i,qc_days[i,]])
+					b2day[i] = mean(b2_days[i,qc_days[i,]])
+					b3day[i] = mean(b3_days[i,qc_days[i,]])
+				}
+			}
+		}	
 	}else if(ebal){
-          ebalday = c()
-          for(i in 1:ndays){
-            xday[i] = mean(x_days[i,])
-            yday[i] = mean(y_days[i,])
+    	ebalday = c()
+    	for(i in 1:ndays){
+        	xday[i] = mean(x_days[i,])
+           	yday[i] = mean(y_days[i,])
             ebalday[i] = sum(x_days[i,]-y_days[i,])
-          }
-        }else{ # not not an energy balance plot and QC data are missing
+        }
+   	}else{ # not not an energy balance plot and QC data are missing
 		for(i in 1:ndays){
 			xday[i] = mean(x_days[i,])
 			yday[i] = mean(y_days[i,])
 		}
+		if(data$bench$exist){
+			if(data$bench$howmany == 1){
+				for(i in 1:ndays){
+					b1day[i] = mean(b1_days[i,])
+				}
+			}else if(data$bench$howmany == 2){
+				for(i in 1:ndays){
+					b1day[i] = mean(b1_days[i,])
+					b2day[i] = mean(b2_days[i,])
+				}
+			}else if(data$bench$howmany == 3){
+				for(i in 1:ndays){
+					b1day[i] = mean(b1_days[i,])
+					b2day[i] = mean(b2_days[i,])
+					b3day[i] = mean(b3_days[i,])
+				}
+			}
+		}
 	}
 	ymax=max(yday,na.rm=TRUE)
 	ymin=min(yday,na.rm=TRUE)
 	xmax=max(xday,na.rm=TRUE)
 	xmin=min(xday,na.rm=TRUE)
-	plot(x=xday,y=yday,main=bquote('Daily average' ~ .(vtext)),col='darkblue',
+	plot(x=xday,y=yday,main=bquote('Daily average' ~ .(vtext)),col=plotcols[1],
 		xlab=xtext,ylab=ytext,type='p',pch='.',cex=3,
 		ylim=c(min(ymin,xmin),max(ymax,xmax)),
 		xlim=c(min(ymin,xmin),max(ymax,xmax)))
+	# Overplot with better looking points:
+	points(x=xday,y=yday,pch=20,col=plotcols[1],cex=0.35)
+	# Define daily regression coefficients:
 	sline = lm(yday~xday,na.action=na.omit)
+	if(data$bench$exist){
+		b1line = lm(b1day~xday,na.action=na.omit)
+		if(data$bench$howmany == 2){
+			b2line = lm(b2day~xday,na.action=na.omit)
+		}else if(data$bench$howmany == 3){
+			b2line = lm(b2day~xday,na.action=na.omit)
+			b3line = lm(b3day~xday,na.action=na.omit)
+		}
+	}
 	# Add 1:1 line:
 	abline(a=0,b=1,col='black',lwd=1)
+	# Add benchmark regresion lines to plot:
+	if(data$bench$exist){
+		abline(a=b1line$coefficients[1],b=b1line$coefficients[2],col=plotcols[2],lwd=4,lty=3)
+		if(data$bench$howmany == 2){
+			abline(a=b2line$coefficients[1],b=b2line$coefficients[2],col=plotcols[3],lwd=4,lty=3)
+		}else if(data$bench$howmany == 3){
+			abline(a=b2line$coefficients[1],b=b2line$coefficients[2],col=plotcols[3],lwd=4,lty=3)
+			abline(a=b3line$coefficients[1],b=b3line$coefficients[2],col=plotcols[4],lwd=4,lty=3)
+		}
+	}
 	# Add regresion line to plot
-	abline(a=sline$coefficients[1],b=sline$coefficients[2],col='darkblue',lwd=4)
+	abline(a=sline$coefficients[1],b=sline$coefficients[2],col=plotcols[1],lwd=5)
 	# Add regression parameter text to plot:
-	intdetail = paste('Intercept:',signif(sline$coefficients[1],3))
-	graddetail = paste('Gradient:',signif(sline$coefficients[2],3))
+	if(data$bench$exist){
+		if(data$bench$howmany == 1){
+			intdetail = paste('Intercept:',signif(sline$coefficients[1],2),',',signif(b1line$coefficients[1],2))
+			graddetail = paste('Gradient:',signif(sline$coefficients[2],2),',',signif(b1line$coefficients[2],2))
+		}else if(data$bench$howmany == 2){
+			intdetail = paste('Intercept:',signif(sline$coefficients[1],2),',',
+				signif(b1line$coefficients[1],2),',',signif(b2line$coefficients[1],2))
+			graddetail = paste('Gradient:',signif(sline$coefficients[2],2),',',
+				signif(b1line$coefficients[2],2),',',signif(b2line$coefficients[2],2))
+		}else if(data$bench$howmany == 3){	
+			intdetail = paste('Intercept:',signif(sline$coefficients[1],2),',',
+				signif(b1line$coefficients[1],2),',',signif(b2line$coefficients[1],2),',',signif(b3line$coefficients[1],2))
+			graddetail = paste('Gradient:',signif(sline$coefficients[2],2),',',
+				signif(b1line$coefficients[2],2),',',signif(b2line$coefficients[2],2),',',signif(b3line$coefficients[2],2))
+		}
+	}else{
+		intdetail = paste('Intercept:',signif(sline$coefficients[1],2))
+		graddetail = paste('Gradient:',signif(sline$coefficients[2],2))
+	}
 	yrange = max(ymax,xmax) - min(ymin,xmin)
 	text(x=min(ymin,xmin),y=max(ymax,xmax),labels=intdetail,pos=4)
 	text(x=min(ymin,xmin),y=(min(ymin,xmin)+0.955*yrange),labels=graddetail,pos=4)
+	# Add legend to plot:
+	legend(x=(max(xmax,ymax)-yrange*0.5),y=(max(ymax,xmax)-yrange*0.85),legendtext,
+		lty=linetypes,col=plotcols,lwd=3,bty="n",yjust=0.8)
 	# If this an energy balance plot, add cumulative total:
 	if(ebal){
 		allebal = sum(ebalday)*timestepsize/3600/ndays # in Watt-hours
@@ -126,12 +305,33 @@ PALSScatter = function(obslabel,y_data,x_data,varname,vtext,
 			labels=alltext,pos=4)
 		text(x=min(ymin,xmin),y=(min(ymin,xmin)+0.85*yrange),
 			labels=avtext,pos=4)
-	}else if(vqcdata[1,1] != -1){
-		text(x=(max(xmax,ymax)-yrange*0.5),y=max(ymax,xmax),labels='Gap-filled data removed',pos=4)
+	}else if(data$obs$qcexists){
+		text(x=(max(xmax,ymax)-yrange*0.45),y=max(ymax,xmax),labels='Gap-filled data removed',pos=4)
 	}
 	# Record metrics:
-	metrics[[3]] = list(name='DailyGrad',model_value=sline$coefficients[2])	
-	metrics[[4]] = list(name='DailyInt',model_value=sline$coefficients[1])	
+	if(data$bench$exist){
+		if(data$bench$howmany == 1){
+			metrics[[3]] = list(name='DailyGrad',model_value=sline$coefficients[2],
+				bench_value=list(bench1=b1line$coefficients[2]))	
+			metrics[[4]] = list(name='DailyInt',model_value=sline$coefficients[1],
+				bench_value=list(bench1=b1line$coefficients[1]))
+		}else if(data$bench$howmany == 2){
+			metrics[[3]] = list(name='DailyGrad',model_value=sline$coefficients[2],
+				bench_value=list(bench1=b1line$coefficients[2],bench2=b2line$coefficients[2]))	
+			metrics[[4]] = list(name='DailyInt',model_value=sline$coefficients[1],
+				bench_value=list(bench1=b1line$coefficients[1],bench2=b2line$coefficients[1]))
+		}else if(data$bench$howmany == 3){
+			metrics[[3]] = list(name='DailyGrad',model_value=sline$coefficients[2],
+				bench_value=list(bench1=b1line$coefficients[2],bench2=b2line$coefficients[2],
+				bench3=b3line$coefficients[2]))	
+			metrics[[4]] = list(name='DailyInt',model_value=sline$coefficients[1],
+				bench_value=list(bench1=b1line$coefficients[1],bench2=b2line$coefficients[1],
+				bench3=b3line$coefficients[1]))
+		}
+	}else{
+		metrics[[3]] = list(name='DailyGrad',model_value=sline$coefficients[2])	
+		metrics[[4]] = list(name='DailyInt',model_value=sline$coefficients[1])	
+	}	
 
 	result=list(err=TRUE,errtext=errtext, metrics=metrics)
 	return(result)
diff --git a/pals/R/1DTaylorDiagram.R b/pals/R/1DTaylorDiagram.R
index e4ce955..d7736d0 100644
--- a/pals/R/1DTaylorDiagram.R
+++ b/pals/R/1DTaylorDiagram.R
@@ -4,24 +4,63 @@
 #
 # Gab Abramowitz UNSW 2014 palshelp@gmail.com
 #
-TaylorDiagram = function(ObsLabel,mod_data,obs_data,varname,xtext,
-	timestepsize,whole,modlabel){
+TaylorDiagram = function(data,varinfo,plotcolours){
 	library(plotrix) # load package with Taylor diagram
 	errtext = 'ok'
+	# Get colours:
+	plotcolours = BenchmarkColours(data$bench,plotobs=TRUE)
+	plot2colours = DesaturateColours(plotcolours,sat=0.3)
+	# Get text for legend:
+	legendtext = LegendText(data,plotobs=TRUE)
+	# First report metrics:
 	metrics = list()
-	metrics[[1]] = list(name='RMSE',model_value=sqrt(mean(((as.vector(mod_data)-as.vector(obs_data))^2),na.rm=TRUE)))
-	taylor.diagram(ref=obs_data,model=mod_data,xlab=xtext,pcex=2,
-		main=paste(varname[1],' Taylor diagram:  Obs - ',ObsLabel,
-		'  Mod - ',modlabel, sep=''),ref.sd=TRUE,show.gamma=TRUE,
-		cex.lab=1.2,cex.axis=1.1,cex.main=1.2,col='blue',
-		sd.arcs=TRUE,mgp = c(2.7,1,0))
-	taylor.diagram(ref=obs_data,model=obs_data,pcex=2,add=TRUE,
-		ref.sd=TRUE,show.gamma=TRUE,col='blue',pch=1)
-	# Calculate daily values:
-	tstepinday=86400/timestepsize # number of time steps in a day
+	if(data$bench$exist){
+		if(data$bench$howmany == 1){
+			metrics[[1]] = list(name='Correlation',
+				model_value=cor(as.vector(data$model$data), as.vector(data$obs$data)),
+				bench_value=list(
+				bench1=cor(as.vector(data$bench[[ data$bench$index[1] ]]$data), as.vector(data$obs$data)) ))	
+		}else if(data$bench$howmany == 2){
+			metrics[[1]] = list(name='Correlation',
+				model_value=cor(as.vector(data$model$data),as.vector(data$obs$data)),
+				bench_value=list(
+				bench1=cor(as.vector(data$bench[[ data$bench$index[1] ]]$data), as.vector(data$obs$data)),
+				bench2=cor(as.vector(data$bench[[ data$bench$index[2] ]]$data), as.vector(data$obs$data)) ))
+		}else if(data$bench$howmany == 3){
+			metrics[[1]] = list(name='Correlation',
+				model_value=cor(as.vector(data$model$data), as.vector(data$obs$data)),
+				bench_value=list(
+				bench1=cor(as.vector(data$bench[[ data$bench$index[1] ]]$data), as.vector(data$obs$data)),
+				bench2=cor(as.vector(data$bench[[ data$bench$index[2] ]]$data), as.vector(data$obs$data)),
+				bench3=cor(as.vector(data$bench[[ data$bench$index[3] ]]$data), as.vector(data$obs$data)) ))	
+		}
+	 }else{
+		metrics[[1]] = list(name='Correlation',
+			model_value=cor(as.vector(data$model$data),as.vector(data$obs$data)))
+	}
+	# Text for x-axis:
+	xtext=bquote(.(varinfo$PlotName) ~ ' (' ~ .(varinfo$UnitsText) ~ ')')
+	# First plot model on Taylor plot:
+	taylor.diagram(ref=as.vector(data$obs$data),model=as.vector(data$model$data),xlab=xtext,pcex=2,
+		main=paste(varinfo$name[1],' Taylor diagram:  Obs - ',data$obs$name,'  Mod - ',data$model$name, sep=''),
+		ref.sd=TRUE,show.gamma=TRUE,cex.lab=1.2,cex.axis=1.1,cex.main=1.2,
+		col=plotcolours[2],sd.arcs=TRUE)	
+	# Then obs reference point on x-axis:
+	taylor.diagram(ref=as.vector(data$obs$data),model=as.vector(data$obs$data),pcex=2,add=TRUE,
+		ref.sd=TRUE,show.gamma=TRUE,col=plotcolours[1])
+	# Then plot benchmarks:
+	if(data$bench$exist){
+		for(b in 1:data$bench$howmany){	
+			taylor.diagram(ref=as.vector(data$obs$data),model=as.vector(data$bench[[ data$bench$index[b] ]]$data),
+				pcex=2,add=TRUE,ref.sd=TRUE,show.gamma=TRUE,col=plotcolours[(b+2)])	
+		}
+	}
+	
+	###### Now calculate daily values ##############
+	tstepinday=86400/data$obs$timing$tstepsize # number of time steps in a day
 	# Reshape data into clolumn hours, day rows:
-	mod_days=matrix(mod_data,ncol=tstepinday,byrow=TRUE) 
-	obs_days=matrix(obs_data,ncol=tstepinday,byrow=TRUE) 
+	mod_days=matrix(data$model$data,ncol=tstepinday,byrow=TRUE) 
+	obs_days=matrix(data$obs$data,ncol=tstepinday,byrow=TRUE) 
 	ndays=length(mod_days[,1]) # find number of days in data set
 	nyears=as.integer(ndays/365) # find # years in data set
 	mavday=c() # initialise
@@ -30,64 +69,68 @@ TaylorDiagram = function(ObsLabel,mod_data,obs_data,varname,xtext,
 		mavday[i] = mean(mod_days[i,])
 		oavday[i] = mean(obs_days[i,])
 	}
-	# Plot daily averages:
+	# Reshape benchmark data:
+	if(data$bench$exist){
+		bavday=matrix(NA,data$bench$howmany,ndays) # initialise
+		for(b in 1:data$bench$howmany){	
+			bench_days=matrix(data$bench[[ data$bench$index[b] ]]$data,ncol=tstepinday,byrow=TRUE)
+			for(i in 1:ndays){
+				bavday[b,i] = mean(bench_days[i,])
+			}
+		}
+	}	
+	# Plot daily averages of model:
 	taylor.diagram(ref=oavday,model=mavday,pcex=2,add=TRUE,
-	ref.sd=TRUE,show.gamma=TRUE,col='red')
+		ref.sd=TRUE,show.gamma=TRUE,col=plot2colours[2],pch=15)
+	# Plot daily averages of obs reference circle:
 	taylor.diagram(ref=oavday,model=oavday,pcex=2,add=TRUE,
-	ref.sd=TRUE,show.gamma=TRUE,col='red',pch=1)
-	# If we have whole years of data, plot monthly values as well
-	if(whole){ 
-		# Calculate monthly values:
-		month_start=c()
-		month_start[1]=1    # Jan
-		month_start[2]=32   # Feb
-		month_start[3]=60   # Mar
-		month_start[4]=91   # Apr
-		month_start[5]=121  # May
-		month_start[6]=152  # Jun
-		month_start[7]=182  # Jul
-		month_start[8]=213  # Aug
-		month_start[9]=244  # Sep
-		month_start[10]=274 # Oct
-		month_start[11]=305 # Nov
-		month_start[12]=335 # Dec
-		month_start[13]=366 # i.e. beginning of next year
-		mod_monthly=c() # initialise monthly averages
-		obs_monthly=c()   # initialise monthly averages
-		# Transform daily means into monthly means:
-		for(l in 1:12){ # for each month
-			month_length=month_start[l+1]-month_start[l]
-			mod_month=0 # initialise
-			obs_month=0   # initialise
-			for(k in 1:nyears){ # for each year of data set
-				# Add all daily averages for a given month
-				# over all data set years:
-				mod_month = mod_month + 
-					sum(mavday[(month_start[l]+(k-1)*365):
-					(month_start[l+1]-1 +(k-1)*365) ] )
-				obs_month = obs_month + 
-					sum(oavday[(month_start[l]+(k-1)*365):
-					(month_start[l+1]-1 +(k-1)*365) ] )
-			}
-			# then divide by the total number of days added above:
-			mod_monthly[l]=mod_month/(month_length*nyears)
-			obs_monthly[l]=obs_month/(month_length*nyears)
+		ref.sd=TRUE,show.gamma=TRUE,col=plot2colours[1],pch=15)
+	# Plot daily averages of benchmarks:
+	if(data$bench$exist){
+		for(b in 1:data$bench$howmany){	
+			taylor.diagram(ref=oavday,model=bavday[b,],pcex=2,add=TRUE,
+				ref.sd=TRUE,show.gamma=TRUE,col=plot2colours[(b+2)],pch=15)	
 		}
-		# Plot monthly averages:
-		taylor.diagram(ref=obs_monthly,model=mod_monthly,pcex=2,add=TRUE,
-		ref.sd=TRUE,show.gamma=TRUE,col='green')
-		taylor.diagram(ref=obs_monthly,model=obs_monthly,pcex=2,add=TRUE,
-		ref.sd=TRUE,show.gamma=TRUE,col='green',pch=1)
-		# Add legend:
-		lpos = 1.2*max(sd(obs_data),sd(mod_data))
-		legend(lpos,1.3*lpos,col=c('blue','red','green'),pch=19,
-			legend=c('per timestep','daily averages','monthly averages'))
-	}else{
-		# Add legend:
-		lpos = 1.2*max(sd(obs_data),sd(mod_data))
-		legend(lpos,1.3*lpos,col=c('blue','red'),pch=19,
-			legend=c('per timestep','daily averages'))
 	}
+	
+	# # If we have whole years of data, plot monthly values as well
+	# if(data$obs$timing$whole){
+		# month = getMonthDays()
+		# mod_monthly=c() # initialise monthly averages
+		# obs_monthly=c()   # initialise monthly averages
+		# # Transform daily means into monthly means:
+		# for(l in 1:12){ # for each month
+			# mod_month=0 # initialise
+			# obs_month=0   # initialise
+			# for(k in 1:nyears){ # for each year of data set
+				# # Add all daily averages for a given month
+				# # over all data set years:
+				# mod_month = mod_month + 
+					# sum(mavday[(month$start[l]+(k-1)*365):
+					# (month$start[l+1]-1 +(k-1)*365) ] )
+				# obs_month = obs_month + 
+					# sum(oavday[(month$start[l]+(k-1)*365):
+					# (month$start[l+1]-1 +(k-1)*365) ] )
+			# }
+			# # then divide by the total number of days added above:
+			# mod_monthly[l]=mod_month/(month$length[l]*nyears)
+			# obs_monthly[l]=obs_month/(month$length[l]*nyears)
+		# }
+		# # Plot monthly averages:
+		# taylor.diagram(ref=obs_monthly,model=mod_monthly,pcex=2,add=TRUE,
+		# ref.sd=TRUE,show.gamma=TRUE,col='green')
+		# taylor.diagram(ref=obs_monthly,model=obs_monthly,pcex=2,add=TRUE,
+		# ref.sd=TRUE,show.gamma=TRUE,col='green',pch=1)
+		# # Add legend:
+		# lpos = 1.2*max(sd(data$obs$data),sd(data$model$data))
+		# legend(lpos,1.3*lpos,col=c('blue','red','green'),pch=19,
+			# legend=c('per timestep','daily averages','monthly averages'))
+	# }else{
+	# Add legend:
+	legend('topright',col=c(plotcolours,plot2colours[1]),
+		pch=c(rep.int(19,(2+data$bench$howmany)),15),
+		inset=c(-0.1,-0.05),legend=c(legendtext,'Daily averages'))		 
+	# }
 	result=list(err=FALSE,errtext=errtext,metrics=metrics)
 	return(result)
 } # End function TaylorDiagram
\ No newline at end of file
diff --git a/pals/R/1DTimeseries.R b/pals/R/1DTimeseries.R
index e4ba47b..7360836 100644
--- a/pals/R/1DTimeseries.R
+++ b/pals/R/1DTimeseries.R
@@ -5,7 +5,7 @@
 # Gab Abramowitz UNSW 2014 (palshelp at gmail dot com)
 
 Timeseries = function(obslabel,tsdata,varname,ytext,legendtext,
-	plotcex,timing,smoothed=FALSE,winsize=1,modlabel='no',
+	plotcex,timing,smoothed=FALSE,winsize=1,plotcolours,modlabel='no',
 	vqcdata=matrix(-1,nrow=1,ncol=1)){
 	#
 	errtext = 'ok'
@@ -15,7 +15,6 @@ Timeseries = function(obslabel,tsdata,varname,ytext,legendtext,
 	tstepinday=86400/timing$tstepsize # number of time steps in a day
 	ndays = ntsteps/tstepinday # number of days in data set
 	nyears=as.integer(ndays/365) # find # years in data set
-		plotcolours=LineColours() 
 	# x-axis labels:
 	xxat=c()
 	xxlab=c()
@@ -66,26 +65,35 @@ Timeseries = function(obslabel,tsdata,varname,ytext,legendtext,
 					sum(abs(mean(tsdata[,1]) - tsdata[,1]))
 			}
 			# Report NME metric:
-			metricname = paste('NME',winsize,'-dayAvs',sep='')
+			metricname = paste('NME',winsize,'day',sep='')
 			if(ncurves==2){ # model only
-				metrics[[1]] = list(name=metricname,model_value=smoothscore[1])	
+				metrics[[1]] = list(name='Bias',model_value=mean(tsdata[,2]-tsdata[,1],na.rm=TRUE))	
 				metrics[[2]] = list(name='NME',model_value=allscore[1])	
+				metrics[[3]] = list(name=metricname,model_value=smoothscore[1])	
 			}else if(ncurves==3){
-				metrics[[1]] = list(name=metricname,model_value=smoothscore[1],
-					bench_value=list(bench1=smoothscore[2]))	
+				metrics[[1]] = list(name='Bias',model_value=mean(tsdata[,2]-tsdata[,1],na.rm=TRUE),
+					bench_value=list(bench1=mean(tsdata[,3]-tsdata[,1],na.rm=TRUE) ))
 				metrics[[2]] = list(name='NME',model_value=allscore[1],bench_value=list(bench1=allscore[2]))	
+				metrics[[3]] = list(name=metricname,model_value=smoothscore[1],
+					bench_value=list(bench1=smoothscore[2]))	
 			}else if(ncurves==4){
-				metrics[[1]] = list(name=metricname,model_value=smoothscore[1],
-					bench_value=list(bench1=smoothscore[2],bench2=smoothscore[3]))	
+				metrics[[1]] = list(name='Bias',model_value=mean(tsdata[,2]-tsdata[,1],na.rm=TRUE),
+					bench_value=list(bench1=mean(tsdata[,3]-tsdata[,1],na.rm=TRUE),
+					bench2=mean(tsdata[,4]-tsdata[,1],na.rm=TRUE) ))
 				metrics[[2]] = list(name='NME',model_value=allscore[1],
 					bench_value=list(bench1=allscore[2],bench2=allscore[3]))
+				metrics[[3]] = list(name=metricname,model_value=smoothscore[1],
+					bench_value=list(bench1=smoothscore[2],bench2=smoothscore[3]))	
 			}else if(ncurves==5){
-				metrics[[1]] = list(name=metricname,model_value=smoothscore[1],
-					bench_value=list(bench1=smoothscore[2],bench2=smoothscore[3],bench3=smoothscore[4]))	
+				metrics[[1]] = list(name='Bias',model_value=mean(tsdata[,2]-tsdata[,1],na.rm=TRUE),
+					bench_value=list(bench1=mean(tsdata[,3]-tsdata[,1],na.rm=TRUE),
+					bench2=mean(tsdata[,4]-tsdata[,1],na.rm=TRUE),
+					bench3=mean(tsdata[,5]-tsdata[,1],na.rm=TRUE) ))
 				metrics[[2]] = list(name='NME',model_value=allscore[1],
 					bench_value=list(bench1=allscore[2],bench2=allscore[3],bench3=allscore[4]))
+				metrics[[3]] = list(name=metricname,model_value=smoothscore[1],
+					bench_value=list(bench1=smoothscore[2],bench2=smoothscore[3],bench3=smoothscore[4]))	
 			}
-			
 		}
 		for(l in 1:nyears){
 			xxat[(2*l-1)] = (l-1)*365 + 1
@@ -139,6 +147,7 @@ Timeseries = function(obslabel,tsdata,varname,ytext,legendtext,
 				labels=paste(qcpc,'% of observed ',varname[1],' is gap-filled:',sep=''))
 		}		
 	}else{
+		# this code not functioning but kept for future modification:
 		yvalmin = signif(min(tsdata),3)
 		yvalmax = signif(max(tsdata),3)
 		datamean = signif(mean(tsdata[,1]),3)
diff --git a/pals/R/2DAus.R b/pals/R/2DAus.R
new file mode 100644
index 0000000..cd6ec44
--- /dev/null
+++ b/pals/R/2DAus.R
@@ -0,0 +1,71 @@
+# 2DAus.R
+#
+# Gab Abramowitz, UNSW, 2015, gabsun at gmail dot com
+#
+PlotAus = function(lon,lat,data,meanval,sdval,varname,unitstxt,longvarname,zrange,zcols,title,textloc,suppressunits=FALSE){
+	# Generates a gridded heatmap style plot for Australia with map, based on input lat/lon.
+	library(maps)
+	library(mapdata)
+	library(fields) # for image.plot
+	errtext = 'ok'
+	# Decide location of plot for text placement:
+	if(textloc=='bottomright'){
+		textloc1 = c((lon[1] + (lon[length(lon)] - lon[1])*0.9), (lat[1] + (lat[length(lat)] - lat[1])*0.19) )
+		textloc2 = c((lon[1] + (lon[length(lon)] - lon[1])*0.9), (lat[1] + (lat[length(lat)] - lat[1])*0.13) )
+		textloc3 = c((lon[1] + (lon[length(lon)] - lon[1])*0.9), (lat[1] + (lat[length(lat)] - lat[1])*0.07) )
+	}else if(textloc=='topleft'){
+		textloc1 = c((lon[1] + (lon[length(lon)] - lon[1])*0.15), (lat[1] + (lat[length(lat)] - lat[1])*0.94) )
+		textloc2 = c((lon[1] + (lon[length(lon)] - lon[1])*0.15), (lat[1] + (lat[length(lat)] - lat[1])*0.88) )
+		textloc3 = c((lon[1] + (lon[length(lon)] - lon[1])*0.15), (lat[1] + (lat[length(lat)] - lat[1])*0.82) )
+	}
+	# Plot:	
+	image.plot(lon,lat,data,xlab='Longitude',ylab='Latitude',col=zcols,zlim=zrange,legend.mar=5.5)
+	map('worldHires',add=TRUE,wrap=TRUE,xlim=c(min(lon),max(lon)),ylim=c(min(lat),max(lat))) # Add map
+	title(title) 
+	if(!suppressunits){
+		text(x=textloc1[1],y=textloc1[2],labels=unitstxt)
+	}
+	text(x=textloc2[1],y=textloc2[2],labels=paste('Mean:',signif(meanval,3)))
+	text(x=textloc3[1],y=textloc3[2],labels=paste('SD:',signif(sdval,2)))
+	return(errtext)
+}
+
+DensityLocationAus = function(npanels){
+	# Simply prescribes locations for density inset on gridded Australia map plots
+	density_location = list()
+	if(npanels == 1){
+		density_location = list()
+		density_location[[1]] = c(0.1,0.5,0.16,0.35)
+	}else if(npanels == 2){
+		density_location = list()
+		density_location[[1]] = c(0.07,0.3,0.17,0.37)
+		density_location[[2]] = c(0.57,0.8,0.17,0.37)
+	}else if(npanels == 4){
+		density_location = list()
+		density_location[[1]] = c(0.06,0.3,0.59,0.68)
+		density_location[[2]] = c(0.56,0.8,0.59,0.68)
+		density_location[[3]] = c(0.06,0.3,0.09,0.18)
+		density_location[[4]] = c(0.56,0.8,0.09,0.18)
+	}else if(npanels==6){
+		density_location = list()
+		density_location[[1]] = c(0.04,0.2,0.58,0.67)
+		density_location[[2]] = c(0.38,0.54,0.58,0.67)
+		density_location[[3]] = c(0.72,0.87,0.58,0.67)
+		density_location[[4]] = c(0.04,0.2,0.08,0.17)
+		density_location[[5]] = c(0.38,0.54,0.08,0.17)
+		density_location[[6]] = c(0.72,0.87,0.08,0.17)
+	}
+	return(density_location)
+}
+
+TextLocationAus = function(npanels){
+	if(npanels<5){
+		text_location = 'bottomright'
+	}else{
+		text_location = 'topleft'
+	}
+}
+
+
+
+
diff --git a/pals/R/2DAusMean.R b/pals/R/2DAusMean.R
deleted file mode 100644
index cbedc80..0000000
--- a/pals/R/2DAusMean.R
+++ /dev/null
@@ -1,445 +0,0 @@
-# 2DAusMean.R
-#
-# Gab Abramowitz, UNSW, 2014, gabsun at gmail dot com
-#
-SpatialAus = function(model,obs,bench,varname,unitstxt,longvarname,metrics,plottype){
-	errtext = 'ok'
-	metrics = list()
-	density_cut = 1/200
-	# Calculate number of map panels:
-	npanels = bench$howmany + 3
-	# Plot layout:
-	if(npanels <= 4){
-		par(mfcol=c(2,2) ,mar=c(3,3,3,0.5),oma=c(0,0,0,1),mgp=c(1.8,0.5,0),ps=15,tcl=-0.4)
-		density_location = DensityLocationAus(4)
-		textloc='bottomright'
-	}else{
-		par(mfcol=c(2,3) ,mar=c(3,3,3,0.5),oma=c(1,0,0.5,1),mgp=c(1.8,0.5,0),ps=18,tcl=-0.2)
-		density_location = list()
-		density_location = DensityLocationAus(6)
-		textloc='topleft'
-	}
-	if(plottype=='TimeMean'){
-		# Calculate time means:
-		modelt = TimeMean(model$data)
-		obst = TimeMean(obs$data) + modelt - modelt # to make sure ocean areas not included
-		modelbias = mean(modelt-obst,na.rm=TRUE)
-	}else if(plottype=='TimeSD'){
-		modelt = TimeSD(model$data)
-		obst = TimeSD(obs$data) + modelt - modelt # to make sure ocean areas not included
-		modelSDbias = mean(modelt-obst,na.rm=TRUE)
-	}else{
-		result = list(errtext = paste('Unknown plot type \'',plottype,'\' requested in function SpatialAus.',sep=''),err=TRUE)
-		return(result)
-	}
-	zrange = c(min(modelt,obst,na.rm=TRUE),max(modelt,obst,na.rm=TRUE))
-	diffrange = c(min(modelt-obst,na.rm=TRUE),max(modelt-obst,na.rm=TRUE))
-	if(plottype=='TimeMean'){
-		if(bench$howmany == 1){
-			# Get benchmark data, noting there may have been other benchmarks that failed (so use bench$index)
-			bench1t = TimeMean(bench[[ bench$index[1] ]]$data)
-			bench1bias = mean(bench1t-obst,na.rm=TRUE)
-			diffrange = c(min(modelt-obst,bench1t-obst,na.rm=TRUE),
-				max(modelt-obst,bench1t-obst,na.rm=TRUE))
-			metrics[[1]] = list(name='TimeSpaceBias',model_value=modelbias,bench_value=list(bench1=bench1bias))	
-		}else if(bench$howmany == 2){
-			# Get benchmark data, noting there may have been other benchmarks that failed (so use bench$index)
-			bench1t = TimeMean(bench[[ bench$index[1] ]]$data)
-			bench2t = TimeMean(bench[[ bench$index[2] ]]$data)
-			bench1bias = mean(bench1t-obst,na.rm=TRUE)
-			bench2bias = mean(bench2t-obst,na.rm=TRUE)
-			diffrange = c(min(modelt-obst,bench1t-obst,bench2t-obst,na.rm=TRUE),
-				max(modelt-obst,bench1t-obst,bench2t-obst,na.rm=TRUE))
-			metrics[[1]] = list(name='TimeSpaceBias',model_value=modelbias,bench_value=list(bench1=bench1bias,bench2=bench2bias))
-		}else if(bench$howmany == 3){
-			bench1t = TimeMean(bench[[ bench$index[1] ]]$data)
-			bench2t = TimeMean(bench[[ bench$index[2] ]]$data)
-			bench3t = TimeMean(bench[[ bench$index[3] ]]$data)
-			bench1bias = mean(bench1t-obst,na.rm=TRUE)
-			bench2bias = mean(bench2t-obst,na.rm=TRUE)
-			bench3bias = mean(bench3t-obst,na.rm=TRUE)
-			diffrange = c(min(modelt-obst,bench1t-obst,bench2t-obst,bench3t-obst,na.rm=TRUE),
-				max(modelt-obst,bench1t-obst,bench2t-obst,bench3t-obst,na.rm=TRUE))
-			metrics[[1]] = list(name='TimeSpaceBias',model_value=modelbias,
-				bench_value=list(bench1=bench1bias,bench2=bench2bias,bench3=bench3bias))
-		}else{
-			# Just save metric list for model:
-			metrics[[1]] = list(name='TimeSpaceBias',model_value=modelbias)	
-		}
-	}else if(plottype=='TimeSD'){
-		if(bench$howmany == 1){
-			# Get benchmark data, noting there may have been other benchmarks that failed (so use bench$index)
-			bench1t = TimeSD(bench[[ bench$index[1] ]]$data)
-			bench1SDbias = mean(bench1t-obst,na.rm=TRUE)
-			diffrange = c(min(modelt-obst,bench1t-obst,na.rm=TRUE),
-				max(modelt-obst,bench1t-obst,na.rm=TRUE))
-			metrics[[1]] = list(name='AvTimeSDbias',model_value=modelSDbias,bench_value=list(bench1=bench1SDbias))	
-		}else if(bench$howmany == 2){
-			# Get benchmark data, noting there may have been other benchmarks that failed (so use bench$index)
-			bench1t = TimeSD(bench[[ bench$index[1] ]]$data)
-			bench2t = TimeSD(bench[[ bench$index[2] ]]$data)
-			bench1SDbias = mean(bench1t-obst,na.rm=TRUE)
-			bench2SDbias = mean(bench2t-obst,na.rm=TRUE)
-			diffrange = c(min(modelt-obst,bench1t-obst,bench2t-obst,na.rm=TRUE),
-				max(modelt-obst,bench1t-obst,bench2t-obst,na.rm=TRUE))
-			metrics[[1]] = list(name='AvTimeSDbias',model_value=modelSDbias,bench_value=list(bench1=bench1SDbias,bench2=bench2SDbias))
-		}else if(bench$howmany == 3){
-			bench1t = TimeSD(bench[[ bench$index[1] ]]$data)
-			bench2t = TimeSD(bench[[ bench$index[2] ]]$data)
-			bench3t = TimeSD(bench[[ bench$index[3] ]]$data)
-			bench1SDbias = mean(bench1t-obst,na.rm=TRUE)
-			bench2SDbias = mean(bench2t-obst,na.rm=TRUE)
-			bench3SDbias = mean(bench3t-obst,na.rm=TRUE)
-			diffrange = c(min(modelt-obst,bench1t-obst,bench2t-obst,bench3t-obst,na.rm=TRUE),
-				max(modelt-obst,bench1t-obst,bench2t-obst,bench3t-obst,na.rm=TRUE))
-			metrics[[1]] = list(name='AvTimeSDbias',model_value=modelSDbias,
-				bench_value=list(bench1=bench1SDbias,bench2=bench2SDbias,bench3=bench3SDbias))
-		}else{
-			# Just save metric list for model:
-			metrics[[1]] = list(name='AvTimeSDbias',model_value=modelSDbias)	
-		}
-	}
-	# Fetch colour scheme:
-	zcols = ChooseColours(zrange,varname,'positive')
-	diffcols = ChooseColours(diffrange,varname,'difference')
-	
-	# First plot: model	mean
-	title = paste(model$name,' ',longvarname,' ',plottype,sep='')
-	errtext = PlotAus(obs$lon,obs$lat,modelt,mean(modelt,na.rm=T),sd(modelt,na.rm=T),
-		varname,unitstxt,longvarname,zrange,zcols,title,textloc)
-	# Second plot: obs mean
-	title = paste(obs$name,' ',longvarname,' ',plottype,sep='')
-	errtext = PlotAus(obs$lon,obs$lat,obst,mean(obst,na.rm=T),sd(obst,na.rm=T),
-		varname,unitstxt,longvarname,zrange,zcols,title,textloc)
-	# Third plot: difference of model, obs mean
-	title = paste('[',model$name,'-',obs$name,'] ',varname,' ',plottype,sep='')
-	errtext = PlotAus(obs$lon,obs$lat,(modelt-obst),mean((modelt-obst),na.rm=T),
-		sd((modelt-obst),na.rm=T),varname,unitstxt,longvarname,diffrange,diffcols,title,textloc)
-	# Plot benchmarks that exist:
-	if(bench$exist){
-		# Fourth plot: difference bench1, obs 
-		title = paste('[',bench[[ bench$index[1] ]]$name,'-',obs$name,'] ',varname,' ',plottype,sep='')
-		errtext = PlotAus(obs$lon,obs$lat,(bench1t - obst),mean((bench1t-obst),na.rm=T),
-			sd((bench1t-obst),na.rm=T),varname,unitstxt,longvarname,diffrange,diffcols,title,textloc)
-		if(bench$howmany >= 2){
-			# Fifth plot: difference bench2, obs 
-			title = paste('[',bench[[ bench$index[2] ]]$name,'-',obs$name,'] ',varname,' ',plottype,sep='')
-			errtext = PlotAus(obs$lon,obs$lat,(bench2t - obst),mean((bench2t-obst),na.rm=T),
-				sd((bench2t-obst),na.rm=T),varname,unitstxt,longvarname,diffrange,diffcols,title,textloc)
-			if(bench$howmany == 3){	
-				# Fifth plot: difference bench3, obs 
-				title = paste('[',bench[[ bench$index[3] ]]$name,'-',obs$name,'] ',varname,' ',plottype,sep='')
-				errtext = PlotAus(obs$lon,obs$lat,(bench3t - obst),mean((bench3t-obst),na.rm=T),
-					sd((bench3t-obst),na.rm=T),varname,unitstxt,longvarname,diffrange,diffcols,title,textloc)
-			}
-		}	
-	}
-	### Add density insets ###
-	mod_den = density(modelt,na.rm=TRUE) # calculate model density estimate 
-	obs_den = density(obst,na.rm=TRUE) # calculate obs density estimate
-	xrange = DensityXrange(list(mod_den,obs_den),density_cut)
-	# Plot pdfs for model and obs
-	in1 = InsetDensity(density_location[[1]],mod_den,xrange)
-	in2 = InsetDensity(density_location[[2]],obs_den,xrange)
-	moderr_den = density((modelt-obst),na.rm=TRUE) # calculate model error density estimate 
-	xrange = DensityXrange(list(moderr_den),density_cut)
-	if(bench$howmany == 1){
-		bench1err_den = density((bench1t-obst),na.rm=TRUE)
-		# If there's a benchmark, overwrite xrange to include benchmark info:
-		xrange = DensityXrange(list(moderr_den,bench1err_den),density_cut)
-		in4 = InsetDensity(density_location[[4]],bench1err_den,xrange)
-	}else if(bench$howmany == 2){
-		bench1err_den = density((bench1t-obst),na.rm=TRUE)
-		bench2err_den = density((bench2t-obst),na.rm=TRUE)
-		xrange = DensityXrange(list(moderr_den,bench1err_den,bench2err_den),density_cut)
-		in4 = InsetDensity(density_location[[4]],bench1err_den,xrange)	
-		in5 = InsetDensity(density_location[[5]],bench2err_den,xrange)	
-	}else if(bench$howmany == 3){
-		bench1err_den = density((bench1t-obst),na.rm=TRUE)
-		bench2err_den = density((bench2t-obst),na.rm=TRUE)
-		bench3err_den = density((bench3t-obst),na.rm=TRUE)
-		xrange = DensityXrange(list(moderr_den,bench1err_den,bench2err_den,bench3err_den),density_cut)
-		in4 = InsetDensity(density_location[[4]],bench1err_den,xrange)	
-		in5 = InsetDensity(density_location[[5]],bench2err_den,xrange)
-		in6 = InsetDensity(density_location[[6]],bench3err_den,xrange)		
-	}
-	in3 = InsetDensity(density_location[[3]],moderr_den,xrange)
-	result = list(errtext=errtext,err=FALSE,metrics=metrics)
-	return(result)
-}
-SpatialAusRelative = function(model,obs,bench,varname,unitstxt,longvarname,metrics,plottype){
-	errtext = 'ok'
-	metrics = list()
-	density_cut = 1/200
-	# Calculate number of map panels:
-	npanels = bench$howmany + 1
-	# Plot layout:
-	if(npanels == 1){
-		density_location = DensityLocationAus(1)
-		textloc='bottomright'
-	}else if(npanels == 2){
-		par(mfcol=c(1,2) ,mar=c(4,4,3,0.5),oma=c(6,0,5,1),mgp=c(2.5,0.7,0),ps=12,tcl=-0.4)
-		density_location = DensityLocationAus(2)
-		textloc='bottomright'
-	}else if(npanels >= 3){
-		par(mfcol=c(2,2) ,mar=c(3,3,3,0.5),oma=c(0,0,0,1),mgp=c(1.8,0.5,0),ps=15,tcl=-0.4)
-		density_location = DensityLocationAus(4)
-		textloc='bottomright'
-	}
-	if(plottype=='TimeRMSE'){
-		# Calculate time means:
-		modelt = TimeRMSE(obs$data, model$data)
-		modelRMSE = sqrt(mean((model$data - obs$data)^2,na.rm=TRUE))
-		supressunits = FALSE
-	}else if(plottype=='TimeCor'){
-		modelt = TimeCor(obs$data, model$data)
-		modelAvTimeCor = mean(modelt,na.rm=TRUE)
-		modelTimeSpaceCor = cor(as.vector(obs$data),as.vector(model$data))
-		supressunits = TRUE
-	}else{
-		result = list(errtext = paste('Unknown plot type \'',plottype,'\' requested in function SpatialAusRelative.',sep=''),err=TRUE)
-		return(result)
-	}
-	zrange = c(min(modelt,na.rm=TRUE),max(modelt,na.rm=TRUE))
-	if(plottype=='TimeRMSE'){
-		if(bench$howmany == 1){
-			# Get benchmark data, noting there may have been other benchmarks that failed (so use bench$index)
-			bench1t = TimeRMSE(obs$data, bench[[ bench$index[1] ]]$data)
-			bench1RMSE = sqrt(mean((bench[[ bench$index[1] ]]$data - obs$data)^2,na.rm=TRUE))
-			metrics[[1]] = list(name='TimeSpaceRMSE',model_value=modelRMSE,bench_value=list(bench1=bench1RMSE))
-			zrange = c(min(modelt,bench1t,na.rm=TRUE),max(modelt,bench1t,na.rm=TRUE))
-		}else if(bench$howmany == 2){
-			# Get benchmark data, noting there may have been other benchmarks that failed (so use bench$index)
-			bench1t = TimeRMSE(obs$data, bench[[ bench$index[1] ]]$data)
-			bench2t = TimeRMSE(obs$data, bench[[ bench$index[2] ]]$data)
-			bench1RMSE = sqrt(mean((bench[[ bench$index[1] ]]$data - obs$data)^2,na.rm=TRUE))
-			bench2RMSE = sqrt(mean((bench[[ bench$index[2] ]]$data - obs$data)^2,na.rm=TRUE))
-			metrics[[1]] = list(name='TimeSpaceRMSE',model_value=modelRMSE,
-				bench_value=list(bench1=bench1RMSE,bench2=bench2RMSE))
-			zrange = c(min(modelt,bench1t,bench2t,na.rm=TRUE),max(modelt,bench1t,bench2t,na.rm=TRUE))
-		}else if(bench$howmany == 3){
-			bench1t = TimeRMSE(obs$data, bench[[ bench$index[1] ]]$data)
-			bench2t = TimeRMSE(obs$data, bench[[ bench$index[2] ]]$data)
-			bench3t = TimeRMSE(obs$data, bench[[ bench$index[3] ]]$data)
-			bench1RMSE = sqrt(mean((bench[[ bench$index[1] ]]$data - obs$data)^2,na.rm=TRUE))
-			bench2RMSE = sqrt(mean((bench[[ bench$index[2] ]]$data - obs$data)^2,na.rm=TRUE))
-			bench3RMSE = sqrt(mean((bench[[ bench$index[3] ]]$data - obs$data)^2,na.rm=TRUE))
-			metrics[[1]] = list(name='TimeSpaceRMSE',model_value=modelRMSE,
-				bench_value=list(bench1=bench1RMSE,bench2=bench2RMSE,bench3=bench3RMSE))
-			zrange = c(min(modelt,bench1t,bench2t,bench3t,na.rm=TRUE),max(modelt,bench1t,bench2t,bench3t,na.rm=TRUE))
-		}else{
-			# Just save metric list for model:
-			metrics[[1]] = list(name='TimeSpaceRMSE',model_value=modelRMSE)
-		}
-	}else if(plottype=='TimeCor'){
-		if(bench$howmany == 1){
-			# Get benchmark data, noting there may have been other benchmarks that failed (so use bench$index)
-			bench1t = TimeCor(obs$data, bench[[ bench$index[1] ]]$data)
-			bench1AvTimeCor = mean(bench1t,na.rm=TRUE)
-			bench1TimeSpaceCor = cor(as.vector(obs$data),as.vector(bench[[ bench$index[1] ]]$data))
-			metrics[[1]] = list(name='AvTimeCor',model_value=modelAvTimeCor,bench_value=list(bench1=bench1AvTimeCor))
-			metrics[[2]] = list(name='TimeSpaceCor',model_value=modelTimeSpaceCor,bench_value=list(bench1=bench1TimeSpaceCor))
-			zrange = c(min(modelt,bench1t,na.rm=TRUE),max(modelt,bench1t,na.rm=TRUE))
-		}else if(bench$howmany == 2){
-			# Get benchmark data, noting there may have been other benchmarks that failed (so use bench$index)
-			bench1t = TimeCor(obs$data, bench[[ bench$index[1] ]]$data)
-			bench2t = TimeCor(obs$data, bench[[ bench$index[2] ]]$data)
-			bench1AvTimeCor = mean(bench1t,na.rm=TRUE)
-			bench1TimeSpaceCor = cor(as.vector(obs$data),as.vector(bench[[ bench$index[1] ]]$data))
-			bench2AvTimeCor = mean(bench2t,na.rm=TRUE)
-			bench2TimeSpaceCor = cor(as.vector(obs$data),as.vector(bench[[ bench$index[2] ]]$data))
-			metrics[[1]] = list(name='AvTimeCor',model_value=modelAvTimeCor,
-				bench_value=list(bench1=bench1AvTimeCor,bench2=bench2AvTimeCor))
-			metrics[[2]] = list(name='TimeSpaceCor',model_value=modelTimeSpaceCor,
-				bench_value=list(bench1=bench1TimeSpaceCor,bench2=bench2TimeSpaceCor))
-			zrange = c(min(modelt,bench1t,bench2t,na.rm=TRUE),max(modelt,bench1t,bench2t,na.rm=TRUE))
-		}else if(bench$howmany == 3){
-			bench1t = TimeCor(obs$data, bench[[ bench$index[1] ]]$data)
-			bench2t = TimeCor(obs$data, bench[[ bench$index[2] ]]$data)
-			bench3t = TimeCor(obs$data, bench[[ bench$index[3] ]]$data)
-			bench1AvTimeCor = mean(bench1t,na.rm=TRUE)
-			bench1TimeSpaceCor = cor(as.vector(obs$data),as.vector(bench[[ bench$index[1] ]]$data))
-			bench2AvTimeCor = mean(bench2t,na.rm=TRUE)
-			bench2TimeSpaceCor = cor(as.vector(obs$data),as.vector(bench[[ bench$index[2] ]]$data))
-			bench3AvTimeCor = mean(bench3t,na.rm=TRUE)
-			bench3TimeSpaceCor = cor(as.vector(obs$data),as.vector(bench[[ bench$index[3] ]]$data))
-			metrics[[1]] = list(name='AvTimeCor',model_value=modelAvTimeCor,
-				bench_value=list(bench1=bench1AvTimeCor,bench2=bench2AvTimeCor,bench3=bench3AvTimeCor))
-			metrics[[2]] = list(name='TimeSpaceCor',model_value=modelTimeSpaceCor,
-				bench_value=list(bench1=bench1TimeSpaceCor,bench2=bench2TimeSpaceCor,bench3=bench3TimeSpaceCor))
-			zrange = c(min(modelt,bench1t,bench2t,bench3t,na.rm=TRUE),max(modelt,bench1t,bench2t,bench3t,na.rm=TRUE))
-		}else{
-			# Just save metric list for model:
-			metrics[[1]] = list(name='AvTimeCor',model_value=modelAvTimeCor)
-			metrics[[2]] = list(name='TimeSpaceCor',model_value=modelTimeSpaceCor)
-		}
-	}
-	# Fetch colour scheme:
-	zcols = ChooseColours(zrange,varname,'positive')
-	
-	# First plot: model
-	title = paste(model$name,' ',longvarname,' ',plottype,sep='')
-	errtext = PlotAus(obs$lon,obs$lat,modelt,mean(modelt,na.rm=T),sd(modelt,na.rm=T),
-		varname,unitstxt,longvarname,zrange,zcols,title,textloc,supressunits)
-	# Plot benchmarks that exist:
-	if(bench$exist){
-		# Second plot: bench1
-		title = paste(bench[[ bench$index[1] ]]$name,' ',longvarname,' ',plottype,sep='')
-		errtext = PlotAus(obs$lon,obs$lat,bench1t,mean(bench1t,na.rm=T),
-			sd(bench1t,na.rm=T),varname,unitstxt,longvarname,zrange,zcols,title,textloc,supressunits)
-		if(bench$howmany >= 2){
-			# Third plot: bench2
-			title = paste(bench[[ bench$index[2] ]]$name,' ',longvarname,' ',plottype,sep='')
-			errtext = PlotAus(obs$lon,obs$lat,bench2t,mean(bench2t,na.rm=T),
-				sd(bench2t,na.rm=T),varname,unitstxt,longvarname,zrange,zcols,title,textloc,supressunits)
-			if(bench$howmany == 3){	
-				# Fourth plot: bench3
-				title = paste(bench[[ bench$index[3] ]]$name,' ',longvarname,' ',plottype,sep='')
-				errtext = PlotAus(obs$lon,obs$lat,bench3t,mean(bench3t,na.rm=T),
-					sd(bench3t,na.rm=T),varname,unitstxt,longvarname,zrange,zcols,title,textloc,supressunits)
-			}
-		}	
-	}
-	### Add density insets ###
-	mod_den = density(modelt,na.rm=TRUE) # calculate model density estimate 
-	xrange = DensityXrange(list(mod_den),density_cut) # truncate x-axis range
-	if(bench$howmany == 1){
-		bench1_den = density(bench1t,na.rm=TRUE)
-		# If there's a benchmark, overwrite xrange to include benchmark info:
-		xrange = DensityXrange(list(mod_den,bench1_den),density_cut)
-		in2 = InsetDensity(density_location[[2]],bench1_den,xrange)
-	}else if(bench$howmany == 2){
-		bench1_den = density(bench1t,na.rm=TRUE)
-		bench2_den = density(bench2t,na.rm=TRUE)
-		xrange = DensityXrange(list(mod_den,bench1_den,bench2_den),density_cut)
-		in2 = InsetDensity(density_location[[2]],bench1_den,xrange)	
-		in3 = InsetDensity(density_location[[3]],bench2_den,xrange)	
-	}else if(bench$howmany == 3){
-		bench1_den = density(bench1t,na.rm=TRUE)
-		bench2_den = density(bench2t,na.rm=TRUE)
-		bench3_den = density(bench3t,na.rm=TRUE)
-		xrange = DensityXrange(list(mod_den,bench1_den,bench2_den,bench3_den),density_cut)
-		in2 = InsetDensity(density_location[[2]],bench1_den,xrange)	
-		in3 = InsetDensity(density_location[[3]],bench2_den,xrange)
-		in4 = InsetDensity(density_location[[4]],bench3_den,xrange)		
-	}
-	# Plot pdfs for model
-	in1 = InsetDensity(density_location[[1]],mod_den,xrange)
-	result = list(errtext=errtext,err=FALSE,metrics=metrics)
-	return(result)
-}
-
-TimeMean = function(threedvar){
-	# Take the time mean of 3D variable
-	twodvar = apply(threedvar,c(1,2),mean)	
-	return(twodvar)
-}
-TimeSD = function(threedvar){
-	# Take the time sd of 3D variable
-	twodvar = apply(threedvar,c(1,2),sd)	
-	return(twodvar)
-}
-TimeRMSE = function(obs3d,model3d){
-	 twodvar = apply((model3d - obs3d),c(1,2),rootmeansquare)
-	 return(twodvar)
-}
-rootmeansquare = function(diffvector){
-	result = sqrt(mean(diffvector^2))
-	return(result)
-}
-TimeCor = function(obs3d,model3d){
-#		omean = TimeMean(obs3d)
-#		mmean = TimeMean(model3d)
-#		scov = apply(((obs3d - omean)*(model3d-mmean)),c(1,2),sum)
-#		twodcor = scov / (TimeSD(obs3d) * TimeSD(model3d))
-		twodcor = matrix(NA,length(obs3d[,1,1]),length(obs3d[1,,1]))
-		for(i in 1:length(obs3d[,1,1])){
-			for(j in 1:length(obs3d[1,,1])){
-				twodcor[i,j] = cor(obs3d[i,j,],model3d[i,j,])
-			}
-		}				
-	 return(twodcor)
-}
-
-PlotAus = function(lon,lat,data,meanval,sdval,varname,unitstxt,longvarname,zrange,zcols,title,textloc,supressunits=FALSE){
-	# Generates a gridded heatmap style plot for Australia with map, based on input lat/lon.
-	library(maps)
-	library(mapdata)
-	library(fields) # for image.plot
-	errtext = 'ok'
-	# Decide location of plot for text placement:
-	if(textloc=='bottomright'){
-		textloc1 = c((obs$lon[1] + (obs$lon[length(obs$lon)] - obs$lon[1])*0.9), (obs$lat[1] + (obs$lat[length(obs$lat)] - obs$lat[1])*0.19) )
-		textloc2 = c((obs$lon[1] + (obs$lon[length(obs$lon)] - obs$lon[1])*0.9), (obs$lat[1] + (obs$lat[length(obs$lat)] - obs$lat[1])*0.13) )
-		textloc3 = c((obs$lon[1] + (obs$lon[length(obs$lon)] - obs$lon[1])*0.9), (obs$lat[1] + (obs$lat[length(obs$lat)] - obs$lat[1])*0.07) )
-	}else if(textloc=='topleft'){
-		textloc1 = c((obs$lon[1] + (obs$lon[length(obs$lon)] - obs$lon[1])*0.15), (obs$lat[1] + (obs$lat[length(obs$lat)] - obs$lat[1])*0.94) )
-		textloc2 = c((obs$lon[1] + (obs$lon[length(obs$lon)] - obs$lon[1])*0.15), (obs$lat[1] + (obs$lat[length(obs$lat)] - obs$lat[1])*0.88) )
-		textloc3 = c((obs$lon[1] + (obs$lon[length(obs$lon)] - obs$lon[1])*0.15), (obs$lat[1] + (obs$lat[length(obs$lat)] - obs$lat[1])*0.82) )
-	}
-	# Plot:	
-	image.plot(lon,lat,data,xlab='Longitude',ylab='Latitude',col=zcols,zlim=zrange,legend.mar=5.5)
-	map('worldHires',add=TRUE,wrap=TRUE,xlim=c(min(lon),max(lon)),ylim=c(min(lat),max(lat))) # Add map
-	title(title) 
-	if(!supressunits){
-		text(x=textloc1[1],y=textloc1[2],labels=unitstxt)
-	}
-	text(x=textloc2[1],y=textloc2[2],labels=paste('Mean:',signif(meanval,3)))
-	text(x=textloc3[1],y=textloc3[2],labels=paste('SD:',signif(sdval,2)))
-	return(errtext)
-}
-
-DensityLocationAus = function(npanels){
-	# Simply prescribes locations for density inset on gridded Australia map plots
-	density_location = list()
-	if(npanels == 1){
-		density_location = list()
-		density_location[[1]] = c(0.1,0.5,0.16,0.35)
-	}else if(npanels == 2){
-		density_location = list()
-		density_location[[1]] = c(0.07,0.3,0.17,0.37)
-		density_location[[2]] = c(0.57,0.8,0.17,0.37)
-	}else if(npanels == 4){
-		density_location = list()
-		density_location[[1]] = c(0.06,0.3,0.59,0.68)
-		density_location[[2]] = c(0.56,0.8,0.59,0.68)
-		density_location[[3]] = c(0.06,0.3,0.09,0.18)
-		density_location[[4]] = c(0.56,0.8,0.09,0.18)
-	}else if(npanels==6){
-		density_location = list()
-		density_location[[1]] = c(0.04,0.2,0.58,0.67)
-		density_location[[2]] = c(0.38,0.54,0.58,0.67)
-		density_location[[3]] = c(0.72,0.87,0.58,0.67)
-		density_location[[4]] = c(0.04,0.2,0.08,0.17)
-		density_location[[5]] = c(0.38,0.54,0.08,0.17)
-		density_location[[6]] = c(0.72,0.87,0.08,0.17)
-	}
-	return(density_location)
-}
-
-InsetDensity = function(location,densitydata,xrange){
-	# Adds an inset density plot
-	par(fig=location,new=T)
-	plot(densitydata,lwd=3,main='',ylab='',xlab='',cex.axis=0.8,bty='n',mgp=c(2,0,0),yaxt='n',xlim=xrange,tcl=-0.2)
-}
-
-DensityXrange = function(density_list,density_cut){
-	# Finds the x-axis range that contains all y values above a threshold
-	# for a list of density functions.
-	ymax = 0 # initialise
-	xmin=NA # initialise
-	xmax=NA # initialise
-	for(d in 1:length(density_list)){
-		ymax = max(ymax,density_list[[d]][[2]])	
-	}
-	# Determine where to truncate x-axis according to density cut threshold:
-	for(d in 1:length(density_list)){
-		xmin = min(xmin, density_list[[d]][[1]][ (density_list[[d]][[2]]>(ymax*density_cut)) ], na.rm=TRUE)
-		xmax = max(xmax, density_list[[d]][[1]][ (density_list[[d]][[2]]>(ymax*density_cut)) ], na.rm=TRUE)
-	}
-	return(c(xmin,xmax))
-}
-
-
-
-
diff --git a/pals/R/2DGlobalMean.R b/pals/R/2DGlobalMean.R
deleted file mode 100644
index eb6fb69..0000000
--- a/pals/R/2DGlobalMean.R
+++ /dev/null
@@ -1,18 +0,0 @@
-# 2DGlobalMean.R
-#
-# Gab Abramowitz, UNSW, 2014, gabsun at gmail dot com
-#
-SpatialGlobalMean = function(lon,lat,moutput,obs=NULL,bench1=NULL,bench2=NULL,bench3=NULL){
-	#library(maps)
-	#library(mapdata)
-	
-	# calculate number of map panels:
-	npanels = 1 + as.numeric((! is.null(obs))) + as.numeric((! is.null(bench1))) +
-		as.numeric((! is.null(bench2))) + as.numeric((! is.null(bench3)))
-	
-	image.plot(lon,lat,meanplot,col=allcols,
-	xlab='Longitude',ylab='Latitude',zlim=zrange)
-	map('worldHires',add=TRUE,wrap=TRUE,xlim=c(-180,180),ylim=c(-90,90))
-	title(paste('(a) dT 2080-99 vs 1980-99 - CMIP5 mean',rcp))
-
-}
\ No newline at end of file
diff --git a/pals/R/2DGlobalPacific.R b/pals/R/2DGlobalPacific.R
new file mode 100644
index 0000000..79df481
--- /dev/null
+++ b/pals/R/2DGlobalPacific.R
@@ -0,0 +1,91 @@
+# 2DGlobal.R
+#
+# Gab Abramowitz, UNSW, 2015, gabsun at gmail dot com
+#
+
+PlotGlobal = function(lon,lat,data,meanval,sdval,varname,unitstxt,longvarname,zrange,zcols,title,textloc,suppressunits=FALSE){
+	# Generates a gridded heatmap style plot for Australia with map, based on input lat/lon.
+	library(maps)
+	library(mapdata)
+	library(fields) # for image.plot
+	errtext = 'ok'
+	# Decide location of plot for text placement:
+	if(textloc=='bottomright'){
+		textloc1 = c((lon[1] + (lon[length(lon)] - lon[1])*0.92), (lat[1] + (lat[length(lat)] - lat[1])*0.30) )
+		textloc2 = c((lon[1] + (lon[length(lon)] - lon[1])*0.92), (lat[1] + (lat[length(lat)] - lat[1])*0.24) )
+		textloc3 = c((lon[1] + (lon[length(lon)] - lon[1])*0.92), (lat[1] + (lat[length(lat)] - lat[1])*0.18) )
+	}else if(textloc=='middle'){
+		textloc1 = c((lon[1] + (lon[length(lon)] - lon[1])*0.53), (lat[1] + (lat[length(lat)] - lat[1])*0.70) )
+		textloc2 = c((lon[1] + (lon[length(lon)] - lon[1])*0.53), (lat[1] + (lat[length(lat)] - lat[1])*0.64) )
+		textloc3 = c((lon[1] + (lon[length(lon)] - lon[1])*0.53), (lat[1] + (lat[length(lat)] - lat[1])*0.58) )
+	}
+	
+	cat('TTT pre plot:',proc.time()[3],'\n')
+	
+	# Set x, y limits to use up more of the plot region:
+	fact = 0.02
+	xmin = min(lon)
+	xmax = max(lon)
+	ymin = min(lat)
+	ymax = max(lat)
+	xrange = c(xmin+fact*(xmax-xmin),xmax-fact*(xmax-xmin))
+	yrange = c(ymin+fact*(ymax-ymin),ymax-fact*(ymax-ymin))
+	
+	# Draw image:	
+	image.plot(lon,lat,data,xlab='Longitude',ylab='Latitude',col=zcols,xlim=xrange,
+		ylim=yrange,zlim=zrange,legend.mar=5.5)
+	
+	cat('TTT pre map:',proc.time()[3],'\n')
+	
+	map('world2Hires',add=TRUE,wrap=TRUE,xlim=c(min(lon),max(lon)),ylim=c(min(lat),max(lat))) # Add map
+	title(title) 
+	if(!suppressunits){
+		text(x=textloc1[1],y=textloc1[2],labels=unitstxt)
+	}
+	text(x=textloc2[1],y=textloc2[2],labels=paste('Mean:',signif(meanval,3)))
+	text(x=textloc3[1],y=textloc3[2],labels=paste('SD:',signif(sdval,2)))
+	
+	cat('TTT post all plotting:',proc.time()[3],'\n')
+	
+	return(errtext)
+}
+
+TextLocationGlobal = function(npanels){
+	if(npanels<5){
+		text_location = 'bottomright'
+	}else{
+		text_location = 'middle'
+	}
+}
+
+DensityLocationGlobal = function(npanels){
+	# Simply prescribes locations for density inset on gridded Australia map plots
+	density_location = list()
+	if(npanels == 1){
+		density_location = list()
+		density_location[[1]] = c(0.45,0.7,0.25,0.4)
+	}else if(npanels == 2){
+		density_location = list()
+		density_location[[1]] = c(0.24,0.37,0.27,0.37)
+		density_location[[2]] = c(0.74,0.86,0.27,0.37)
+	}else if(npanels == 4){
+		density_location = list()
+		density_location[[1]] = c(0.24,0.37,0.63,0.71)
+		density_location[[2]] = c(0.74,0.86,0.63,0.71)
+		density_location[[3]] = c(0.24,0.37,0.13,0.21)
+		density_location[[4]] = c(0.74,0.86,0.13,0.21)
+	}else if(npanels==6){
+		density_location = list()
+		density_location[[1]] = c(0.155,0.245,0.62,0.7)
+		density_location[[2]] = c(0.49,0.58,0.62,0.7)
+		density_location[[3]] = c(0.825,0.915,0.62,0.7)
+		density_location[[4]] = c(0.155,0.245,0.12,0.2)
+		density_location[[5]] = c(0.49,0.58,0.12,0.2)
+		density_location[[6]] = c(0.825,0.915,0.12,0.2)
+	}
+	return(density_location)
+}
+
+
+
+
diff --git a/pals/R/2DMetrics.R b/pals/R/2DMetrics.R
new file mode 100644
index 0000000..537e49e
--- /dev/null
+++ b/pals/R/2DMetrics.R
@@ -0,0 +1,177 @@
+# 2DMetrics.R
+#
+# Gab Abramowitz, UNSW, 2015, gabsun at gmail dot com
+#
+TimeMeanAll = function(model,obs,bench,variable,plottype,cl){
+	# Calculates time means for each grid point in obs, model and benchmark data  
+	metrics = list()
+	benchm = list()
+	benchbias = c()
+	# Calculate time means for plotting:
+	modelm = TimeMean(model$data,cl)
+	obsm = TimeMean(obs$data,cl) + modelm - modelm # to make sure ocean areas not included
+	# Ranges of obs, model metric values (benchmarks only appear in difference plots):
+	zrange = c(min(modelm,obsm,na.rm=TRUE),max(modelm,obsm,na.rm=TRUE))
+	# Scalar metric for reporting:
+	modelbias = mean(modelm-obsm,na.rm=TRUE)
+	# Initial ranges for difference plots:
+	dmax = max(modelm-obsm,na.rm=TRUE)
+	dmin = min(modelm-obsm,na.rm=TRUE)
+	if(bench$exist){
+		for(b in 1:bench$howmany){
+			# Get benchmark metric data, noting there may have been other benchmarks 
+			# that failed (so use bench$index)
+			benchm[[b]] = TimeMean(bench[[ bench$index[b] ]]$data,cl)
+			benchbias[b] = mean(benchm[[b]]-obsm,na.rm=TRUE)
+			dmax = max(dmax,(benchm[[b]]-obsm),na.rm=TRUE)
+			dmin = min(dmin,(benchm[[b]]-obsm),na.rm=TRUE)
+		}
+	}
+	metrics[[1]] = list(name='TimeSpaceBias',model_value=modelbias,bench_value=benchbias)	
+	result = list(modelm = modelm, obsm=obsm, benchm=benchm, diffrange = c(dmin, dmax), 
+		zrange=zrange,metrics=metrics)
+	return(result)
+}
+
+TimeSDAll = function(model,obs,bench,variable,plottype,cl){
+	# Calculates standard deviation for each grid point in obs, model and benchmark data
+	metrics = list()
+	benchm = list()
+	benchSDbias = c()
+	# Calculate time means:
+	modelm = TimeSD(model$data,cl)
+	obsm = TimeSD(obs$data,cl) + modelm - modelm # to make sure ocean areas not included
+	# Ranges of obs, model metric values (benchmarks only appear in difference plots):
+	zrange = c(min(modelm,obsm,na.rm=TRUE),max(modelm,obsm,na.rm=TRUE))
+	# Scalar metric for reporting:
+	modelSDbias = mean(modelm-obsm,na.rm=TRUE)
+	# Initial ranges for difference plots:
+	dmax = max(modelm-obsm,na.rm=TRUE)
+	dmin = min(modelm-obsm,na.rm=TRUE)
+	if(bench$exist){
+		for(b in 1:bench$howmany){
+			# Get benchmark metric data, noting there may have been other benchmarks 
+			# that failed (so use bench$index)
+			benchm[[b]] = TimeSD(bench[[ bench$index[b] ]]$data,cl)
+			benchSDbias[b] = mean(benchm[[b]]-obsm,na.rm=TRUE)
+			dmax = max(dmax,(benchm[[b]]-obsm),na.rm=TRUE)
+			dmin = min(dmin,(benchm[[b]]-obsm),na.rm=TRUE)
+		}
+	}
+	metrics[[1]] = list(name='AvTimeSDbias',model_value=modelSDbias,bench_value=benchSDbias)	
+	result = list(modelm = modelm, obsm=obsm, benchm=benchm, diffrange = c(dmin, dmax), 
+		zrange=zrange,metrics=metrics)
+	return(result)
+}
+
+TimeRMSEAll = function(model,obs,bench,variable,plottype,cl){
+	# Calculates root mean square error for each grid point for model and benchmark data
+	metrics = list()
+	benchm = list()
+	benchRMSE = c()
+	suppressunits = FALSE # i.e. RMSE has units - unlike, e.g. correlation
+	# Calculate time RMSE for plotting:
+	modelm = TimeRMSE(obs$data, model$data,cl)
+	modelRMSE = sqrt(mean((model$data - obs$data)^2,na.rm=TRUE)) # scalar reporting metric
+	# Initial ranges:
+	rmax = max(modelm,na.rm=TRUE)
+	rmin = min(modelm,na.rm=TRUE)
+	if(bench$exist){
+		for(b in 1:bench$howmany){
+			# Get benchmark metric data, noting there may have been other benchmarks 
+			# that failed (so use bench$index)
+			benchm[[b]] = TimeRMSE(obs$data, bench[[ bench$index[b] ]]$data,cl)
+			benchRMSE[b] = sqrt(mean((bench[[ bench$index[b] ]]$data - obs$data)^2,na.rm=TRUE))
+			rmax = max(rmax,benchm[[b]],na.rm=TRUE)
+			rmin = min(rmin,benchm[[b]],na.rm=TRUE)
+		}
+	}
+	metrics[[1]] = list(name='TimeSpaceRMSE',model_value=modelRMSE,bench_value=benchRMSE)	
+	result = list(modelm = modelm, benchm=benchm, zrange = c(rmin, rmax), 
+		metrics=metrics, suppressunits = suppressunits)
+	return(result)
+}
+
+TimeCorAll = function(model,obs,bench,variable,plottype,cl){
+	# Calculates correlation for each grid point for model,obs and benchmark,obs data
+	metrics = list()  
+	benchm = list()
+	benchAvTimeCor = c()
+	benchTimeSpaceCor = c()
+	suppressunits = TRUE # i.e. correlation has no units
+	# Calculate time correlation for plotting:
+	modelm = TimeCor(obs$data, model$data,cl)
+	# Two scalar metrics:
+	modelAvTimeCor = mean(modelm,na.rm=TRUE) # Average of time correlation
+	modelTimeSpaceCor = cor(as.vector(obs$data),as.vector(model$data)) # Cor over time and space
+	# initial ranges:
+	rmax = max(modelm,na.rm=TRUE)
+	rmin = min(modelm,na.rm=TRUE)
+	if(bench$exist){
+		for(b in 1:bench$howmany){
+			# Get benchmark metric data, noting there may have been other benchmarks 
+			# that failed (so use bench$index)
+			benchm[[b]] = TimeCor(obs$data, bench[[ bench$index[b] ]]$data,cl)
+			benchAvTimeCor[b] = mean(benchm[[b]],na.rm=TRUE)
+			benchTimeSpaceCor[b] = cor(as.vector(obs$data),as.vector(bench[[ bench$index[b] ]]$data))
+			rmax = max(rmax,benchm[[b]],na.rm=TRUE)
+			rmin = min(rmin,benchm[[b]],na.rm=TRUE)
+		}
+	}
+	metrics[[1]] = list(name='AvTimeCor',model_value=modelAvTimeCor,bench_value=benchAvTimeCor)
+	metrics[[2]] = list(name='TimeSpaceCor',model_value=modelTimeSpaceCor,bench_value=benchTimeSpaceCor)
+	result = list(modelm = modelm, benchm=benchm, zrange = c(rmin, rmax), 
+		metrics=metrics, suppressunits = suppressunits)
+	return(result)
+}
+
+TimeMean = function(threedvar,cl){
+	# Take the time mean of 3D variable
+	if(is.null(cl)){
+		twodvar = apply(threedvar,c(1,2),mean)	
+	}else{
+		twodvar = parApply(cl,threedvar,c(1,2),mean)
+	}
+	return(twodvar)
+}
+
+TimeSD = function(threedvar, cl){
+	# Take the time sd of 3D variable
+	if(is.null(cl)){
+		twodvar = apply(threedvar,c(1,2),sd)
+	}else{
+		twodvar = parApply(cl,threedvar,c(1,2),sd)	
+	}
+	return(twodvar)
+}
+
+TimeRMSE = function(obs3d,model3d,cl){
+	if(is.null(cl)){
+		twodvar = apply((model3d - obs3d),c(1,2),rootmeansquare)
+	}else{
+		twodvar = parApply(cl,(model3d - obs3d),c(1,2),rootmeansquare)
+	}
+	return(twodvar)
+}
+
+rootmeansquare = function(diffvector){
+	result = sqrt(mean(diffvector^2))
+	return(result)
+}
+
+TimeCor = function(obs3d,model3d,cl){
+	spacedim = dim(obs3d[,,1])
+	indx = array(NA,dim=c(spacedim,2))
+	indx[,,1] = matrix(1:spacedim[1],nrow=spacedim[1],ncol=spacedim[2])
+	indx[,,2] = matrix(1:spacedim[2],nrow=spacedim[1],ncol=spacedim[2],byrow=TRUE)
+	if(is.null(cl)){
+		twodcor = apply(indx,c(1,2),ApplyCor,obs3d,model3d)
+	}else{
+		twodcor = parApply(cl,indx,c(1,2),ApplyCor,obs3d,model3d)
+	}
+	return(twodcor)
+}
+
+ApplyCor = function(index,obs3d,model3d){
+	scalarcor = cor(obs3d[index[1],index[2],],model3d[index[1],index[2],])
+}
diff --git a/pals/R/2DPlotLayout.R b/pals/R/2DPlotLayout.R
new file mode 100644
index 0000000..8312d1f
--- /dev/null
+++ b/pals/R/2DPlotLayout.R
@@ -0,0 +1,211 @@
+# 2DPlotLayout.R
+#
+# Gab Abramowitz, UNSW, 2015, gabsun at gmail dot com
+#
+SpatialPlotAbsolute = function(model,obs,bench,md,variable,plottype,region){
+	# Layout of plots that show obs panel, model panel, model-obs panel and 
+	# up to 3 benchmark-obs panels (i.e. max 6 panels).
+	# Density plots are always in the lower left, mean and SD values are either bottom
+	# right <= 4 panels, or top left for > 4 panels
+	errtext = 'ok'
+	varname = variable[['Name']][1]
+	unitstxt = variable[['UnitsText']]
+	longvarname = variable[['PlotName']]
+	density_cut = 1/200 # truncate denisty plot x-axis at this fraction of max y value 
+	npanels = bench$howmany + 3 # Number of map panels in plot
+	
+	# Plot layout:
+	if(npanels <= 4){
+		par(mfcol=c(2,2) ,mar=c(3,3,3,0.5),oma=c(0,0,0,1),mgp=c(1.8,0.5,0),ps=15,tcl=-0.4)
+		density_location = DensityLocation(region,4)
+		textloc = TextLocation(region,4)
+	}else{
+		par(mfcol=c(2,3) ,mar=c(3,3,3,0.5),oma=c(1,0,0.5,1),mgp=c(1.8,0.5,0),ps=18,tcl=-0.2)
+		density_location = DensityLocation(region,6)
+		textloc = TextLocation(region,6)
+	}
+	
+	# Fetch colour scheme:
+	zcols = ChooseColours(md$zrange,varname,'positive')
+	diffcols = ChooseColours(md$diffrange,varname,'difference')
+	
+	### Draw plot panels ###
+	# First plot: model	
+	title = paste(model$name,' ',longvarname,' ',plottype,sep='')
+	errtext = DrawPlot(region,obs$grid$lon,obs$grid$lat,md$modelm,mean(md$modelm,na.rm=T),sd(md$modelm,na.rm=T),
+		varname,unitstxt,longvarname,md$zrange,zcols,title,textloc)
+	# Second plot: obs 
+	title = paste(obs$name,' ',longvarname,' ',plottype,sep='')
+	errtext = DrawPlot(region,obs$grid$lon,obs$grid$lat,md$obsm,mean(md$obsm,na.rm=T),sd(md$obsm,na.rm=T),
+		varname,unitstxt,longvarname,md$zrange,zcols,title,textloc)
+	# Third plot: difference of model, obs
+	title = paste('[',model$name,'-',obs$name,'] ',varname,' ',plottype,sep='')
+	errtext = DrawPlot(region,obs$grid$lon,obs$grid$lat,(md$modelm-md$obsm),mean((md$modelm-md$obsm),na.rm=T),
+		sd((md$modelm-md$obsm),na.rm=T),varname,unitstxt,longvarname,md$diffrange,diffcols,title,textloc)
+	# Plot benchmark obs differences that exist (plots 4-6):
+	if(bench$exist){
+		for(b in 1:bench$howmany){
+			title = paste('[',bench[[ bench$index[b] ]]$name,'-',obs$name,'] ',varname,' ',plottype,sep='')
+			errtext = DrawPlot(region,obs$grid$lon,obs$grid$lat,(md$benchm[[b]] - md$obsm),mean((md$benchm[[b]]-md$obsm),na.rm=T),
+				sd((md$benchm[[b]]-md$obsm),na.rm=T),varname,unitstxt,longvarname,md$diffrange,diffcols,title,textloc)
+		}
+	}	
+		
+	### Add density insets ###
+	mod_den = density(md$modelm,na.rm=TRUE) # calculate model density estimate 
+	obs_den = density(md$obsm,na.rm=TRUE) # calculate obs density estimate
+	xrange = DensityXrange(list(mod_den,obs_den),density_cut)
+	# Plot pdfs for model and obs
+	in1 = InsetDensity(density_location[[1]],mod_den,xrange)
+	in2 = InsetDensity(density_location[[2]],obs_den,xrange)
+	moderr_den = density((md$modelm-md$obsm),na.rm=TRUE) # calculate model error density estimate 
+	xrange = DensityXrange(list(moderr_den),density_cut)
+	if(bench$exist){
+		bencherr_den = list()
+		density_range_list = list(moderr_den)
+		for(b in 1:bench$howmany){
+			bencherr_den[[b]] = density((md$benchm[[b]]-md$obsm),na.rm=TRUE)
+			density_range_list[[b+1]] = bencherr_den[[b]]
+		}
+		# If there's a benchmark, overwrite xrange to include benchmark info:
+		xrange = DensityXrange(density_range_list,density_cut)
+		for(b in 1:bench$howmany){
+			inb = InsetDensity(density_location[[(b+3)]],bencherr_den[[b]],xrange)
+		}
+	}
+	in3 = InsetDensity(density_location[[3]],moderr_den,xrange)
+	
+	result = list(errtext=errtext,err=FALSE,metrics=md$metrics)
+	return(result)
+}
+
+SpatialPlotRelative = function(model,obs,bench,md,variable,plottype,region){
+	# Layout of plots that show metrics that use model and obs together, e.g. RMSE model panel 
+	# and RMSE panels for up to 3 benchmark-obs (i.e. max 4 panels).
+	# Density plots are always in the lower left, mean and SD values are on the bottom right
+	errtext = 'ok'
+	varname = variable[['Name']][1]
+	unitstxt = variable[['UnitsText']]
+	longvarname = variable[['PlotName']]
+	density_cut = 1/200 # truncate denisty plot x-axis at this fraction of max y value
+	npanels = bench$howmany + 1 # Number of map panels
+	
+	# Plot layout:
+	if(npanels == 1){
+		density_location = DensityLocation(region,1)
+		textloc=TextLocation(region,1)
+	}else if(npanels == 2){
+		par(mfcol=c(1,2) ,mar=c(4,4,3,0.5),oma=c(6,0,5,1),mgp=c(2.5,0.7,0),ps=12,tcl=-0.4)
+		density_location = DensityLocation(region,2)
+		textloc=TextLocation(region,2)
+	}else if(npanels >= 3){
+		par(mfcol=c(2,2) ,mar=c(3,3,3,0.5),oma=c(0,0,0,1),mgp=c(1.8,0.5,0),ps=15,tcl=-0.4)
+		density_location = DensityLocation(region,4)
+		textloc=TextLocation(region,4)
+	}
+	
+	# Fetch colour scheme:
+	zcols = ChooseColours(md$zrange,varname,'positive')
+	
+	### Draw plot panels ###
+	# First plot: model
+	title = paste(model$name,' ',longvarname,' ',plottype,sep='')
+	errtext = DrawPlot(region,obs$grid$lon,obs$grid$lat,md$modelm,mean(md$modelm,na.rm=T),sd(md$modelm,na.rm=T),
+		varname,unitstxt,longvarname,md$zrange,zcols,title,textloc,md$suppressunits)
+	# Plot benchmarks that exist:
+	if(bench$exist){
+		for(b in 1:bench$howmany){
+			title = paste(bench[[ bench$index[b] ]]$name,' ',longvarname,' ',plottype,sep='')
+			errtext = DrawPlot(region,obs$grid$lon,obs$grid$lat,md$benchm[[b]],mean(md$benchm[[b]],na.rm=T),
+				sd(md$benchm[[b]],na.rm=T),varname,unitstxt,longvarname,md$zrange,zcols,
+				title,textloc,md$suppressunits)
+		}
+	}	
+		
+	### Add density insets ###
+	mod_den = density(md$modelm,na.rm=TRUE) # calculate model density estimate 
+	xrange = DensityXrange(list(mod_den),density_cut) # truncate x-axis range
+	if(bench$exist){
+		bench_den = list()
+		density_range_list = list(mod_den)
+		for(b in 1:bench$howmany){
+			bench_den[[b]] = density(md$benchm[[b]],na.rm=TRUE)
+			density_range_list[[b+1]] = bench_den[[b]]
+		}
+		# If there's a benchmark, overwrite xrange to include benchmark info:
+		xrange = DensityXrange(density_range_list,density_cut)
+		for(b in 1:bench$howmany){
+			inb = InsetDensity(density_location[[(b+1)]],bench_den[[b]],xrange)
+		}
+	}
+	# Plot density inset for model metric
+	in1 = InsetDensity(density_location[[1]],mod_den,xrange)	
+		
+	result = list(errtext=errtext,err=FALSE,metrics=md$metrics)
+	return(result)
+}
+
+DensityXrange = function(density_list,density_cut){
+	# Finds the x-axis range that contains all y values above a threshold
+	# for a list of density functions.
+	ymax = 0 # initialise
+	xmin=NA # initialise
+	xmax=NA # initialise
+	for(d in 1:length(density_list)){
+		ymax = max(ymax,density_list[[d]][[2]])	
+	}
+	# Determine where to truncate x-axis according to density cut threshold:
+	for(d in 1:length(density_list)){
+		xmin = min(xmin, density_list[[d]][[1]][ (density_list[[d]][[2]]>(ymax*density_cut)) ], na.rm=TRUE)
+		xmax = max(xmax, density_list[[d]][[1]][ (density_list[[d]][[2]]>(ymax*density_cut)) ], na.rm=TRUE)
+	}
+	return(c(xmin,xmax))
+}
+
+InsetDensity = function(location,densitydata,xrange){
+	# Adds an inset density plot
+	par(fig=location,new=T)
+	plot(densitydata,lwd=3,main='',ylab='',xlab='',cex.axis=0.8,bty='n',mgp=c(2,0,0),yaxt='n',xlim=xrange,tcl=-0.2)
+}
+
+DensityLocation = function(region,maxpanels){
+	# Returns location within plot tiles of density functions
+	if(region=='Australia'){
+		result = DensityLocationAus(maxpanels)	
+	}else if(region=='Global'){
+		result = DensityLocationGlobal(maxpanels)
+	}else{
+		stop('region unknown!')	
+	}
+	return(result)
+}
+
+TextLocation = function(region,maxpanels){
+	# Returns location within plot tiles of density functions
+	if(region=='Australia'){
+		result = TextLocationAus(maxpanels)	
+	}else if(region=='Global'){
+		result = TextLocationGlobal(maxpanels)
+	}else{
+		stop('region unknown!')	
+	}
+	return(result)
+}
+
+DrawPlot = function(region,lon,lat,data,meanval,sdval,varname,unitstxt,longvarname,
+	zrange,zcols,title,textloc,suppressunits=FALSE){
+	
+	if(region=='Australia'){
+		result = PlotAus(lon,lat,data,meanval,sdval,varname,unitstxt,longvarname,
+			zrange,zcols,title,textloc,suppressunits)	
+	}else if(region=='Global'){
+		result = PlotGlobal(lon,lat,data,meanval,sdval,varname,unitstxt,longvarname,
+			zrange,zcols,title,textloc,suppressunits)
+	}
+	else{
+		stop('region unknown!')	
+	}
+	return(result)			
+}
+
+
diff --git a/pals/R/Colours.R b/pals/R/Colours.R
new file mode 100644
index 0000000..c038927
--- /dev/null
+++ b/pals/R/Colours.R
@@ -0,0 +1,138 @@
+# Colours.R
+#
+# Functions deciding PALS plot colours
+#
+# Gab Abramowitz UNSW 2014 (palshelp at gmail dot com)
+#
+LineColours = function() {
+	# For line plots:
+	plotcolours=c('black','blue2','indianred3','gold2','yellowgreen')
+}
+LineColours2 = function() {
+	# For line plots:
+	plotcolours2=DesaturateColours(LineColours(),sat=0.5)
+}
+
+TableWinColour = function() {
+	'lightgreen'	
+}
+TableModelColour = function() {
+	'gray87'	
+}
+TableBenchmarkColour = function() {
+	'gray93'	
+}
+BenchmarkColours = function(bench,plotobs=TRUE){
+	# Returns vector of colours for lines in a plot.
+	# If no obs line in the plot (e.g. error plot), first index will be model, else obs
+	# If for some reason a benchmark failed (e.g. missing variable), colours are adjusted to make them 
+	# consistent across different plots 
+	plotcolours = c()
+	if(plotobs){ # i.e. obs line will be part of the plot
+		plotcolours[1] = LineColours()[1]
+	}
+	plotcolours = c(plotcolours, LineColours()[2])
+	if(bench$exist){
+		plotcolours = c(plotcolours, LineColours()[ (bench$index[1] + 2) ])
+		if(bench$howmany == 2){
+			plotcolours = c(plotcolours, LineColours()[ (bench$index[2] + 2) ])
+		}else if(bench$howmany == 3){
+			plotcolours = c(plotcolours, LineColours()[ (bench$index[2] + 2) ])
+			plotcolours = c(plotcolours, LineColours()[ (bench$index[3] + 2) ])
+		}
+	}
+	return(plotcolours)
+}
+ChooseColours = function(range,variablename,plottype,diffthreshold=NULL){
+	# Returns a colour range for gridded plots
+	library(colorRamps)
+	
+	# Full / most range:
+	red2blue = colorRampPalette(c('red','orange','yellow','green','blue'))
+	yellow2purpleCool = colorRampPalette(c('yellow','green3','blue','darkorchid4'))
+	yellow2purpleWarm = colorRampPalette(c('yellow','red','magenta'))
+	purple2yellowWarm = colorRampPalette(c('magenta','red','yellow'))
+	iceblue2green = colorRampPalette(c('slategray1','midnightblue','blue','green3','green'))
+	green2iceblue = colorRampPalette(c('green','green3','blue','midnightblue','slategray1'))
+	
+	# Half range:
+	green2darkblue = colorRampPalette(c('green','green4','blue','midnightblue'))
+	darkblue2green = colorRampPalette(c('midnightblue','blue','green4','green'))
+	darkred2yellow = colorRampPalette(c('red4','red','orange','yellow'))
+	yellow2darkred = colorRampPalette(c('yellow','orange','red','red4'))
+	
+	# Small range:
+	yellow2red = colorRampPalette(c('yellow','red'))
+	red2yellow = colorRampPalette(c('red','yellow'))
+	green2blue = colorRampPalette(c('green','blue'))
+	blue2green = colorRampPalette(c('blue','green'))
+	
+	coolvars = c('Qle','Evap')
+	warmvars = c('Tair','Qh','Rnet','SWdown','SWnet')
+	colourres = 36 # approximately how many colours in a plot (will control size of white space if diff plot)
+	
+	# If no difference threshold has been specified, use 5%:
+	if(is.null(diffthreshold)){
+		diffthreshold = (range[2] - range[1]) / 20
+	}
+	
+	# Assess cases where colours for a difference plot are requested first:
+	if(plottype=='difference'){
+		# i.e. the plot will contain a zero that we want coloured white
+		# First check that we really do need a difference plot:
+		if(range[1] > (-1*diffthreshold)){
+			# Just use a positive scale
+			plottype = 'positive'
+		}else if(range[2]<diffthreshold){
+			# Just use a negative scale
+			plottype = 'negative'
+		}
+		# Find fraction of range below 0
+		lowfrac = abs(range[1]) / (abs(range[1]) + range[2])
+		lownum = floor(lowfrac * colourres)
+		# Find fraction of range above 0
+		highfrac = range[2] / (abs(range[1]) + range[2])
+		highnum = floor(highfrac * colourres)
+		# Decide colour range:
+		if(any(warmvars == variablename)){ # For variables warm colours when positive
+			if(lowfrac/highfrac > 2){ # most fo the range is below 0
+				colours = c(iceblue2green(lownum),'#FFFFFF','#FFFFFF',yellow2red(highnum))
+			}else if(lowfrac/highfrac < 1/2){ # most of the range is above 0
+				colours = c(blue2green(lownum),'#FFFFFF','#FFFFFF',yellow2purpleWarm(highnum))
+			}else{
+				colours = c(darkblue2green(lownum),'#FFFFFF','#FFFFFF',yellow2darkred(highnum))
+			}
+		}else if(any(coolvars == variablename)){ # For variables cool colours when positive
+			if(lowfrac/highfrac > 2){ # most fo the range is below 0
+				colours = c(purple2yellowWarm(lownum),'#FFFFFF','#FFFFFF',green2blue(highnum))
+			}else if(lowfrac/highfrac < 1/2){ # most of the range is above 0
+				colours = c(red2yellow(lownum),'#FFFFFF','#FFFFFF',green2iceblue(highnum))
+			}else{
+				colours = c(darkred2yellow(lownum),'#FFFFFF','#FFFFFF',green2darkblue(highnum))
+			}
+		}
+	}
+	
+	# Now assess cases where just a positive or negative scale is required:
+	if((plottype=='positive') && (any(coolvars == variablename))){
+		colours = yellow2purpleCool(colourres)
+	}else if((plottype=='positive') && (any(warmvars == variablename))){
+		colours = yellow2purpleWarm(colourres)
+	}else if((plottype=='negative') && (any(coolvars == variablename))){
+		colours = purple2yellowWarm(colourres)
+	}else if((plottype=='negative') && (any(warmvars == variablename))){
+		colours = iceblue2green(colourres)
+	}
+	
+	return(colours)
+}
+
+#library(RColorBrewer) ## For some example colors
+
+# Function for desaturating colors by specified proportion
+DesaturateColours = function(colours, sat) {
+	library(colorspace) 
+    desat_hsv = diag(c(1, sat, 1)) %*% rgb2hsv(col2rgb(colours))
+    result = hsv(desat_hsv[1,], desat_hsv[2,], desat_hsv[3,])
+    return(result)
+}
diff --git a/pals/R/Constants.R b/pals/R/Constants.R
index 3cd71d9..b7ca43d 100644
--- a/pals/R/Constants.R
+++ b/pals/R/Constants.R
@@ -54,244 +54,3 @@ varIndex = function(varname,templateVersion){
 	return(idx)	
 }
 
-# Acceptable variable ranges
-GetVariableRanges = function(){
-	SWdown = c(0,1360) # surface incident shortwave rad [W/m^2]
-	LWdown = c(0,750)  # surface incident longwave rad [W/m^2]
-	Tair = c(200,333)  # near surface air temperature [K]
-	Qair = c(0,0.1)   # near surface specific humidity [kg/kg]
-	Rainf = c(0,0.05)  # rainfall rate [mm/s]
-	Snowf = c(0,0.03)  # snowfall rate [mm/s]
-	PSurf = c(50000,110000) # surface air pressure [Pa]
-	CO2air = c(160,2000) # near surface CO2 concentration [ppmv]
-	Wind = c(0,75)     # scalar windspeed [m/s]
-	Qle = c(-1000,1000) # latent heat flux [W/m^2]
-	Qh = c(-1000,1000) # sensible heat flux [W/m^2]
-	Qg = c(-1000,1000) # ground heat flux [W/m^2]
-	NEE = c(-100,100)    # met ecosystem exchange CO2 [umol/m^2/s]
-	GPP = c(-100,100)    # met ecosystem exchange CO2 [umol/m^2/s]
-	SWup = c(0,1350)   # reflected SW rad [W/m^2]
-	Rnet = c(-500,1250)# net absorbed radiation [W/m^2]
-	range = list(SWdown=SWdown,LWdown=LWdown,Tair=Tair,
-		Qair=Qair,Rainf=Rainf,Snowf=Snowf,PSurf=PSurf,
-		CO2air=CO2air,Wind=Wind,Qle=Qle,Qh=Qh,NEE=NEE,
-		Rnet=Rnet)
-	return(range)	
-}
-
-GetVariableDetails = function(request_names){
-	variable_list = list()
-	for(v in 1:length(request_names)){
-		var_details = list()
-		# Dimensional variables
-		if(request_names[v] == 'lat'){
-			var_details[['Name']] = c('y','lat','latitude','Lat','Latitude')
-			var_details[['UnitsName']] = c('degrees_north')
-			var_details[['Multiplier']] = c(1)
-			var_details[['Addition']] = c(0)
-			var_details[['UnitsText']] = ''
-			var_details[['PlotName']] = 'Latitude'
-		}else if(request_names[v] == 'lon'){
-			var_details[['Name']] = c('x','lon','longitude','Lon','Longitude')
-			var_details[['UnitsName']] = c('degrees_east')
-			var_details[['Multiplier']] = c(1)
-			var_details[['Addition']] = c(0)
-			var_details[['UnitsText']] = ''
-			var_details[['PlotName']] = 'Longitude'
-		# Met variables
-		}else if(request_names[v] == 'SWdown'){
-			var_details[['Name']] = c('SWdown')
-			var_details[['UnitsName']] = c('W/m2','W/m^2','Wm^{-2}','wm2','watt/m^2')
-			var_details[['Multiplier']] = c(1,1,1,1,1)
-			var_details[['Addition']] = c(0,0,0,0,0)
-			var_details[['UnitsText']] = 'W/'~m^{2}
-			var_details[['PlotName']] = 'Downward shortwave radiation'
-		}else if(request_names[v] == 'LWdown'){
-			var_details[['Name']] = c('LWdown')
-			var_details[['UnitsName']] = c('W/m2','W/m^2','Wm^{-2}','wm2','watt/m^2')
-			var_details[['Multiplier']] = c(1,1,1,1,1)
-			var_details[['Addition']] = c(0,0,0,0,0)
-			var_details[['UnitsText']] = 'W/'~m^{2}
-			var_details[['PlotName']] = 'Downward longwave radiation'
-		}else if(request_names[v] == 'Tair'){
-			var_details[['Name']] = c('Tair')
-			var_details[['UnitsName']] = c('K')
-			var_details[['Multiplier']] = c(1)
-			var_details[['Addition']] = c(0)
-			var_details[['UnitsText']] = 'K'
-			var_details[['PlotName']] = 'Surface air temperature'
-		}else if(request_names[v] == 'Qair'){
-			var_details[['Name']] = c('Qair')
-			var_details[['UnitsName']] = c('kg/kg','g/g')
-			var_details[['Multiplier']] = c(1,1)
-			var_details[['Addition']] = c(0,0)
-			var_details[['UnitsText']] = 'kg/kg'
-			var_details[['PlotName']] = 'Specific humidity'
-		}else if(request_names[v] == 'PSurf'){
-			var_details[['Name']] = c('PSurf')
-			var_details[['UnitsName']] = c('Pa','pa','hPa')
-			var_details[['Multiplier']] = c(1,1,100)
-			var_details[['Addition']] = c(0,0,0)
-			var_details[['UnitsText']] = 'Pa'
-			var_details[['PlotName']] = 'Surface air pressure'
-		}else if(request_names[v] == 'Rainf'){
-			var_details[['Name']] = c('Rainf')
-			var_details[['UnitsName']] = c('mm/s','kg/m^2/s', 'kg/m^2s', 'mm/day')
-			var_details[['Multiplier']] = c(1,1,1,86400)
-			var_details[['Addition']] = c(0,0,0,0)
-			var_details[['UnitsText']] = 'mm/s'
-			var_details[['PlotName']] = 'Specific humidity'
-		}else if(request_names[v] == 'Wind'){
-			var_details[['Name']] = c('Wind')
-			var_details[['UnitsName']] = c('m/s')
-			var_details[['Multiplier']] = c(1)
-			var_details[['Addition']] = c(0)
-			var_details[['UnitsText']] = 'm/s'
-			var_details[['PlotName']] = 'Specific humidity'
-		# Flux variables
-		}else if(request_names[v] == 'Qh'){
-			var_details[['Name']] = c('Qh','FSH') # FSH->CLM
-			var_details[['UnitsName']] = c('W/m2','W/m^2','Wm^{-2}','wm2','watt/m^2')
-			var_details[['Multiplier']] = c(1,1,1,1,1)
-			var_details[['Addition']] = c(0,0,0,0,0)
-			var_details[['UnitsText']] = 'W/'~m^{2}
-			var_details[['PlotName']] = 'Sensible heat flux'
-		}else if(request_names[v] == 'Qle'){
-			var_details[['Name']] = c('Qle','FCEV') # NB: FCEV is only PART of Qle for CLM
-			var_details[['UnitsName']] = c('W/m2','W/m^2','Wm^{-2}','wm2','watt/m^2')
-			var_details[['Multiplier']] = c(1,1,1,1,1)
-			var_details[['Addition']] = c(0,0,0,0,0)
-			var_details[['UnitsText']] = 'W/'~m^{2}
-			var_details[['PlotName']] = 'Latent heat flux'
-		}else if(request_names[v] == 'Rnet'){
-			var_details[['Name']] = c('Rnet')
-			var_details[['UnitsName']] = c('W/m2','W/m^2','Wm^{-2}','wm2','watt/m^2')
-			var_details[['Multiplier']] = c(1,1,1,1,1)
-			var_details[['Addition']] = c(0,0,0,0,0)
-			var_details[['UnitsText']] = 'W/'~m^{2}
-			var_details[['PlotName']] = 'Net radiation'
-		}else if(request_names[v] == 'SWnet'){
-			var_details[['Name']] = c('SWnet')
-			var_details[['UnitsName']] = c('W/m2','W/m^2','Wm^{-2}','wm2','watt/m^2')
-			var_details[['Multiplier']] = c(1,1,1,1,1)
-			var_details[['Addition']] = c(0,0,0,0,0)
-			var_details[['UnitsText']] = 'W/'~m^{2}
-			var_details[['PlotName']] = 'Net shortwave radiation'
-		}else if(request_names[v] == 'LWnet'){
-			var_details[['Name']] = c('LWnet')
-			var_details[['UnitsName']] = c('W/m2','W/m^2','Wm^{-2}','wm2','watt/m^2')
-			var_details[['Multiplier']] = c(1,1,1,1,1)
-			var_details[['Addition']] = c(0,0,0,0,0)
-			var_details[['UnitsText']] = 'W/'~m^{2}
-			var_details[['PlotName']] = 'Net longwave radiation'
-		}else if(request_names[v] == 'Qg'){
-			var_details[['Name']] = c('Qg')
-			var_details[['UnitsName']] = c('W/m2','W/m^2','Wm^{-2}','wm2','watt/m^2')
-			var_details[['Multiplier']] = c(1,1,1,1,1)
-			var_details[['Addition']] = c(0,0,0,0,0)
-			var_details[['UnitsText']] = 'W/'~m^{2}
-			var_details[['PlotName']] = 'Ground heat flux'
-		}else if(request_names[v] == 'NEE'){
-			var_details[['Name']] = c('NEE','FCO2') # FCO2->CLM
-			var_details[['UnitsName']] = c('umol/m2/s','mumol/m2/s','umol/m^2/s',
-				'umol/m2s','gC/m^2/s','gC/m2/s')
-			var_details[['Multiplier']] = c(1,1,1,1,1/(1.201E-5),1/(1.201E-5))
-			var_details[['Addition']] = c(0,0,0,0,0)
-			var_details[['UnitsText']] = ~mu~'mol/'~m^{2}~'/s' # default PALS units for plots
-			var_details[['PlotName']] = 'Net ecosystem exchange'
-		}else if(request_names[v] == 'GPP'){
-			var_details[['Name']] = c('GPP')
-			var_details[['UnitsName']] = c('umol/m2/s','mumol/m2/s','umol/m^2/s',
-				'umol/m2s','gC/m^2/s','gC/m2/s')
-			var_details[['Multiplier']] = c(1,1,1,1,1/(1.201E-5),1/(1.201E-5))
-			var_details[['Addition']] = c(0,0,0,0,0)
-			var_details[['UnitsText']] = ~mu~'mol/'~m^{2}~'/s' # default PALS units for plots
-			var_details[['PlotName']] = 'Gross primary production'
-		}
-		variable_list[[v]] = var_details
-	}
-	return(variable_list)
-}
-GetVariableIndex = function(vars,varname){
-	# Returns the index of a requested variable in a
-	# vars variable from a call to GetVariableDetails
-	vindex=NA
-	for(v in 1:length(vars)){
-		if(vars[[v]][['Name']][1] == varname){
-			vindex = v
-			break
-		}
-	}
-	return(vindex)	
-}
-
-# This is a list of alternative names, units and 
-# units conversions for variables of interest.
-SWdownNames = c('SWdown')
-QhNames = c('Qh','FSH') # FSH->CLM
-QleNames = c('Qle','FCEV') # NB: FCEV is only PART of Qle for CLM
-NEENames = c('NEE','FCO2') # FCO2->CLM
-GPPNames = c('GPP')
-RnetNames = c('Rnet')
-SWnetNames = c('SWnet')
-LWnetNames = c('LWnet')
-SoilMoistNames = c('SoilMoist')
-QgNames = c('Qg')
-# First listed units are default; conversions are relative to
-# order of listed units to convert to first units type.
-# Downward shortwave radiation:
-SWdownUnitsName = c('W/m2','W/m^2','Wm^{-2}','wm2','watt/m^2')
-SWdownMultiplier = c(1,1,1,1,1)
-SWdownAddition = c(0,0,0,0,0)
-SWdownUnits = list(name=SWdownUnitsName,multiplier=SWdownMultiplier,
-	addition=SWdownAddition)
-# Sensible heat:
-QhUnitsName = c('W/m2','W/m^2','Wm^{-2}','wm2','watt/m^2')
-QhMultiplier = c(1,1,1,1,1)
-QhAddition = c(0,0,0,0,0)
-QhUnits = list(name=QhUnitsName,multiplier=QhMultiplier,
-	addition=QhAddition)
-# Latent heat:
-QleUnitsName = c('W/m2','W/m^2','Wm^{-2}','wm2','watt/m^2')
-QleMultiplier = c(1,1,1,1,1)
-QleAddition = c(0,0,0,0,0)
-QleUnits = list(name=QleUnitsName,multiplier=QleMultiplier,
-	addition=QleAddition)
-# Net Ecosystem Exchange of CO2:
-NEEUnitsName = c('umol/m2/s','mumol/m2/s','umol/m^2/s',
-	'umol/m2s','gC/m^2/s','gC/m2/s')
-NEEMultiplier = c(1,1,1,1,1/(1.201E-5),1/(1.201E-5))
-NEEAddition = c(0,0,0,0,0,0)
-NEEUnits = list(name=NEEUnitsName,multiplier=NEEMultiplier,
-	addition=NEEAddition)
-# Gross Primary Productivity of CO2:
-GPPUnitsName = c('umol/m2/s','mumol/m2/s','umol/m^2/s',
-	'umol/m2s','gC/m^2/s','gC/m2/s')
-GPPMultiplier = c(1,1,1,1,1/(1.201E-5),1/(1.201E-5))
-GPPAddition = c(0,0,0,0,0,0)
-GPPUnits = list(name=GPPUnitsName,multiplier=GPPMultiplier,
-	addition=GPPAddition)
-# Net absorbed radiation:
-RnetUnitsName = c('W/m2','W/m^2','Wm^{-2}','wm2','watt/m^2')
-RnetMultiplier = c(1,1,1,1,1)
-RnetAddition = c(0,0,0,0,0)
-RnetUnits = list(name=RnetUnitsName,multiplier=RnetMultiplier,
-	addition=RnetAddition)
-# Net absorbed shortwave radiation:
-SWnetUnitsName = c('W/m2','W/m^2','Wm^{-2}','wm2','watt/m^2')
-SWnetMultiplier = c(1,1,1,1,1)
-SWnetAddition = c(0,0,0,0,0)
-SWnetUnits = list(name=SWnetUnitsName,multiplier=SWnetMultiplier,
-	addition=SWnetAddition)
-# Net absorbed longwave radiation:
-LWnetUnitsName = c('W/m2','W/m^2','Wm^{-2}','wm2','watt/m^2')
-LWnetMultiplier = c(1,1,1,1,1)
-LWnetAddition = c(0,0,0,0,0)
-LWnetUnits = list(name=LWnetUnitsName,multiplier=LWnetMultiplier,
-	addition=LWnetAddition)
-# Ground heat flux:
-QgUnitsName = c('W/m2','W/m^2','Wm^{-2}','wm2','watt/m^2')
-QgMultiplier = c(1,1,1,1,1)
-QgAddition = c(0,0,0,0,0)
-QgUnits = list(name=QgUnitsName,multiplier=QgMultiplier,
-	addition=QgAddition)
diff --git a/pals/R/DistributeAnalyses.R b/pals/R/DistributeAnalyses.R
index f53f5eb..0dcb783 100644
--- a/pals/R/DistributeAnalyses.R
+++ b/pals/R/DistributeAnalyses.R
@@ -2,31 +2,22 @@
 #
 # Functions to distribute analyses across multiple cores
 #
-# Gab Abramowitz, UNSW, 2014 (palshelp at gmail dot com)
+# Gab Abramowitz, UNSW, 2015 (palshelp at gmail dot com)
 
-DistributeGriddedAnalyses = function(Analysis,vars,obs,model,bench){
-	# Each call to this function will generate a single plot and its statistics
+DistributeGriddedAnalyses = function(Analysis,vars,obs,model,bench,region,cl){
+	# Each call to this function will generate a single plot and a list of metrics associate with it
 	
 	# Create outfilename:
 	outfile = setOutput('default')
 	
-	# For now, assumes analysis will be for a single variable.
-	# The name of this variable
-	varname = vars[[Analysis$vindex]][['Name']][1]
-	# Units expression:
-	unitstxt = vars[[Analysis$vindex]][['UnitsText']]
-	# Longer variable name for plots:
-	longvarname = vars[[Analysis$vindex]][['PlotName']]
-	# File name for graphics file:
-	filestring = paste(getwd(),outfile,sep = "/")
 	# Analysis identifier for javascript:
-	outfiletype = paste(varname,tolower(Analysis$type))
+	outfiletype = paste(vars[[Analysis$vindex]][['Name']][1],tolower(Analysis$type))
 	
-	# Check obs or model aren't missing variable data and that their timing is compatible:
+	# Check obs or model aren't missing variable, their timing is compatible and grids match:
 	errcheck = CanAnalysisProceed(obs, model)
 	if(errcheck$err){
-		result = list(type=outfiletype,filename=filestring,mimetype="image/png",
-			error=errcheck$errtext,bencherror=bench$errtext,metrics=list(first=list(name='fail',model_value=NA)))
+		result = list(type=outfiletype,filename=paste(getwd(),outfile,sep = "/"),mimetype="image/png",
+			error=errcheck$errtext,bencherror=bench$errtext,metrics=list(first=list(name='failed',model_value=NA)))
 		return(result)
 	}
 	
@@ -35,37 +26,43 @@ DistributeGriddedAnalyses = function(Analysis,vars,obs,model,bench){
 	
 	# Call analysis function:	
 	if(Analysis$type == 'TimeMean'){
-		bencherrtext = bench$errtext
-		areturn = SpatialAus(model,obs,bench,varname,unitstxt,longvarname,metrics,plottype=Analysis$type)				
+		metrics_data = TimeMeanAll(model,obs,bench,variable=vars[[Analysis$vindex]],plottype=Analysis$type,cl)
+		areturn = SpatialPlotAbsolute(model,obs,bench,metrics_data,
+			variable=vars[[Analysis$vindex]],plottype=Analysis$type,region)
 	}else if(Analysis$type == 'TimeSD'){
-		bencherrtext = bench$errtext
-		areturn = SpatialAus(model,obs,bench,varname,unitstxt,longvarname,metrics,plottype=Analysis$type)
+		metrics_data = TimeSDAll(model,obs,bench,variable=vars[[Analysis$vindex]],plottype=Analysis$type,cl)
+		areturn = SpatialPlotAbsolute(model,obs,bench,metrics_data,
+			variable=vars[[Analysis$vindex]],plottype=Analysis$type,region)
 	}else if(Analysis$type == 'TimeRMSE'){
-		bencherrtext = bench$errtext
-		areturn = SpatialAusRelative(model,obs,bench,varname,unitstxt,longvarname,metrics,plottype=Analysis$type)		
+		metrics_data = TimeRMSEAll(model,obs,bench,variable=vars[[Analysis$vindex]],plottype=Analysis$type,cl)
+		areturn = SpatialPlotRelative(model,obs,bench,metrics_data,
+			variable=vars[[Analysis$vindex]],plottype=Analysis$type,region)		
 	}else if(Analysis$type == 'TimeCor'){
-		bencherrtext = bench$errtext
-		areturn = SpatialAusRelative(model,obs,bench,varname,unitstxt,longvarname,metrics,plottype=Analysis$type)	
+		metrics_data = TimeCorAll(model,obs,bench,variable=vars[[Analysis$vindex]],plottype=Analysis$type,cl)
+		areturn = SpatialPlotRelative(model,obs,bench,metrics_data,
+			variable=vars[[Analysis$vindex]],plottype=Analysis$type,region)	
+	}else{
+		result = list(errtext = paste('Unknown analysis type \'',Analysis$type,
+			'\' requested in function DistributeGriddedAnalyses.',sep=''),err=TRUE)
+		return(result)
 	}
 	
 	if(areturn$errtext=='ok'){	
 		result = list(type=outfiletype,filename=paste(getwd(),outfile,sep = "/"),mimetype="image/png",
-			metrics = areturn$metrics,analysistype=Analysis$type, variablename=varname,bencherror=bencherrtext)
+			metrics = areturn$metrics,analysistype=Analysis$type, 
+			variablename=vars[[Analysis$vindex]][['Name']][1],bencherror=bench$errtext)
 	}else{
 		cat('\n###',areturn$errtext,'###\n')
 		result = list(type=outfiletype,filename=paste(getwd(),outfile,sep = "/"),mimetype="image/png",
-			metrics = areturn$metrics,analysistype=Analysis$type, variablename=varname,
-			error=areturn$errtext,bencherror=bencherrtext)
+			metrics = areturn$metrics,analysistype=Analysis$type, 
+			variablename=vars[[Analysis$vindex]][['Name']][1],error=areturn$errtext,bencherror=bench$errtext)
 	}	
 	
 	return(result)
 }
 
 DistributeSingleSiteAnalyses = function(Analysis,data,vars){
-	
-	# These will be metrics passed unless overwritten below
-	metrics = list(nme = 0, rmse=0,correlation=1)
-	
+		
 	# Create outfilename:
 	outfile = setOutput('ModelAnalysis')
 	
@@ -97,14 +94,14 @@ DistributeSingleSiteAnalyses = function(Analysis,data,vars){
 		# Check obs or model aren't missing variable data and and that their timing is compatible:
 		errcheck = CanAnalysisProceed(data[[Analysis$vindex]]$obs,data[[Analysis$vindex]]$model)
 		if(errcheck$err){
-			result = list(type=outfiletype,filename=filestring,mimetype="image/png",
+			result = list(type=outfiletype,filename=filestring,mimetype="image/png",analysistype=Analysis$type,
 				error=errcheck$errtext,bencherror=data[[Analysis$vindex]]$bench$errtext,
-				metrics=list(first=list(name='fail',model_value=NA)))
+				metrics=list(first=list(name='failed',model_value=NA)),variablename=varname)
 			return(result)
 		}
 		
 		# Test benchmark timing compatibility, and remove any benchmarks if necessary:
-		bench = PruneBenchmarks(data[[Analysis$vindex]]$obs,data[[Analysis$vindex]]$bench)
+		data[[Analysis$vindex]]$bench = PruneBenchmarks(data[[Analysis$vindex]]$obs,data[[Analysis$vindex]]$bench)
 		
 		# Create data matrix to send to analysis function:
 		adata=matrix(NA,length(data[[Analysis$vindex]]$obs$data),(2+data[[Analysis$vindex]]$bench$howmany))
@@ -134,7 +131,9 @@ DistributeSingleSiteAnalyses = function(Analysis,data,vars){
 				benchnames[b] = data[[Analysis$vindex]]$bench[[ data[[Analysis$vindex]]$bench$index[b] ]]$name
 			}
 		}
-		legendtext = c('Observed',moname,benchnames)
+		
+		legendtext = LegendText(data[[Analysis$vindex]],plotobs=TRUE)
+		plotcolours = BenchmarkColours(data[[Analysis$vindex]]$bench,plotobs=TRUE)
 		
 		# Call analysis function:	
 		if(Analysis$type == 'Timeseries'){
@@ -143,43 +142,32 @@ DistributeSingleSiteAnalyses = function(Analysis,data,vars){
 			winsize = 14
 			ytext=bquote('Smoothed'~.(tolower(longvarname)) ~ ' (' ~ .(unitstxt) ~ ')')
 			areturn = Timeseries(obsname,adata,varname,ytext,legendtext,plotcex,
-				data[[Analysis$vindex]]$obs$timing,smoothed=TRUE,winsize,moname,vqcdata=vqcdata)				
+				data[[Analysis$vindex]]$obs$timing,smoothed=TRUE,winsize,plotcolours,
+				moname,vqcdata=vqcdata)				
 		}else if(Analysis$type == 'AnnualCycle'){
 			bencherrtext = data[[Analysis$vindex]]$bench$errtext
 			ytext = bquote('Average'~.(tolower(longvarname)) ~ ' (' ~ .(unitstxt) ~ ')')
 			areturn = AnnualCycle(obsname,adata,varname,ytext,legendtext,
 				data[[Analysis$vindex]]$obs$timing$tstepsize,
-				data[[Analysis$vindex]]$obs$timing$whole,moname)
+				data[[Analysis$vindex]]$obs$timing$whole,plotcolours,moname)
 		}else if(Analysis$type == 'DiurnalCycle'){
 			bencherrtext = data[[Analysis$vindex]]$bench$errtext
 			ytext=bquote('Average'~.(varname) ~ ' (' ~.(unitstxt) ~ ')')
 			areturn = DiurnalCycle(obsname,adata,varname,ytext,legendtext,
 				data[[Analysis$vindex]]$obs$timing$tstepsize,
-				data[[Analysis$vindex]]$obs$timing$whole,moname,vqcdata=vqcdata)
+				data[[Analysis$vindex]]$obs$timing$whole,plotcolours,moname,vqcdata=vqcdata)
 		}else if(Analysis$type == 'PDF'){
 			bencherrtext = data[[Analysis$vindex]]$bench$errtext
 			nbins=500
 			xtext=bquote(.(longvarname) ~ ' (' ~ .(unitstxt) ~ ')')
 			areturn = PALSPdf(obsname,adata,varname,xtext,legendtext,
-				data[[Analysis$vindex]]$obs$timing,nbins,moname,vqcdata=vqcdata)
+				data[[Analysis$vindex]]$obs$timing,nbins,plotcolours,moname,vqcdata=vqcdata)
 		}else if(Analysis$type == 'Scatter'){
-			# Not a benhcmark plot for the moment:
-			bencherrtext = 'Benchmark analysis not available for this analysis type'
-			vtext = bquote(.(tolower(longvarname)) ~ ' (' ~.(unitstxt) ~ ')')
-			xytext = c('Observed','Modelled')
-			areturn = PALSScatter(obsname,data[[Analysis$vindex]]$model$data,
-				data[[Analysis$vindex]]$obs$data,varname,vtext,
-				xytext,data[[Analysis$vindex]]$obs$timing$tstepsize,
-				data[[Analysis$vindex]]$obs$timing$whole,ebal=FALSE,
-				modlabel=moname,vqcdata=vqcdata)
+			bencherrtext = data[[Analysis$vindex]]$bench$errtext
+			areturn = PALSScatter(data[[Analysis$vindex]],vars[[Analysis$vindex]],ebal=FALSE)
 		}else if(Analysis$type == 'Taylor'){
-			# Not a benhcmark plot for the moment:
-			bencherrtext = 'Benchmark analysis not available for this analysis type'
-			xtext=bquote(.(longvarname) ~ ' (' ~ .(unitstxt) ~ ')')
-			areturn = TaylorDiagram(obsname,data[[Analysis$vindex]]$model$data,
-				data[[Analysis$vindex]]$obs$data,varname,xtext,
-				data[[Analysis$vindex]]$obs$timing$tstepsize,
-				data[[Analysis$vindex]]$obs$timing$whole,moname)
+			bencherrtext = data[[Analysis$vindex]]$bench$errtext
+			areturn = TaylorDiagram(data[[Analysis$vindex]],vars[[Analysis$vindex]],plotcolours)
 		}else if(Analysis$type == 'AvWindow'){
 			# Not a benhcmark plot for the moment:
 			bencherrtext = 'Benchmark analysis not available for this analysis type'
@@ -190,90 +178,18 @@ DistributeSingleSiteAnalyses = function(Analysis,data,vars){
 		}
 		
 	}
-
+	# Don't return errtext in output list unless there is an error - as requested by Eden
 	if(areturn$errtext=='ok'){	
 		result = list(type=outfiletype,filename=paste(getwd(),outfile,sep = "/"),mimetype="image/png",
-			metrics = areturn$metrics,analysistype=Analysis$type, variablename=varname,bencherror=bencherrtext)
+			metrics = areturn$metrics,analysistype=Analysis$type, variablename=varname,
+			bencherror=bencherrtext,obsname=obsname,moname=moname,benchnames=benchnames)
 	}else{
 		cat('\n###',areturn$errtext,'###\n')
 		result = list(type=outfiletype,filename=paste(getwd(),outfile,sep = "/"),mimetype="image/png",
 			metrics = areturn$metrics,analysistype=Analysis$type, variablename=varname,
-			error=areturn$errtext,bencherror=bencherrtext)
+			error=areturn$errtext,bencherror=bencherrtext,obsname=obsname,moname=moname,benchnames=benchnames)
 	}	
 
 	return(result)
 }
 
-CanAnalysisProceed = function(obs,model){
-	# Checks obs, model variables were found and timing is appropriate.
-	
-	# Check for obs or model aren't missing variable data:
-	readcheck = CheckDataRead(obs$err,obs$errtext,model$err,model$errtext)
-	# Don't proceed and report error if there's an issue:		
-	if(readcheck$err){
-		return(readcheck)
-	}	
-	# Check model, obs timing consistency
-	tcheck = CheckTiming(obs$timing,model$timing)	
-	# Don't proceed and report error if there's an issue:	
-	if(tcheck$err){
-		return(tcheck)
-	}
-	proceed = list(err=FALSE)
-	return(proceed)
-}
-
-CheckDataRead = function(obserr,obserrtext,moderr,moderrtext){
-	# Simply reports whether there was a read error from either
-	# obs or model output reading (no file, appropriate variable etc)
-	errtext = 'ok'
-	err = FALSE
-	if(obserr){
-		err = TRUE
-		errtext = obserrtext
-	}else if(moderr){
-		err = TRUE
-		errtext = moderrtext
-	}	
-	result = list(err = err,errtext = errtext)
-	return(result)
-}
-
-PruneBenchmarks = function(obs,bench){
-	# Test benchmark timing compatibility, and remove benchmark(s) if necessary:
-	if(bench$exist){
-		# We'll need to use bench$index inside the for loop and also want to 
-		# modify it for future use, so modify new_benchindex instead, then overwrite bench$index:
-		new_benchindex = bench$index
-		for(b in 1: bench$howmany){
-			# Check benchmark and obs timing are compatible:
-			tcheck = CheckTiming(obs$timing,bench[[ bench$index[b] ]]$timing,benchmark_timing=TRUE)
-			if(tcheck$err){
-				# Report error with benchmark
-				bench$errtext = paste(bench$errtext,'Benchmark',bench$index[b],':',tcheck$errtext)
-				bench[[bench$index[b]]]$errtext = tcheck$errtext
-				# Remove benchmark from benchmark list:
-				bench$howmany = bench$howmany - 1
-				if(bench$howmany == 0){
-					# If that was the only benchmark, note there no longer any:
-					bench$exist = FALSE
-				}else{
-					# Change index of appropriate benchmarks:
-					oldlength = length(new_benchindex)
-					if(b==1){
-						new_benchindex = new_benchindex[2:oldlength]
-					}else if(b==oldlength){	
-						new_benchindex = new_benchindex[1:(oldlength-1)]
-					}else{
-						new_benchindex = 
-							c(new_benchindex[1:(b-1)],new_benchindex[(b+1):oldlength])
-					}
-				}
-			}
-		}
-		# Overwrite bench$index with values that account for any benchmarks that failed:
-		bench$index = new_benchindex
-	}
-	cat('Remaining benchmarks:',bench$howmany,' \n')
-	return(bench)	
-}
\ No newline at end of file
diff --git a/pals/R/GetFluxnet.R b/pals/R/GetFluxnet.R
index 406de0a..cd6a47d 100644
--- a/pals/R/GetFluxnet.R
+++ b/pals/R/GetFluxnet.R
@@ -48,8 +48,9 @@ GetFluxnetVariable = function(variable,filedetails,flagonly=FALSE){
 	if(! flagonly){ # if this function call is actually about fetching data:
 		timing = GetTimingNcfile(fid)
 		data=ncvar_get(fid,variable[['Name']][1])   # read observed variable data
+		grid = GetGrid(fid)
 		obs=list(data=data,timing=timing,qc=qc,qcexists=vexists$qc,name=filedetails$name,
-			err=FALSE,errtext=errtext)
+			grid=grid,err=FALSE,errtext=errtext)
 	}else{
 		obs=list(qcexists=vexists$qc,qc=qc,name=filedetails$name,
 			err=FALSE,errtext=errtext)
diff --git a/pals/R/GetGLEAM.R b/pals/R/GetGLEAM.R
new file mode 100644
index 0000000..0ca88b0
--- /dev/null
+++ b/pals/R/GetGLEAM.R
@@ -0,0 +1,137 @@
+# GetGLEAM.R
+#
+# This script fetches GLEAM ET data
+#
+# Gab Abramowitz, UNSW, 2015, gabsun at gmail dot com
+#
+
+GetGLEAM_Aus = function(variable,filelist,force_interval='no',dsetversion='default'){
+	library(ncdf4) # load package
+	# Should really access complete GLEAM data files and pull out Aus,
+	# but for now...
+	errtext='ok'	
+	if((variable[['Name']][1] != 'Qle') && (variable[['Name']][1] != 'Evap')){
+		errtext = 'Request for non-Qle, non-Evap variable to GetGLEAM_Aus read routine.'
+		obs = list(err=TRUE,errtext=errtext)
+		return(obs)
+	}
+	nyears = length(filelist)
+	year = c() 
+	for(f in 1:nyears){ # For each file sent by js
+		# Establish which year the file contains:
+		year[f] = as.numeric(substr(filelist[[f]][['path']], 
+			(nchar(filelist[[f]][['path']])-6), (nchar(filelist[[f]][['path']])-3) ) )
+	}
+	# Define the order to read files:
+	fileorder = order(year)	
+	# Define number of days in total:
+	yds = Yeardays(min(year),nyears*366)
+	daysvector = yds$daysperyear
+	ndays = sum(daysvector[1:nyears])
+	dayctr = 1 # initialise
+	if((force_interval == 'no') | (force_interval == 'daily')){
+		interval = 'daily'
+		tsteps = ndays
+		ET = array(NA,dim=c(208,143,ndays))	# Initialise data array:
+	}else if(force_interval == 'monthly'){
+		interval = 'monthly'
+		tsteps = nyears*12
+		ET = array(NA,dim=c(208,143,tsteps))	# Initialise data array:
+	}else{
+		errtext = paste('GLEAM_Aus requested to force to unknown interval:',force_interval)
+		obs = list(err=TRUE,errtext=errtext)
+		return(obs)	
+	}
+	# Get data:
+	for(f in 1:nyears){ # For each file sent by js		
+		# Open file:
+		fid = nc_open(filelist[[ fileorder[f] ]][['path']],write=FALSE,readunlim=FALSE)
+		# Read GLEAM data for this year:
+		if(interval == 'daily'){
+			ET[,,dayctr:(dayctr + daysvector[f] - 1)] = ncvar_get(fid, 'EVAP' ) # read model output data
+		}else if(interval == 'monthly'){
+			tmp = ncvar_get(fid, 'EVAP' ) # read model output data
+			ET[,, ((f-1)*12+1) : ((f-1)*12+12)] = DailyToMonthly(tmp,year[fileorder[f]],daysvector[f])
+		}	
+		# Increment counter:
+		dayctr = dayctr + daysvector[f]
+		# Close netcdf file for this year:
+		nc_close(fid)
+	}
+	# Reopen first file to fetch lat and lon:
+	fid = nc_open(filelist[[ 1 ]][['path']],write=FALSE,readunlim=FALSE)
+	# Then get spatial grid structure from first model output file:
+	grid = GetGrid(fid)
+	if(grid$err){	
+		obs = list(err=TRUE,errtext=grid$errtext)
+		nc_close(fid) # Close netcdf file
+		return(obs)
+	}
+	nc_close(fid) # Close netcdf file
+	
+	if(variable[['Name']][1] == 'Qle'){
+		ET = ET*28.4 # convert from mm/day to W/m^2
+	}
+	
+	timing = list(interval=interval,tsteps=tsteps)
+	
+	# Return result
+	obs = list(err=FALSE,errtext=errtext,data=ET,grid=grid,timing=timing,name='GLEAM ET')
+	return(obs)	
+}
+
+GetGLEAM_Global = function(variable,filelist,force_interval='no',dsetversion='default'){
+	library(ncdf4) # load package
+	errtext='ok'	
+	if((variable[['Name']][1] != 'Qle') && (variable[['Name']][1] != 'Evap')){
+		errtext = 'Request for non-Qle, non-Evap variable to GetGLEAM_Global read routine.'
+		obs = list(err=TRUE,errtext=errtext)
+		return(obs)
+	}
+	nyears = length(filelist)
+	year = c() 
+	for(f in 1:nyears){ # For each file sent by js
+		# Establish which year the file contains:
+		year[f] = as.numeric(substr(filelist[[f]][['path']], 
+			(nchar(filelist[[f]][['path']])-6), (nchar(filelist[[f]][['path']])-3) ) )
+	}
+	# Define the order to read files:
+	fileorder = order(year)	
+	# Define number of days in total:
+	if((force_interval == 'no') | (force_interval == 'monthly')){
+		interval = 'monthly'
+		tsteps = nyears*12
+		ET = array(NA,dim=c(720,360,tsteps))	# Initialise data array:
+	}else{
+		errtext = paste('GetGLEAM_Global requested to force to unknown interval:',force_interval)
+		obs = list(err=TRUE,errtext=errtext)
+		return(obs)	
+	}
+	# Get data:
+	for(f in 1:nyears){ # For each file sent by js		
+		# Open file:
+		fid = nc_open(filelist[[ fileorder[f] ]][['path']],write=FALSE,readunlim=FALSE)
+		# Read GLEAM data for this year:
+		if(interval == 'monthly'){
+			ET[,, ((f-1)*12+1) : ((f-1)*12+12)] = ncvar_get(fid, 'et' ) # read model output data
+		}	
+		# Close netcdf file for this year:
+		nc_close(fid)
+	}
+	# Reopen first file to fetch lat and lon:
+	fid = nc_open(filelist[[ 1 ]][['path']],write=FALSE,readunlim=FALSE)
+	# Then get spatial grid structure from first model output file:
+	grid = GetGrid(fid)
+	if(grid$err){	
+		obs = list(err=TRUE,errtext=grid$errtext)
+		nc_close(fid) # Close netcdf file
+		return(obs)
+	}
+	nc_close(fid)
+	
+	timing = list(interval=interval,tsteps=tsteps)
+	
+	# Return result
+	obs = list(err=FALSE,errtext=errtext,data=ET,grid=grid,timing=timing,name='GLEAM ET')
+	return(obs)	
+}
\ No newline at end of file
diff --git a/pals/R/GetGLEAM_Aus.R b/pals/R/GetGLEAM_Aus.R
deleted file mode 100644
index a61fb77..0000000
--- a/pals/R/GetGLEAM_Aus.R
+++ /dev/null
@@ -1,83 +0,0 @@
-# GetGLEAM_Aus.R
-#
-# This script fetches GLEAM ET data masked for Australia
-#
-# Gab Abramowitz, UNSW, 2013, gabsun at gmail dot com
-#
-# Inputs:
-# missing_threshold = what percentage of missing data in a grid
-# 					cell is acceptable?
-#
-# Default values are those specified in function arguments:
-
-GetGLEAM_Aus = function(variable,filelist,force_interval='no',dsetversion='default'){
-	library(ncdf4) # load package
-	# Should really access complete GLEAM data files and pull out Aus,
-	# but for now...
-	errtext='ok'	
-	if((variable[['Name']][1] != 'Qle') && (variable[['Name']][1] != 'Evap')){
-		errtext = 'Request for non-Qle, non-Evap variable to GLEAM_Aus read routine.'
-		obs = list(err=TRUE,errtext=errtext)
-		return(obs)
-	}
-	nyears = length(filelist)
-	year = c() 
-	for(f in 1:nyears){ # For each file sent by js
-		# Establish which year the file contains:
-		year[f] = as.numeric(substr(filelist[[f]][['path']], 
-			(nchar(filelist[[f]][['path']])-6), (nchar(filelist[[f]][['path']])-3) ) )
-	}
-	# Define the order to read files:
-	fileorder = order(year)	
-	# Define number of days in total:
-	yds = Yeardays(min(year),nyears*366)
-	daysvector = yds$daysperyear
-	ndays = sum(daysvector[1:nyears])
-	dayctr = 1 # initialise
-	if((force_interval == 'no') | (force_interval == 'daily')){
-		interval = 'daily'
-		tsteps = ndays
-		ET = array(NA,dim=c(208,143,ndays))	# Initialise data array:
-	}else if(force_interval == 'monthly'){
-		interval = 'monthly'
-		tsteps = nyears*12
-		ET = array(NA,dim=c(208,143,tsteps))	# Initialise data array:
-	}else{
-		errtext = paste('GLEAM_Aus requested to force to unknown interval:',force_interval)
-		obs = list(err=TRUE,errtext=errtext)
-		return(obs)	
-	}
-	# Get data:
-	for(f in 1:nyears){ # For each file sent by js		
-		# Open file:
-		fid = nc_open(filelist[[ fileorder[f] ]][['path']],write=FALSE,readunlim=FALSE)
-		# Read GLEAM data for this year:
-		if(interval == 'daily'){
-			ET[,,dayctr:(dayctr + daysvector[f] - 1)] = ncvar_get(fid, 'EVAP' ) # read model output data
-		}else if(interval == 'monthly'){
-			tmp = ncvar_get(fid, 'EVAP' ) # read model output data
-			ET[,, ((f-1)*12+1) : ((f-1)*12+12)] = DailyToMonthly(tmp,year[fileorder[f]],daysvector[f])
-		}	
-		# Increment counter:
-		dayctr = dayctr + daysvector[f]
-		# Close netcdf file for this year:
-		nc_close(fid)
-	}
-	# Reopen first file to fetch lat and lon:
-	fid = nc_open(filelist[[ 1 ]][['path']],write=FALSE,readunlim=FALSE)
-	# Fetch lat and lon:
-	lat=ncvar_get(fid, 'lat' ) # read latitude 
-	lon=ncvar_get(fid, 'lon' ) # read latitude 
-	# Close netcdf file :
-	nc_close(fid)
-	
-	if(variable[['Name']][1] == 'Qle'){
-		ET = ET*28.4 # convert from mm/day to W/m^2
-	}
-	
-	timing = list(interval=interval,tsteps=tsteps)
-	
-	# Return result
-	obs = list(err=FALSE,errtext=errtext,data=ET,lat=lat,lon=lon,timing=timing,name='GLEAM')
-	return(obs)	
-}
\ No newline at end of file
diff --git a/pals/R/GetModelOutput.R b/pals/R/GetModelOutput.R
index 4b3a880..853d803 100644
--- a/pals/R/GetModelOutput.R
+++ b/pals/R/GetModelOutput.R
@@ -3,7 +3,7 @@
 # Reads a variable from model output netcdf file
 # Gab Abramowitz UNSW 2014 (palshelp at gmail dot com)
 #
-GetModelOutput = function(variable,filelist){
+GetModelOutput = function(variable,filelist,forcegrid='no'){
 	library(ncdf4) # load netcdf library
 	errtext='ok'
 	
@@ -15,7 +15,7 @@ GetModelOutput = function(variable,filelist){
 	for(f in 1:length(filelist)){ # For each file of this MO:
 		# Check file exists:
 		if(!file.exists(filelist[[f]][['path']])){
-			errtext = paste('M4: Model output file',filelist[[f]][['path']],'does not exist.')
+			errtext = paste('Model output file',filelist[[f]][['path']],'does not exist.')
 			model=list(errtext=errtext,err=TRUE)
 			return(model)	
 		}
@@ -24,7 +24,7 @@ GetModelOutput = function(variable,filelist){
 		# Check that requested variable exists:
 		exists = AnyNcvarExists(mfid[[f]],variable[['Name']])
 		if( ! exists$var){
-			errtext = paste('M3: Requested variable',variable[['Name']][1],
+			errtext = paste('Requested variable',variable[['Name']][1],
 				'does not appear to exist in Model Ouput:', filelist[[f]][['name']])
 			model=list(errtext=errtext,err=TRUE)
 			mfid = lapply(mfid, nc_close)
@@ -47,6 +47,13 @@ GetModelOutput = function(variable,filelist){
 			return(model)
 		}
 	}
+	# Then get spatial grid structure from first model output file:
+	grid = GetGrid(mfid[[1]])
+	if(grid$err){	
+		model = list(err=TRUE,errtext=grid$errtext)
+		mfid = lapply(mfid, nc_close)
+		return(model)
+	}
 	
 	# Now try to ascertain how the data in these files needs to be assembled:
 	allintervals = c()
@@ -57,17 +64,12 @@ GetModelOutput = function(variable,filelist){
 		alltsteps[f] = modeltiming[[f]]$tsteps
 		allyear[f] = modeltiming[[f]]$syear
 	}
+	########### Monthly data in one year files ####################
 	if((all(allintervals == 'monthly')) && (all(alltsteps == 12))){		
 		# i.e. all are monthly data in year length files
 		nyears = length(filelist)
 		ntsteps = 12 * nyears # total number of time steps
-		# Get lat / lon from file
-		latlon = GetLatLon(mfid[[1]])
-		if(latlon$err){	
-			model = list(err=TRUE,errtext=latlon$errtext)
-			mfid = lapply(mfid, nc_close)
-			return(model)
-		}
+		
 		# Check that MO files don't repeat years:
 		if(length(unique(allyear)) != length(allyear)){ # i.e. a repeated year has been removed
 			errtext='Model output has two files with the same starting year.'
@@ -77,7 +79,7 @@ GetModelOutput = function(variable,filelist){
 		}
 		
 		# Allocate space for data:
-		vdata = array(NA,dim=c(latlon$lonlen,latlon$latlen,ntsteps) )
+		vdata = array(NA,dim=c(grid$lonlen,grid$latlen,ntsteps) )
 		# Define the order to read files:
 		fileorder = order(allyear)
 		
@@ -101,19 +103,14 @@ GetModelOutput = function(variable,filelist){
 		modeltimingall = list(tstepsize=modeltiming[[1]]$tstepsize,tsteps=ntsteps,
 			syear=modeltiming[[1]]$syear,smonth=modeltiming[[1]]$smonth,sdoy=modeltiming[[1]]$sdoy,
 			interval=modeltiming[[1]]$interval)
+	########### Model time step data with same number of time steps in each file ####################
 	}else if((all(allintervals == 'timestep')) && (all(alltsteps == alltsteps[1]))){
 		# i.e. all are per time step data with the same number of time steps in each file
 		ntsteps = alltsteps[1] * length(filelist) # total number of time steps
 		tsteps1 = alltsteps[1]
-		# Get lat / lon from file
-		latlon = GetLatLon(mfid[[1]])
-		if(latlon$err){	
-			model = list(err=TRUE,errtext=latlon$errtext)
-			mfid = lapply(mfid, nc_close)
-			return(model)
-		}
+		
 		# Allocate space for data:
-		vdata = array(NA,dim=c(latlon$lonlen,latlon$latlen,ntsteps) )
+		vdata = array(NA,dim=c(grid$lonlen,grid$latlen,ntsteps) )
 		# Define the order to read files:
 		fileorder = order(allyear)
 		# Get data:
@@ -132,12 +129,12 @@ GetModelOutput = function(variable,filelist){
 			for(f in 1:length(filelist)){ # For each file sent by js
 				vdata_tmp = ncvar_get(mfid[[ fileorder[f] ]],variable[['Name']][exists$index],collapse_degen=FALSE)
 
-				if ((variable[['Name']][1]=='NEE') & (length(vdata_tmp) != (ntsteps*latlon$lonlen*latlon$latlen))) {
+				if ((variable[['Name']][1]=='NEE') & (length(vdata_tmp) != (ntsteps*grid$lonlen*grid$latlen))) {
 					# likely an ORCHIDEE file where NEE has dim (x,y,t,vegtype), in which case sum over
 					# vegtype dim - NEE values are already weighted by vegtype fraction:
 					vdata_tmp = apply(vdata_tmp,c(1,2,4),sum)
 				}
-				if(length(vdata_tmp) != (ntsteps*latlon$lonlen*latlon$latlen)){
+				if(length(vdata_tmp) != (ntsteps*grid$lonlen*grid$latlen)){
 					errtext = paste('Requested variable',variable[['Name']][1],
 							'has more dimensions than expected in Model Ouput:', filelist[[f]][['name']])
 					model = list(err=TRUE,errtext=errtext)
@@ -164,102 +161,39 @@ GetModelOutput = function(variable,filelist){
 	
 	# Apply any units changes:
 	vdata = vdata*units$multiplier + units$addition
+	
+	# Assess nature of grid here
+	
+	# If not equal to requested grid structure (if forcegrid != 'no'), reshape data here
+	
 	# Create list to return from function:	
-	model=list(data=vdata,timing = modeltimingall,name=filelist[[1]]$name,
+	model=list(data=vdata,timing = modeltimingall,name=filelist[[1]]$name,grid=grid,
 		err=FALSE,errtext=errtext)
 	return(model)
 }
 
-GetLatLon = function(mfid){
-	# Gets the latitude and longitide dimensions from a model output netcdf file
-	errtext='ok'
-	# Fetch accepted names for latitude:
-	lat_det = GetVariableDetails('lat')
-	# Check an acceptable latitude variable exists:
-	exists_lat = AnyNcvarExists(mfid,lat_det[[1]][['Name']])
-	if(! exists_lat$var){
-		errtext = paste('Could not find latitude variable in',stripFilename(mfid$filename))
-		latlon = list(err=TRUE,errtext=errtext)
-		return(latlon)
-	}
-	# Now fetch latitude information:
-	if(exists_lat$dimvar){
-		# either as a dimension variable, if that was the found acceptable name:
-		lat = mfid$dim[[exists_lat$dimnum]]$vals
-		latlen = mfid$dim[[exists_lat$dimnum]]$len
-	}else{
-		# or as a normal variable:
-		if(mfid$var[[exists_lat$varnum]]$ndims != 1){ # If not a 1D latitude variable
-			# Check that it's still a lat-lon grid:
-			lat = ncvar_get(mfid,lat_det[[1]][['Name']][exists_lat$index]) # get lat data
-			latlen = length(unique(as.vector(lat))) # i.e. the number of unique latitude values
-			if(any(mfid$var[exists_lat$index]$size != latlen)){
-				# i.e. number of unique lat vals is not equal to the length of either of the dims
-				# on which latitude depends
-				errtext = 'PALS cannot yet cope with non lat-lon grids yet.'	
-				latlon = list(err=TRUE,errtext=errtext)
-				return(latlon)
-			}
-		}
-	}
-	
-	# Repeat the process for longitude:
-	# Fetch accepted names for longitude:
-	lon_det = GetVariableDetails('lon')
-	# Check an acceptable longitude variable exists:
-	exists_lon = AnyNcvarExists(mfid,lon_det[[1]][['Name']])
-	if(! exists_lon$var){
-		errtext = paste('Could not find longitude variable in',stripFilename(mfid$filename))
-		latlon = list(err=TRUE,errtext=errtext)
-		return(latlon)
-	}
-	# Now fetch longitude information:
-	if(exists_lon$dimvar){
-		# either as a dimension variable, if that was the found acceptable name:
-		lon = mfid$dim[[exists_lon$dimnum]]$vals
-		lonlen = mfid$dim[[exists_lon$dimnum]]$len
-	}else{
-		# or as a normal variable:
-		if(mfid$var[[exists_lon$varnum]]$ndims != 1){ # If not a 1D latitude variable
-			# Check that it's still a lat-lon grid:
-			lon = ncvar_get(mfid,lon_det[[1]][['Name']][exists_lon$index]) # get lon data
-			lonlen = length(unique(as.vector(lon))) # i.e. the number of unique latitude values
-			if(any(mfid$var[exists_lat$index]$size != latlen)){
-				# i.e. number of unique lat vals is not equal to the length of either of the dims
-				# on which latitude depends
-				errtext = 'PALS cannot yet cope with non lat-lon grids yet.'	
-				latlon = list(err=TRUE,errtext=errtext)
-				return(latlon)
-			}
-		}
-	}
-	# Return result:
-	latlon = list(err=FALSE,errtext=errtext,lat=lat,lon=lon,latlen=latlen,lonlen=lonlen)
-	return(latlon)
-}
-
-GetBenchmarks = function(variable,filelist,nBench){
+GetBenchmarks = function(variable,BenchmarkFiles,BenchInfo){
 	# Collects the model outputs that form benchmark simulations as a list
 	errtext = 'ok'
-	bench = list()
-	if(nBench$number > 0){
+	bench = list() # initialise benchmark data llist
+	if(BenchInfo$number > 0){
 		# Save for nBench$number positions in list for benchmark data;
 		# they'll be filled shortly:
-		for(b in 1:nBench$number){
+		for(b in 1:BenchInfo$number){
 			bench[[b]] = 0
 		}
     	bench[['exist']] = TRUE
     	bench[['howmany']] = 0
     	bench[['index']] = c()
     	bench[['errtext']] = ' '
-    	for(b in 1:nBench$number){
+    	for(b in 1:BenchInfo$number){
     		bench[['howmany']] = bench[['howmany']] + 1
     		# Select those files in file list that correspond to benchmark 'b'
     		thisbenchfiles = list()
-    		for(f in 1:length(nBench$benchfiles[[b]])){
-    			thisbenchfiles[[f]] = filelist[[ nBench$benchfiles[[b]][f] ]]
+    		for(f in 1:length(BenchInfo$benchfiles[[b]])){
+    			thisbenchfiles[[f]] = BenchmarkFiles[[ BenchInfo$benchfiles[[b]][f] ]]
     		}
-	    	bench[[b]] = GetModelOutput(variable,thisbenchfiles)
+	    	bench[[b]] = GetModelOutput(variable,thisbenchfiles)	    	
 	    	if(bench[[b]]$err){ # i.e. there was an error of some sort retrieving benchmark
 	    		# Add to this benchamrk's errtext field - note it's a benchmark:
 	    		bench[[b]]$errtext = paste('Benchmark error: ',bench[[b]]$errtext,sep='')
@@ -283,91 +217,40 @@ GetBenchmarks = function(variable,filelist,nBench){
 	return(bench)
 }
 
-NcvarExists = function(fid,varname){
-	# Checks that variable exists in netcdf file, and additionally 
-	# checks whether quality control couterpart variable exists:
-	exists_var = FALSE
-	exists_qc = FALSE
-	varnum = NULL
-	dimnum = NULL
-	dimvar = FALSE
-	for (v in 1:fid$nvars){ # Search through all variables in netcdf file
-		if(fid$var[[v]]$name==varname){
-			exists_var=TRUE
-			varnum = v
-			if(exists_var & exists_qc) {break}
-		}else if(fid$var[[v]]$name==paste(varname,'_qc',sep='')){
-			exists_qc=TRUE
-			if(exists_var & exists_qc) {break}	
-		}
-	}
-	# If not found as a variable, it may be a dimension variable:
-	if(! exists_var){
-		for(d in 1:length(fid$dim)){
-			# If a dimension name & dimension variable
-			if(fid$dim[[d]]$name == varname){
-				if(fid$dim[[d]]$dimvarid$id != -1){# i.e. dim var exists
-					exists_var=TRUE
-					dimvar = TRUE
-					dimnum = d
+PruneBenchmarks = function(obs,bench){
+	# Test benchmark timing compatibility, and remove benchmark(s) if necessary:
+	if(bench$exist){
+		# We'll need to use bench$index inside the for loop and also want to 
+		# modify it for future use, so modify new_benchindex instead, then overwrite bench$index:
+		new_benchindex = bench$index
+		for(b in 1: bench$howmany){
+			# Check benchmark and obs timing are compatible:
+			tcheck = CheckTiming(obs$timing,bench[[ bench$index[b] ]]$timing,benchmark_timing=TRUE)
+			if(tcheck$err){
+				# Report error with benchmark
+				bench$errtext = paste(bench$errtext,'Benchmark',bench$index[b],':',tcheck$errtext)
+				bench[[bench$index[b]]]$errtext = tcheck$errtext
+				# Remove benchmark from benchmark list:
+				bench$howmany = bench$howmany - 1
+				if(bench$howmany == 0){
+					# If that was the only benchmark, note there no longer any:
+					bench$exist = FALSE
+				}else{
+					# Change index of appropriate benchmarks:
+					oldlength = length(new_benchindex)
+					if(b==1){
+						new_benchindex = new_benchindex[2:oldlength]
+					}else if(b==oldlength){	
+						new_benchindex = new_benchindex[1:(oldlength-1)]
+					}else{
+						new_benchindex = 
+							c(new_benchindex[1:(b-1)],new_benchindex[(b+1):oldlength])
+					}
 				}
 			}
 		}
+		# Overwrite bench$index with values that account for any benchmarks that failed:
+		bench$index = new_benchindex
 	}
-	ncvexists = list(var=exists_var,qc=exists_qc,dimvar=dimvar,dimnum=dimnum,varnum=varnum)
-	return(ncvexists)	
-}
-
-AnyNcvarExists = function(fid,varnames){
-	# Checks whether a variable exists as any of its alternative standard names.
-	exists_var = FALSE
-	exists_qc = FALSE
-	vindex = 0
-	for(v in 1:length(varnames)){
-		exists = NcvarExists(fid,varnames[v])
-		if(exists$var){
-			exists_var=TRUE
-			exists_qc = exists$qc
-			vindex = v
-			break
-		}	
-	}
-	ncvexists = list(var=exists_var,qc=exists_qc,index=vindex,dimvar=exists$dimvar,
-		dimnum=exists$dimnum,varnum=exists$varnum)
-	return(ncvexists)	
-}
-
-CheckNcvarUnits = function(fid,vname,variable,file){
-	# Check whether units exist for a variable, and if so, whether
-	# they correspond to any standard version of units text:
-	# Get units for variable
-	errtext='ok'
-	mvunits=ncatt_get(fid,vname,'units')
-	if(! mvunits$hasatt){
-		errtext = paste('Variable',vname,'in file',
-			file,'does not have a units attribute.')
-		units = list(errtext = errtext, err=TRUE)
-		return(units)
-	}
-	UnitsExist = TRUE
-	UnitsMatch = FALSE
-	for (u in 1:length(variable$UnitsName)){
-		if(mvunits$value==variable$UnitsName[u]){ # i.e. found units match
-			UnitsMatch = TRUE
-			# Set units adjustments appropriately:
-			multiplier = variable$Multiplier[u]
-			addition = variable$Addition[u]
-			break # out of units name for loop
-		}
-	}
-	if( ! UnitsMatch){ # i.e. didn't recognise variable units
-		units=mvunits$value
-		errtext = paste('M3: Did not recognise units',mvunits$value,'for',
-			fid$var[[v]]$name,'in model output file',file)
-		units = list(errtext = errtext, err=TRUE)
-		return(units)	
-	}
-	units = list(exist = UnitsExist, match=UnitsMatch, value = mvunits$value,
-		multiplier=multiplier, addition=addition,err=FALSE,errtext=errtext)
-	return(units)
+	return(bench)	
 }
\ No newline at end of file
diff --git a/pals/R/MetricsTable.R b/pals/R/MetricsTable.R
new file mode 100644
index 0000000..8562ea5
--- /dev/null
+++ b/pals/R/MetricsTable.R
@@ -0,0 +1,310 @@
+# Collates metrics from analysis script into a PNG table
+# Gab Abramowitz, UNSW, 2014 (palshelp at gmail dot com)
+
+MetricTableSingleSite = function(outinfo,BenchInfo){
+	# Draws in a single image:
+	# summary table: just number of 'wins' for model and all benchmarks
+	# relationship analysis table: e.g. evapFrac; water use efficiency
+	# variable analysis table: metrics that are repeated for each analysis type
+	library(plotrix)
+	
+	# First check that all previous analyses didn't fail 
+	# (in which case don't bother with table):
+	if(CheckIfAllFailed(outinfo)){
+		return(outinfo)
+	}
+	
+	# Create analysis name for javascript:
+	thisanalysistype = 'Summary'
+	
+	# Number of significant figures in table data:
+	sfig = 2
+	# Table spacing:
+	xpadding = 0.2
+	ypadding = 1.0
+	# Table text sizing:
+	tablecex = 0.95
+	# Table grid:
+	xpos=c(-0.3,5,8)
+	ypos=c(5.3,0)
+	
+	# First get obs, model and benchmark names:
+	obsname = outinfo[[1]]$obsname
+	moname = outinfo[[1]]$moname
+	##### ************ get names from javascript read in instead!
+	benchnames = BenchInfo$names
+	
+	# Then collect model and benchmark metric value vectors for all analyses into a data frame.
+	# Initialisations:
+	ctr = 0
+	analysistype = c()
+	variable = c()
+	metricname = c()
+	model = c()
+	bench1 = c()
+	bench2 = c()
+	bench3 = c()
+	winner = c()
+	# Collate metric data:
+	for(i in 1: length(outinfo)){ # for each analysis performed
+		if(outinfo[[i]]$variablename != 'multiple'){ # i.e. this is a single variable analysis		
+			nmetrics = length(outinfo[[i]]$metrics) # Number of metrics for this analysis
+			for(m in 1:nmetrics){				
+				if(outinfo[[i]]$metrics[[m]]$name != 'failed'){ # if analysis was successful
+					# First get obs DS, MO name for plot title (need only happen once, overwritten):
+					obsname = outinfo[[i]]$obsname
+					moname = outinfo[[i]]$moname
+					
+					ctr = ctr + 1 # start index for this analysis
+					analysistype[ctr] = outinfo[[i]]$analysistype
+					variable[ctr] = outinfo[[i]]$variablename
+					metricname[ctr] = outinfo[[i]]$metrics[[m]]$name
+					model[ctr] = outinfo[[i]]$metrics[[m]]$model_value
+					# Initialise benchmark metric values for this analysis:
+					bench1[ctr] = NA
+					bench2[ctr] = NA
+					bench3[ctr] = NA
+					# Note that bench 1 in any given analysis in outinfo may not correspond to the user's
+					# first nominated benchmark - a relevant variable may be missing from the benchmark MO,
+					# or there may have been an error reading. The AlignBenchmarks function below deals with this:
+					if(BenchInfo$number >= 1){
+						bench1[ctr] = 
+							AlignBenchmarks(benchnames,outinfo[[i]]$benchnames,1,outinfo[[i]]$metrics[[m]]$bench_value)
+					}
+					if(BenchInfo$number >= 2){
+						bench2[ctr] = 
+							AlignBenchmarks(benchnames,outinfo[[i]]$benchnames,2,outinfo[[i]]$metrics[[m]]$bench_value)
+					}
+					if(BenchInfo$number >= 3){
+						bench3[ctr] = 
+							AlignBenchmarks(benchnames,outinfo[[i]]$benchnames,3,outinfo[[i]]$metrics[[m]]$bench_value)
+					}
+					# Note whether model or bench 1, 2 or 3 'won' this metric:
+					winner[ctr] = BestInMetric(metricname[ctr],model[ctr],bench1[ctr],bench2[ctr],bench3[ctr])
+				}
+			}
+		}	
+	}
+	# Create data frame with vars*analyses rows:
+	dfall = data.frame(analysis=analysistype,variable=variable,metric=metricname,
+		mvalue=model,b1value=bench1,b2value=bench2,b3value=bench3,winner=winner)
+		
+	# Create vector of unique variable names:
+	vars = unique(dfall[,'variable'])
+	halfvars = ceiling(length(vars)/2) # used to split table
+	split_table = FALSE
+	# Define row names:
+	initnames = paste(dfall[variable==vars[1],][,3],' (',dfall[variable==vars[1],][,1],')',sep='')
+	dfrownames = PadRowNames(initnames)
+	ctr = 1
+	vctr = c() # initialise index of columns for each variable in final data frame
+	
+	# First (column) entry in final table data frame is model metric values (for all analyses) for var 1:
+	dffinal = data.frame(a=signif(dfall[variable==vars[1],][,4],sfig))
+	dfcolumns = c(paste('Mod',vars[1]))
+	# Create colours data frame and add first vector:
+	dfcolours = data.frame(row1 = ifelse(c((dfall[variable==vars[1],][,8]) == 1),TableWinColour(),TableModelColour()))
+	if(!all(is.na(dfall[variable==vars[1],][,5]))){ # if any bench1 values for var 1 (for all analyses) exist
+		ctr = ctr + 1
+		dffinal[ctr] = signif(dfall[variable==vars[1],][,5],sfig) # benchmark 1 metric values for var 1
+		dfcolumns = c(dfcolumns,'B1')
+		dfcolours[ctr] = ifelse(c((dfall[variable==vars[1],][,8]) == 2),TableWinColour(),TableBenchmarkColour())
+	}
+	if(!all(is.na(dfall[variable==vars[1],][,6]))){
+		ctr = ctr + 1
+		dffinal[ctr] = signif(dfall[variable==vars[1],][,6],sfig) # benchmark 2 metric values for var 1
+		dfcolumns = c(dfcolumns,'B2')
+		dfcolours[ctr] = ifelse(c((dfall[variable==vars[1],][,8]) == 3),TableWinColour(),TableBenchmarkColour())
+	}
+	if(!all(is.na(dfall[variable==vars[1],][,7]))){
+		ctr = ctr + 1
+		dffinal[ctr] = signif(dfall[variable==vars[1],][,7],sfig) # benchmark 3 metric values for var 1
+		dfcolumns = c(dfcolumns,'B3')
+		dfcolours[ctr] = ifelse(c((dfall[variable==vars[1],][,8]) == 4),TableWinColour(),TableBenchmarkColour())
+	}
+	vctr[1] = ctr
+	if(length(vars) > 1){
+		for(v in 2:length(vars)){
+			ctr = ctr + 1
+			dffinal[ctr] = signif(dfall[variable==vars[v],][,4],sfig)
+			dfcolumns = c(dfcolumns,paste('Mod',vars[v]))
+			dfcolours[ctr] = ifelse(c((dfall[variable==vars[v],][,8]) == 1),TableWinColour(),TableModelColour())
+			if(!all(is.na(dfall[variable==vars[v],][,5]))){
+				ctr = ctr + 1
+				dffinal[ctr] = signif(dfall[variable==vars[v],][,5],sfig)
+				dfcolumns = c(dfcolumns,'B1')
+				dfcolours[ctr] = ifelse(c((dfall[variable==vars[v],][,8]) == 2),TableWinColour(),TableBenchmarkColour())
+			}
+			if(!all(is.na(dfall[variable==vars[v],][,6]))){
+				ctr = ctr + 1
+				dffinal[ctr] = signif(dfall[variable==vars[v],][,6],sfig)
+				dfcolumns = c(dfcolumns,'B2')
+				dfcolours[ctr] = ifelse(c((dfall[variable==vars[v],][,8]) == 3),TableWinColour(),TableBenchmarkColour())
+			}
+			if(!all(is.na(dfall[variable==vars[v],][,7]))){
+				ctr = ctr + 1
+				dffinal[ctr] = signif(dfall[variable==vars[v],][,7],sfig)
+				dfcolumns = c(dfcolumns,'B3')
+				dfcolours[ctr] = ifelse(c((dfall[variable==vars[v],][,8]) == 4),TableWinColour(),TableBenchmarkColour())
+			}
+			vctr[v] = ctr
+			# Decide whether table should be split:
+			if((v==halfvars) && (ctr > 5)){
+				halftable = ctr + 1	
+				split_table = TRUE
+			}
+		}
+	}
+	# Add column and row names to data frame:
+	colnames(dffinal) = dfcolumns
+	rownames(dffinal) = dfrownames
+	# Produce table grob...
+	outfile = setOutput('ModelAnalysis') # First set output
+	# Create blank plot with size:
+	plot(1:10, axes = FALSE, xlab = "", ylab = "", type = "n")
+	# Title and subtitle:
+	if(BenchInfo$number>0){
+		title_line2 = paste('Model: ',moname,'    Benchmarks: ',paste(benchnames,collapse=', '))
+	}else{
+		title_line2 = paste('Model:',moname)
+	}
+	title(main=paste(obsname,'metric summary'))
+	mtext(title_line2,cex=1.1,col='blue4')
+	# Draw table:
+	if(split_table){
+		addtable2plot(xpos[1],ypos[1],dffinal[1:vctr[1]],bty='o',hlines=TRUE,vlines=TRUE,cex=tablecex,
+			display.rownames=TRUE,bg=as.matrix(dfcolours)[,1:vctr[1]],xpad=xpadding,ypad=ypadding)
+		for(v in 2:3){
+			addtable2plot(xpos[v],ypos[1],dffinal[(vctr[v-1]+1):(vctr[v])],bty='o',hlines=TRUE,vlines=TRUE,cex=tablecex,
+				display.rownames=FALSE,bg=as.matrix(dfcolours)[,(vctr[v-1]+1):(vctr[v])],xpad=xpadding,ypad=ypadding)	
+		}
+		addtable2plot(xpos[1],ypos[2],dffinal[(vctr[3]+1):vctr[4]],bty='o',hlines=TRUE,vlines=TRUE,cex=tablecex,
+			display.rownames=TRUE,bg=as.matrix(dfcolours)[,(vctr[3]+1):vctr[4]],xpad=xpadding,ypad=ypadding)
+		for(v in 5:length(vars)){
+			addtable2plot(xpos[v%%3],ypos[2],dffinal[(vctr[v-1]+1):(vctr[v])],bty='o',hlines=TRUE,vlines=TRUE,cex=tablecex,
+				display.rownames=FALSE,bg=as.matrix(dfcolours)[,(vctr[v-1]+1):(vctr[v])],xpad=xpadding,ypad=ypadding)	
+		}
+	}else{
+		addtable2plot(xpos[1],7,dffinal[1:vctr[1]],bty='o',hlines=TRUE,vlines=TRUE,cex=tablecex,
+			display.rownames=TRUE,bg=as.matrix(dfcolours)[,1:vctr[1]],xpad=xpadding,ypad=ypadding)
+		for(v in 2:length(vars)){
+			addtable2plot(xpos[v],7,dffinal[(vctr[v-1]+1):(vctr[v])],bty='o',hlines=TRUE,vlines=TRUE,cex=tablecex,
+				display.rownames=FALSE,bg=as.matrix(dfcolours)[,(vctr[v-1]+1):(vctr[v])],xpad=xpadding,ypad=ypadding)	
+		}
+		
+		addtable2plot(-0.2,7,dffinal,bty='o',hlines=TRUE,vlines=TRUE,display.rownames=TRUE,cex=tablecex,
+			bg=as.matrix(dfcolours),xpad=xpadding,ypad=ypadding)
+	}
+		
+	# Return analysis result by adding to outinfo list of analyses results.		
+	result = list(type=thisanalysistype,filename=paste(getwd(),outfile,sep = "/"),
+		mimetype="image/png",analysistype=thisanalysistype)
+	outinfo[[(length(outinfo)+1)]] = result
+	return(outinfo)
+}
+
+AlignBenchmarks = function(UserBenchNames,AnalysisBenchNames,CurrentUserBench,BenchValues){
+	# For a given user nominated benchmark - CurrentUserBench - this function determines what the
+	# metric value of the benchmark is. In particular, when a user-nominated benchmark is not utlised for
+	# an analysis (variable missing / error), the user-nominated list of benchmarks becomes out of sync
+	# with the analysis list of benchmarks. This function ensures that the user-nominated benchmark metric 
+	# value is allocated appropriately in such cases.
+	
+	#print(UserBenchNames)
+	#print(AnalysisBenchNames)
+	#print(CurrentUserBench)
+	#print(BenchValues)
+	
+	if(CurrentUserBench == 1){
+		# Either user nominate benchamrk is analysis benchmark 1, or it failed:
+		if(is.na(AnalysisBenchNames[1])){ # i.e. all analyses failed (no 1st benchmark in analysis)
+			result = NA
+		}else if(UserBenchNames[1] == AnalysisBenchNames[1]){
+			result = BenchValues$bench1
+		}else{ # just the first user nominated benchmark failed
+			result = NA
+		}
+	}else if(CurrentUserBench == 2){
+		if(is.na(AnalysisBenchNames[1])){ # i.e. all analyses failed (no 1st benchmark in analysis)
+			result = NA
+		}else if(UserBenchNames[2] == AnalysisBenchNames[1]){
+			result = BenchValues$bench1
+		}else if(is.na(AnalysisBenchNames[2])){ # i.e. no 2nd benchmark in analysis
+			result = NA
+		}else if(UserBenchNames[2] == AnalysisBenchNames[2]){
+			result = BenchValues$bench2
+		}else{ # the 2nd user nominated benchmark failed
+			result = NA
+		}
+	}else if(CurrentUserBench == 3){
+		if(is.na(AnalysisBenchNames[1])){ # i.e. all analyses failed (no 1st benchmark in analysis)
+			result = NA
+		}else if(UserBenchNames[3] == AnalysisBenchNames[1]){
+			result = BenchValues$bench1
+		}else if(is.na(AnalysisBenchNames[2])){ # i.e. no 2nd benchmark in analysis
+			result = NA
+		}else if(UserBenchNames[3] == AnalysisBenchNames[2]){
+			result = BenchValues$bench2
+		}else if(is.na(AnalysisBenchNames[3])){ # i.e. no 3rd benchmark in analysis
+			result = NA	
+		}else if(UserBenchNames[3] == AnalysisBenchNames[3]){
+			result = BenchValues$bench3
+		}
+	}
+	
+	if(is.null(result)){ # e.g. if analysis does not return a value for a benchmark
+		result = NA	
+	}
+	
+	return(result)
+}
+
+BestInMetric = function(name,model,bench1,bench2,bench3){
+	# Returns an integer that indicates whether the submitting LSM (1), 
+	# the first (2), second (3), or third (4) nominated benchmark was
+	# best in a particular metric.
+	# If there are no benchmark values, or no metric value is given, there is no winner and 0 is returned
+	if(is.na(bench1) | grepl('nometric',tolower(name))){
+		return(0)
+	}
+	if(grepl('nme',tolower(name)) | grepl('rmse',tolower(name))){
+		result = which.min(c(model,bench1,bench2,bench3))
+	}else if(grepl('grad',tolower(name)) | grepl('gradient',tolower(name))){
+		result = which.min(abs(1-c(model,bench1,bench2,bench3)))
+	}else if(grepl('pdf',tolower(name)) | grepl('overlap',tolower(name)) | 
+		grepl('correlation',tolower(name)) | grepl('cor',tolower(name))){
+		result = which.max(c(model,bench1,bench2,bench3))
+	}else if(grepl('int',tolower(name)) | grepl('intercept',tolower(name))){
+		result = which.min(abs(c(model,bench1,bench2,bench3)))
+	}else if(grepl('bias',tolower(name))){
+		result = which.min(abs(c(model,bench1,bench2,bench3)))
+	}
+	return(result)
+}
+
+PadRowNames = function(initnames){
+	len = max(nchar(initnames))
+	result = initnames
+	for(n in 1:length(initnames)){
+		addchar = len - nchar(initnames[n])
+		if(addchar != 0){
+			result[n] = paste(' ',initnames[n],sep='')
+			if(addchar > 1){
+				for(a in 1:length(addchar)){
+					result[n] = paste(' ',result[n],sep='')
+				}
+			}
+		}
+	}
+	return(result)
+}
+
+# make metrics into matrix / data frame
+
+# use grid.table with data frame
+
+# or heatmap.2
+
+# or addtable2plot
+
+# or color2D.matplot
\ No newline at end of file
diff --git a/pals/R/NetcdfCompatibility.R b/pals/R/NetcdfCompatibility.R
new file mode 100644
index 0000000..a83cc0d
--- /dev/null
+++ b/pals/R/NetcdfCompatibility.R
@@ -0,0 +1,377 @@
+# NetcdfCompatibility.R
+#
+# Functions that check netcdf compatibility, obs-model compatibility
+# Gab Abramowitz UNSW 2014 (palshelp at gmail dot com)
+#
+CanAnalysisProceed = function(obs,model){
+	# Checks obs, model variables were found, timing is appropriate, and spatial grids match
+		
+	# Check for obs or model aren't missing variable data:
+	readcheck = CheckDataRead(obs$err,obs$errtext,model$err,model$errtext)
+	# Don't proceed and report error if there's an issue:		
+	if(readcheck$err){
+		return(readcheck)
+	}	
+	# Check model, obs timing consistency
+	tcheck = CheckTiming(obs$timing,model$timing)	
+	# Don't proceed and report error if there's an issue:	
+	if(tcheck$err){
+		return(tcheck)
+	}
+	# Check model, obs grid consistency
+	gcheck = CheckGrids(obs$grid,model$grid)	
+	# Don't proceed and report error if there's an issue:	
+	if(gcheck$err){
+		return(gcheck)
+	}
+	proceed = list(err=FALSE)
+	return(proceed)
+}
+
+GetGrid = function(fid){
+	# Gets spatial grid information from a netcdf file
+	errtext='ok'
+	# Fetch accepted names for latitude:
+	lat_det = GetVariableDetails('lat')
+	# Check an acceptable latitude variable exists:
+	exists_lat = AnyNcvarExists(fid,lat_det[[1]][['Name']])
+	if(! exists_lat$var){
+		errtext = paste('Could not find latitude variable in',stripFilename(fid$filename))
+		latlon = list(err=TRUE,errtext=errtext)
+		return(latlon)
+	}
+	# Now fetch latitude information:
+	if(exists_lat$dimvar){
+		# either as a dimension variable, if that was the found acceptable name:
+		lat = fid$dim[[exists_lat$dimnum]]$vals
+		latlen = fid$dim[[exists_lat$dimnum]]$len
+	}else{
+		# or as a normal variable:
+		if(fid$var[[exists_lat$varnum]]$ndims != 1){ # If not a 1D latitude variable
+			# Check that it's still a lat-lon grid:
+			lat = ncvar_get(fid,lat_det[[1]][['Name']][exists_lat$index]) # get lat data
+			latlen = length(unique(as.vector(lat))) # i.e. the number of unique latitude values
+			if(any(fid$var[exists_lat$index]$size != latlen)){
+				# i.e. number of unique lat vals is not equal to the length of either of the dims
+				# on which latitude depends
+				errtext = 'PALS cannot yet cope with non rectilinear grids yet.'	
+				latlon = list(err=TRUE,errtext=errtext)
+				return(latlon)
+			}
+		}
+	}
+	
+	# Repeat the process for longitude:
+	# Fetch accepted names for longitude:
+	lon_det = GetVariableDetails('lon')
+	# Check an acceptable longitude variable exists:
+	exists_lon = AnyNcvarExists(fid,lon_det[[1]][['Name']])
+	if(! exists_lon$var){
+		errtext = paste('Could not find longitude variable in',stripFilename(fid$filename))
+		latlon = list(err=TRUE,errtext=errtext)
+		return(latlon)
+	}
+	# Now fetch longitude information:
+	if(exists_lon$dimvar){
+		# either as a dimension variable, if that was the found acceptable name:
+		lon = fid$dim[[exists_lon$dimnum]]$vals
+		lonlen = fid$dim[[exists_lon$dimnum]]$len
+	}else{
+		# or as a normal variable:
+		if(fid$var[[exists_lon$varnum]]$ndims != 1){ # If not a 1D latitude variable
+			# Check that it's still a lat-lon grid:
+			lon = ncvar_get(fid,lon_det[[1]][['Name']][exists_lon$index]) # get lon data
+			lonlen = length(unique(as.vector(lon))) # i.e. the number of unique latitude values
+			if(any(fid$var[exists_lat$index]$size != latlen)){
+				# i.e. number of unique lat vals is not equal to the length of either of the dims
+				# on which latitude depends
+				errtext = 'PALS cannot yet cope with non lat-lon grids yet.'	
+				latlon = list(err=TRUE,errtext=errtext)
+				return(latlon)
+			}
+		}
+	}
+	# Return result:
+	latlon = list(err=FALSE,errtext=errtext,lat=lat,lon=lon,latlen=latlen,lonlen=lonlen)
+	return(latlon)
+}
+
+CheckGrids = function(obs,mod){
+	# Checks compatibility of observed and model output spatial grids
+	errtext = 'ok'
+	err = FALSE
+#	if(any(obs$lat != mod$lat) | any(obs$lon != mod$lon)){
+#		err=TRUE
+#		errtext = 'Spatial grids incompatible between obs and modelled data.'	
+#	}
+	result = list(err = err,errtext = errtext)
+	return(result)
+}
+
+CheckDataRead = function(obserr,obserrtext,moderr,moderrtext){
+	# Simply reports whether there was a read error from either
+	# obs or model output reading (no file, appropriate variable etc)
+	errtext = 'ok'
+	err = FALSE
+	if(obserr){
+		err = TRUE
+		errtext = obserrtext
+	}else if(moderr){
+		err = TRUE
+		errtext = moderrtext
+	}	
+	result = list(err = err,errtext = errtext)
+	return(result)
+}
+
+NcvarExists = function(fid,varname){
+	# Checks that variable exists in netcdf file, and additionally 
+	# checks whether quality control couterpart variable exists:
+	exists_var = FALSE
+	exists_qc = FALSE
+	varnum = NULL
+	dimnum = NULL
+	dimvar = FALSE
+	for (v in 1:fid$nvars){ # Search through all variables in netcdf file
+		if(fid$var[[v]]$name==varname){
+			exists_var=TRUE
+			varnum = v
+			if(exists_var & exists_qc) {break}
+		}else if(fid$var[[v]]$name==paste(varname,'_qc',sep='')){
+			exists_qc=TRUE
+			if(exists_var & exists_qc) {break}	
+		}
+	}
+	# If not found as a variable, it may be a dimension variable:
+	if(! exists_var){
+		for(d in 1:length(fid$dim)){
+			# If a dimension name & dimension variable
+			if(fid$dim[[d]]$name == varname){
+				if(fid$dim[[d]]$dimvarid$id != -1){# i.e. dim var exists
+					exists_var=TRUE
+					dimvar = TRUE
+					dimnum = d
+				}
+			}
+		}
+	}
+	ncvexists = list(var=exists_var,qc=exists_qc,dimvar=dimvar,dimnum=dimnum,varnum=varnum)
+	return(ncvexists)	
+}
+
+AnyNcvarExists = function(fid,varnames){
+	# Checks whether a variable exists as any of its alternative standard names.
+	exists_var = FALSE
+	exists_qc = FALSE
+	vindex = 0
+	for(v in 1:length(varnames)){
+		exists = NcvarExists(fid,varnames[v])
+		if(exists$var){
+			exists_var=TRUE
+			exists_qc = exists$qc
+			vindex = v
+			break
+		}	
+	}
+	ncvexists = list(var=exists_var,qc=exists_qc,index=vindex,dimvar=exists$dimvar,
+		dimnum=exists$dimnum,varnum=exists$varnum)
+	return(ncvexists)	
+}
+
+CheckNcvarUnits = function(fid,vname,variable,file){
+	# Check whether units exist for a variable, and if so, whether
+	# they correspond to any standard version of units text:
+	# Get units for variable
+	errtext='ok'
+	mvunits=ncatt_get(fid,vname,'units')
+	if(! mvunits$hasatt){
+		errtext = paste('Variable',vname,'in file',
+			file,'does not have a units attribute.')
+		units = list(errtext = errtext, err=TRUE)
+		return(units)
+	}
+	UnitsExist = TRUE
+	UnitsMatch = FALSE
+	for (u in 1:length(variable$UnitsName)){
+		if(mvunits$value==variable$UnitsName[u]){ # i.e. found units match
+			UnitsMatch = TRUE
+			# Set units adjustments appropriately:
+			multiplier = variable$Multiplier[u]
+			addition = variable$Addition[u]
+			break # out of units name for loop
+		}
+	}
+	if( ! UnitsMatch){ # i.e. didn't recognise variable units
+		units=mvunits$value
+		errtext = paste('Did not recognise units',mvunits$value,'for',
+			fid$var[[v]]$name,'in model output file',file)
+		units = list(errtext = errtext, err=TRUE)
+		return(units)	
+	}
+	units = list(exist = UnitsExist, match=UnitsMatch, value = mvunits$value,
+		multiplier=multiplier, addition=addition,err=FALSE,errtext=errtext)
+	return(units)
+}
+
+GetVariableDetails = function(request_names){
+	variable_list = list()
+	for(v in 1:length(request_names)){
+		var_details = list()
+		# Dimensional variables
+		if(request_names[v] == 'lat'){
+			var_details[['Name']] = c('y','lat','latitude','Lat','Latitude')
+			var_details[['UnitsName']] = c('degrees_north')
+			var_details[['Multiplier']] = c(1)
+			var_details[['Addition']] = c(0)
+			var_details[['UnitsText']] = ''
+			var_details[['PlotName']] = 'Latitude'
+			var_details[['Range']] = c(-90,180)
+		}else if(request_names[v] == 'lon'){
+			var_details[['Name']] = c('x','lon','longitude','Lon','Longitude')
+			var_details[['UnitsName']] = c('degrees_east')
+			var_details[['Multiplier']] = c(1)
+			var_details[['Addition']] = c(0)
+			var_details[['UnitsText']] = ''
+			var_details[['PlotName']] = 'Longitude'
+			var_details[['Range']] = c(-180,360)
+		# Met variables
+		}else if(request_names[v] == 'SWdown'){
+			var_details[['Name']] = c('SWdown')
+			var_details[['UnitsName']] = c('W/m2','W/m^2','Wm^{-2}','wm2','watt/m^2')
+			var_details[['Multiplier']] = c(1,1,1,1,1)
+			var_details[['Addition']] = c(0,0,0,0,0)
+			var_details[['UnitsText']] = 'W/'~m^{2}
+			var_details[['PlotName']] = 'Downward shortwave radiation'
+			var_details[['Range']] = c(0,1360)
+		}else if(request_names[v] == 'LWdown'){
+			var_details[['Name']] = c('LWdown')
+			var_details[['UnitsName']] = c('W/m2','W/m^2','Wm^{-2}','wm2','watt/m^2')
+			var_details[['Multiplier']] = c(1,1,1,1,1)
+			var_details[['Addition']] = c(0,0,0,0,0)
+			var_details[['UnitsText']] = 'W/'~m^{2}
+			var_details[['PlotName']] = 'Downward longwave radiation'
+			var_details[['Range']] = c(0,750)
+		}else if(request_names[v] == 'Tair'){
+			var_details[['Name']] = c('Tair')
+			var_details[['UnitsName']] = c('K')
+			var_details[['Multiplier']] = c(1)
+			var_details[['Addition']] = c(0)
+			var_details[['UnitsText']] = 'K'
+			var_details[['PlotName']] = 'Surface air temperature'
+			var_details[['Range']] = c(200,333)
+		}else if(request_names[v] == 'Qair'){
+			var_details[['Name']] = c('Qair')
+			var_details[['UnitsName']] = c('kg/kg','g/g')
+			var_details[['Multiplier']] = c(1,1)
+			var_details[['Addition']] = c(0,0)
+			var_details[['UnitsText']] = 'kg/kg'
+			var_details[['PlotName']] = 'Specific humidity'
+			var_details[['Range']] = c(0,0.1)
+		}else if(request_names[v] == 'PSurf'){
+			var_details[['Name']] = c('PSurf')
+			var_details[['UnitsName']] = c('Pa','pa','hPa')
+			var_details[['Multiplier']] = c(1,1,100) # when using first listed units
+			var_details[['Addition']] = c(0,0,0)
+			var_details[['UnitsText']] = 'Pa'
+			var_details[['PlotName']] = 'Surface air pressure'
+			var_details[['Range']] = c(50000,110000) # when using first listed units
+		}else if(request_names[v] == 'Rainf'){
+			var_details[['Name']] = c('Rainf')
+			var_details[['UnitsName']] = c('mm/s','kg/m^2/s', 'kg/m^2s', 'mm/day')
+			var_details[['Multiplier']] = c(1,1,1,86400)
+			var_details[['Addition']] = c(0,0,0,0)
+			var_details[['UnitsText']] = 'mm/s'
+			var_details[['PlotName']] = 'Precipitation'
+			var_details[['Range']] = c(0,0.05) # when using first listed units
+		}else if(request_names[v] == 'Wind'){
+			var_details[['Name']] = c('Wind')
+			var_details[['UnitsName']] = c('m/s')
+			var_details[['Multiplier']] = c(1)
+			var_details[['Addition']] = c(0)
+			var_details[['UnitsText']] = 'm/s'
+			var_details[['PlotName']] = 'Windspeed'
+			var_details[['Range']] = c(0,75)
+		# Flux variables
+		}else if(request_names[v] == 'Qh'){
+			var_details[['Name']] = c('Qh','FSH') # FSH->CLM
+			var_details[['UnitsName']] = c('W/m2','W/m^2','Wm^{-2}','wm2','watt/m^2')
+			var_details[['Multiplier']] = c(1,1,1,1,1)
+			var_details[['Addition']] = c(0,0,0,0,0)
+			var_details[['UnitsText']] = 'W/'~m^{2}
+			var_details[['PlotName']] = 'Sensible heat flux'
+			var_details[['Range']] = c(-1000,1000)
+		}else if(request_names[v] == 'Qle'){
+			var_details[['Name']] = c('Qle','FCEV') # NB: FCEV is only PART of Qle for CLM
+			var_details[['UnitsName']] = c('W/m2','W/m^2','Wm^{-2}','wm2','watt/m^2')
+			var_details[['Multiplier']] = c(1,1,1,1,1)
+			var_details[['Addition']] = c(0,0,0,0,0)
+			var_details[['UnitsText']] = 'W/'~m^{2}
+			var_details[['PlotName']] = 'Latent heat flux'
+			var_details[['Range']] = c(-1000,1000)
+		}else if(request_names[v] == 'Rnet'){
+			var_details[['Name']] = c('Rnet')
+			var_details[['UnitsName']] = c('W/m2','W/m^2','Wm^{-2}','wm2','watt/m^2')
+			var_details[['Multiplier']] = c(1,1,1,1,1)
+			var_details[['Addition']] = c(0,0,0,0,0)
+			var_details[['UnitsText']] = 'W/'~m^{2}
+			var_details[['PlotName']] = 'Net radiation'
+			var_details[['Range']] =  c(-1000,1360)
+		}else if(request_names[v] == 'SWnet'){
+			var_details[['Name']] = c('SWnet')
+			var_details[['UnitsName']] = c('W/m2','W/m^2','Wm^{-2}','wm2','watt/m^2')
+			var_details[['Multiplier']] = c(1,1,1,1,1)
+			var_details[['Addition']] = c(0,0,0,0,0)
+			var_details[['UnitsText']] = 'W/'~m^{2}
+			var_details[['PlotName']] = 'Net shortwave radiation'
+			var_details[['Range']] =  c(0,1360)
+		}else if(request_names[v] == 'LWnet'){
+			var_details[['Name']] = c('LWnet')
+			var_details[['UnitsName']] = c('W/m2','W/m^2','Wm^{-2}','wm2','watt/m^2')
+			var_details[['Multiplier']] = c(1,1,1,1,1)
+			var_details[['Addition']] = c(0,0,0,0,0)
+			var_details[['UnitsText']] = 'W/'~m^{2}
+			var_details[['PlotName']] = 'Net longwave radiation'
+			var_details[['Range']] = c(0,750)
+		}else if(request_names[v] == 'Qg'){
+			var_details[['Name']] = c('Qg')
+			var_details[['UnitsName']] = c('W/m2','W/m^2','Wm^{-2}','wm2','watt/m^2')
+			var_details[['Multiplier']] = c(1,1,1,1,1)
+			var_details[['Addition']] = c(0,0,0,0,0)
+			var_details[['UnitsText']] = 'W/'~m^{2}
+			var_details[['PlotName']] = 'Ground heat flux'
+			var_details[['Range']] = c(-1000,1000)
+		}else if(request_names[v] == 'NEE'){
+			var_details[['Name']] = c('NEE','FCO2') # FCO2->CLM
+			var_details[['UnitsName']] = c('umol/m2/s','mumol/m2/s','umol/m^2/s',
+				'umol/m2s','gC/m^2/s','gC/m2/s')
+			var_details[['Multiplier']] = c(1,1,1,1,1/(1.201E-5),1/(1.201E-5))
+			var_details[['Addition']] = c(0,0,0,0,0)
+			var_details[['UnitsText']] = ~mu~'mol/'~m^{2}~'/s' # default PALS units for plots
+			var_details[['PlotName']] = 'Net ecosystem exchange'
+			var_details[['Range']] = c(-100,100)
+		}else if(request_names[v] == 'GPP'){
+			var_details[['Name']] = c('GPP')
+			var_details[['UnitsName']] = c('umol/m2/s','mumol/m2/s','umol/m^2/s',
+				'umol/m2s','gC/m^2/s','gC/m2/s')
+			var_details[['Multiplier']] = c(1,1,1,1,1/(1.201E-5),1/(1.201E-5))
+			var_details[['Addition']] = c(0,0,0,0,0)
+			var_details[['UnitsText']] = ~mu~'mol/'~m^{2}~'/s' # default PALS units for plots
+			var_details[['PlotName']] = 'Gross primary production'
+			var_details[['Range']] = c(-100,100)
+		}
+		variable_list[[v]] = var_details
+	}
+	return(variable_list)
+}
+GetVariableIndex = function(vars,varname){
+	# Returns the index of a requested variable in a
+	# vars variable from a call to GetVariableDetails
+	vindex=NA
+	for(v in 1:length(vars)){
+		if(vars[[v]][['Name']][1] == varname){
+			vindex = v
+			break
+		}
+	}
+	return(vindex)	
+}
\ No newline at end of file
diff --git a/pals/R/SpatialFunctions.R b/pals/R/SpatialFunctions.R
index 4e9694e..8b7d720 100644
--- a/pals/R/SpatialFunctions.R
+++ b/pals/R/SpatialFunctions.R
@@ -10,13 +10,13 @@ GridAreaVector = function(obs,mcount){
 	cat('Creating gridcell surface area vector')
 	if(sum(obs$missingmask) == length(obs$missingmask)){ # no missing data in 20C obs
 		cat(' without concat')
-		surf_areas=c(aperm(array(area$cell,dim=c(length(obs$lon),length(obs$lat),mcount)),c(3,2,1)))
+		surf_areas=c(aperm(array(area$cell,dim=c(length(obs$grid$lon),length(obs$grid$lat),mcount)),c(3,2,1)))
 	}else{
 		gdctr=0 # initialise grid cell counter
 		surf_areas = c()
 		#gdctr_ref = matrix(NA,obs$gdctr,2)
-		for(i in 1:length(obs$lon)){
-			for(j in 1:length(obs$lat)){
+		for(i in 1:length(obs$grid$lon)){
+			for(j in 1:length(obs$grid$lat)){
 				# First check if there is enough obs data for result:
 				if(obs$missingmask[i,j]){ # where there IS enough data
 					gdctr = gdctr + 1 # increment grid cell counter
@@ -37,14 +37,14 @@ GridAreaVector = function(obs,mcount){
 #
 # Calculates surface area of Earth and grid cells:
 EarthArea = function(obs){
-	cellwidth =  obs$lat[2] - obs$lat[1] # assume n degree by n degree
+	cellwidth =  obs$grid$lat[2] - obs$grid$lat[1] # assume n degree by n degree
 	earthradius=6371
 	earth = 4*pi*earthradius^2
-	cell = matrix(NA,length(obs$lon),length(obs$lat))
-	for(i in 1:length(obs$lon)){
-		for(j in 1:length(obs$lat)){
+	cell = matrix(NA,length(obs$grid$lon),length(obs$grid$lat))
+	for(i in 1:length(obs$grid$lon)){
+		for(j in 1:length(obs$grid$lat)){
 			cell[i,j] = earthradius^2 * cellwidth/360*2*pi * 
-				(sin((obs$lat[j]+cellwidth/2)/90*pi/2) - sin((obs$lat[j]-cellwidth/2)/90*pi/2))
+				(sin((obs$grid$lat[j]+cellwidth/2)/90*pi/2) - sin((obs$grid$lat[j]-cellwidth/2)/90*pi/2))
 		}
 	}
 	
diff --git a/pals/R/UtilityFunctions.R b/pals/R/UtilityFunctions.R
index 65656cf..f6f9d8d 100644
--- a/pals/R/UtilityFunctions.R
+++ b/pals/R/UtilityFunctions.R
@@ -4,93 +4,6 @@
 #
 # Gab Abramowitz UNSW 2014 (palshelp at gmail dot com)
 #
-LineColours = function() {
-	# For line plots:
-	plotcolours=c('black','blue2','indianred3','gold2','yellowgreen')
-}
-ChooseColours = function(range,variablename,plottype,diffthreshold=NULL){
-	# Returns a colour range for gridded plots
-	library(colorRamps)
-	
-	# Full / most range:
-	red2blue = colorRampPalette(c('red','orange','yellow','green','blue'))
-	yellow2purpleCool = colorRampPalette(c('yellow','green3','blue','darkorchid4'))
-	yellow2purpleWarm = colorRampPalette(c('yellow','red','magenta'))
-	purple2yellowWarm = colorRampPalette(c('magenta','red','yellow'))
-	iceblue2green = colorRampPalette(c('slategray1','midnightblue','blue','green3','green'))
-	green2iceblue = colorRampPalette(c('green','green3','blue','midnightblue','slategray1'))
-	
-	# Half range:
-	green2darkblue = colorRampPalette(c('green','green4','blue','midnightblue'))
-	darkblue2green = colorRampPalette(c('midnightblue','blue','green4','green'))
-	darkred2yellow = colorRampPalette(c('red4','red','orange','yellow'))
-	yellow2darkred = colorRampPalette(c('yellow','orange','red','red4'))
-	
-	# Small range:
-	yellow2red = colorRampPalette(c('yellow','red'))
-	red2yellow = colorRampPalette(c('red','yellow'))
-	green2blue = colorRampPalette(c('green','blue'))
-	blue2green = colorRampPalette(c('blue','green'))
-	
-	coolvars = c('Qle','Evap')
-	warmvars = c('Tair','Qh','Rnet','SWdown','SWnet')
-	colourres = 36 # approximately how many colours in a plot (will control size of white space if diff plot)
-	
-	# If no difference threshold has been specified, use 5%:
-	if(is.null(diffthreshold)){
-		diffthreshold = (range[2] - range[1]) / 20
-	}
-	
-	# Assess cases where colours for a difference plot are requested first:
-	if(plottype=='difference'){
-		# i.e. the plot will contain a zero that we want coloured white
-		# First check that we really do need a difference plot:
-		if(range[1] > (-1*diffthreshold)){
-			# Just use a positive scale
-			plottype = 'positive'
-		}else if(range[2]<diffthreshold){
-			# Just use a negative scale
-			plottype = 'negative'
-		}
-		# Find fraction of range below 0
-		lowfrac = abs(range[1]) / (abs(range[1]) + range[2])
-		lownum = floor(lowfrac * colourres)
-		# Find fraction of range above 0
-		highfrac = range[2] / (abs(range[1]) + range[2])
-		highnum = floor(highfrac * colourres)
-		# Decide colour range:
-		if(any(warmvars == variablename)){ # For variables warm colours when positive
-			if(lowfrac/highfrac > 2){ # most fo the range is below 0
-				colours = c(iceblue2green(lownum),'#FFFFFF','#FFFFFF',yellow2red(highnum))
-			}else if(lowfrac/highfrac < 1/2){ # most of the range is above 0
-				colours = c(blue2green(lownum),'#FFFFFF','#FFFFFF',yellow2purpleWarm(highnum))
-			}else{
-				colours = c(darkblue2green(lownum),'#FFFFFF','#FFFFFF',yellow2darkred(highnum))
-			}
-		}else if(any(coolvars == variablename)){ # For variables cool colours when positive
-			if(lowfrac/highfrac > 2){ # most fo the range is below 0
-				colours = c(purple2yellowWarm(lownum),'#FFFFFF','#FFFFFF',green2blue(highnum))
-			}else if(lowfrac/highfrac < 1/2){ # most of the range is above 0
-				colours = c(red2yellow(lownum),'#FFFFFF','#FFFFFF',green2iceblue(highnum))
-			}else{
-				colours = c(darkred2yellow(lownum),'#FFFFFF','#FFFFFF',green2darkblue(highnum))
-			}
-		}
-	}
-	
-	# Now assess cases where just a positive or negative scale is required:
-	if((plottype=='positive') && (any(coolvars == variablename))){
-		colours = yellow2purpleCool(colourres)
-	}else if((plottype=='positive') && (any(warmvars == variablename))){
-		colours = yellow2purpleWarm(colourres)
-	}else if((plottype=='negative') && (any(coolvars == variablename))){
-		colours = purple2yellowWarm(colourres)
-	}else if((plottype=='negative') && (any(warmvars == variablename))){
-		colours = iceblue2green(colourres)
-	}
-	
-	return(colours)
-}
 
 # Function for crashing semi-gracefully:
 CheckError = function(errtext,errcode='U1:'){
@@ -102,7 +15,37 @@ CheckError = function(errtext,errcode='U1:'){
 		cat(alltext,' ^ \n',file=stderr()); stop(alltext,call. = FALSE)
 	}
 }
-
+CheckIfAllFailed = function(outinfo){
+	# Checks if all requested analyses failed (e.g. in which case don't run summary table function)
+	allfail = TRUE
+	for(a in 1:length(outinfo)){
+		if(is.null(outinfo[[a]]$error)){
+			allfail = FALSE
+		}
+	}
+	return(allfail)
+}
+LegendText = function(data,plotobs=TRUE){
+	# Returns text vector of legend names for a plot.
+	# If no obs line in the plot (e.g. error plot), first index will be model, else obs
+	# If for some reason a benchmark failed (e.g. missing variable), colours are adjusted to make them 
+	# consistent across different plots 
+	legendtext = c()
+	if(plotobs){ # i.e. obs line will be part of the plot
+		legendtext[1] = data$obs$name
+	}
+	legendtext = c(legendtext, data$model$name)
+	if(data$bench$exist){
+		legendtext = c(legendtext, data$bench[[ data$bench$index[1] ]]$name)
+		if(data$bench$howmany == 2){
+			legendtext = c(legendtext, data$bench[[ data$bench$index[2] ]]$name)
+		}else if(data$bench$howmany == 3){
+			legendtext = c(legendtext, data$bench[[ data$bench$index[2] ]]$name)
+			legendtext = c(legendtext, data$bench[[ data$bench$index[3] ]]$name)
+		}
+	}
+	return(legendtext)
+}
 FindRangeViolation = function(varin,varrange){
 	offendingValue=0 # init
 	for(i in 1:length(varin)){
@@ -136,34 +79,39 @@ CheckVersionCompatibility = function(filepath1,filepath2){
 	}
 }
 
-NumberOfBenchmarks = function(bench,Bctr){
+BenchmarkInfo = function(BenchmarkFiles,Bctr){
 	# Determines the number of user nominated benchmarks in a call to an 
 	# experiment script, as well as the number of files associated with each.
 	# Bctr - total number of benchmark files
-	# bench - contains data for each file
-	if(Bctr == 0){
+	# BenchmarkFiles - contains data for each file
+	# nBench - number of user nominated benchmarks
+	# nBenchfiles - a list of vectors of file indices for each benchmark
+	if(Bctr == 0){ # no benchmark files sent by javascript
 		nBench = 0
-		benchfiles = NA
+		nBenchfiles = NA
+		benchnames = NA
 	}else{
 		nBench = 1
-		# Determine number of user nominated benchmarks:
+		# Determine number of user nominated benchmarks, and the name of each:
+		benchnames = c()
 		for(b in 1:Bctr){
-			nBench = max(nBench, as.integer(bench[[Bctr]]$number))
+			nBench = max(nBench, as.integer(BenchmarkFiles[[b]]$number))
+			benchnames[as.integer(BenchmarkFiles[[b]]$number)] = BenchmarkFiles[[b]]$name
 		}
 		# Store which files belong to which benchmark:
-		benchfiles = list()
+		nBenchfiles = list()
 		bexists = c(0)
 		for(b in 1:Bctr){
-			benchnumber = as.integer(bench[[b]]$number)
+			benchnumber = as.integer(BenchmarkFiles[[b]]$number)
 			if(any(bexists == benchnumber)){
-				benchfiles[[benchnumber]] = c(benchfiles[[benchnumber]] , b)
+				nBenchfiles[[benchnumber]] = c(nBenchfiles[[benchnumber]] , b)
 			}else{
-				benchfiles[[benchnumber]] = c(b)
+				nBenchfiles[[benchnumber]] = c(b)
 				bexists = c(bexists,benchnumber)
 			}
 		}
 	}
-	result = list(number = nBench, benchfiles=benchfiles)
+	result = list(number = nBench, benchfiles=nBenchfiles, names=benchnames)
 	return(result)
 }
 #
diff --git a/scripts/GlobalGSWP30.5Experiment.R b/scripts/GlobalGSWP30.5Experiment.R
new file mode 100644
index 0000000..56388b1
--- /dev/null
+++ b/scripts/GlobalGSWP30.5Experiment.R
@@ -0,0 +1,88 @@
+# Master script for Global 0.5 x 0.5 degree experiment based on GSWP3 forcing
+# Gab Abramowitz, UNSW, 2015 (palshelp at gmail dot com)
+
+library(pals)
+library(parallel)
+
+print(paste('ID:',input["_id"]))
+files <- input[["files"]]
+
+# Retrieve model output, forcing and evaluation data set and benchmark location and 
+# meta data from javascript input list: 
+ModelOutputFiles = list()
+ForcingDataSetFiles = list()
+EvalDataSetFiles = list()
+BenchmarkFiles = list()
+MOctr = 0
+FDSctr = 0
+EDSctr = 0
+Bctr = 0
+for (i in 1:(length(files))  ) {
+    file <- files[[i]]
+    if( file[['type']] == "ModelOutput" ) {
+    	MOctr = MOctr + 1
+        ModelOutputFiles[[MOctr]] = list(path=file[['path']],mimetype=file[['mimetype']],
+        	name=file[['name']])
+    }else if( (file[['type']] == "DataSet")) {
+    	EDSctr = EDSctr + 1
+        EvalDataSetFiles[[EDSctr]] = list(path=file[['path']],mimetype=file[['mimetype']],
+        	name=file[['name']])
+    }else if( file[['type']] == "Benchmark") {
+    	Bctr = Bctr + 1
+        BenchmarkFiles[[Bctr]] = list(path=file[['path']],mimetype=file[['mimetype']],
+        	name=file[['name']],number=file[['number']]) # user rank of benchmark
+    }
+}
+
+region = 'Global'
+
+# Nominate variables to analyse here (use ALMA standard names) - fetches
+# alternate names, units, units transformations etc:
+vars = GetVariableDetails(c('Qle'))
+
+# Analyses that can apply to any variable:
+genAnalysis = c('TimeMean','TimeSD','TimeRMSE','TimeCor') #'PDFall','PDF2D','Taylor')
+
+# Determine number of user-nominated benchmarks:
+nBench = BenchmarkInfo(BenchmarkFiles,Bctr)
+
+# Set up analysis data and analysis list so we can use lapply or parlapply:
+AnalysisList = list()
+
+# Load all variables from obs and model output
+for(v in 1:length(vars)){
+	obs = GetGLEAM_Global(vars[[v]],EvalDataSetFiles,force_interval='monthly')
+    model = GetModelOutput(vars[[v]],ModelOutputFiles)  
+    bench = GetBenchmarks(vars[[v]],BenchmarkFiles,nBench)
+	# Add those analyses that are equally applicable to any variable to analysis list:
+	for(a in 1:length(genAnalysis)){
+		analysis_number = (v-1)*length(genAnalysis) + a
+		AnalysisList[[analysis_number]] = list(vindex=v, type=genAnalysis[a])
+	}
+}
+
+# Create cluster:
+cl = makeCluster(getOption('cl.cores', detectCores()))
+
+	OutInfo = lapply(AnalysisList,DistributeGriddedAnalyses,vars=vars,
+		obs=obs,model=model,bench=bench,region=region,cl)
+#	OutInfo = parLapply(cl=cl,AnalysisList,DistributeGriddedAnalyses,vars=vars,
+#		obs=obs,model=model,bench=bench,region=region,NULL)
+
+stopCluster(cl) # stop cluster
+
+# Write outinfo to output list for javascript:
+output = list(files=OutInfo);
+
+for(i in 1: length(output[["files"]])){
+	cat('Output ',i,': \n')
+	cat('  type:',output[["files"]][[i]]$type,'\n')
+	if(!is.null(output[["files"]][[i]]$error)){
+		cat('  ERROR: ',output[["files"]][[i]]$error,'\n')
+	}else{
+		cat('  filename:',output[["files"]][[i]]$filename,'\n')
+		cat('  bench error:',output[["files"]][[i]]$bencherror,'\n')
+		cat('  first metric for model - ',output[["files"]][[i]]$metrics[[1]]$name,':',
+			output[["files"]][[i]]$metrics[[1]]$model_value,'\n')
+	}
+}
\ No newline at end of file
diff --git a/scripts/RegionalAus0.25Experiment.R b/scripts/RegionalAus0.25Experiment.R
index b5ee593..253ffb8 100644
--- a/scripts/RegionalAus0.25Experiment.R
+++ b/scripts/RegionalAus0.25Experiment.R
@@ -9,10 +9,10 @@ files <- input[["files"]]
 
 # Retrieve model output, forcing and evaluation data set and benchmark location and 
 # meta data from javascript input list: 
-ModelOutputs = list()
-ForcingDataSets = list()
-EvalDataSets = list()
-Benchmarks = list()
+ModelOutputFiles = list()
+ForcingDataSetFiles = list()
+EvalDataSetFiles = list()
+BenchmarkFiles = list()
 MOctr = 0
 FDSctr = 0
 EDSctr = 0
@@ -21,29 +21,20 @@ for (i in 1:(length(files))  ) {
     file <- files[[i]]
     if( file[['type']] == "ModelOutput" ) {
     	MOctr = MOctr + 1
-        ModelOutputs[[MOctr]] = list(path=file[['path']],mimetype=file[['mimetype']],
+        ModelOutputFiles[[MOctr]] = list(path=file[['path']],mimetype=file[['mimetype']],
         	name=file[['name']])
     }else if( (file[['type']] == "DataSet")) {
     	EDSctr = EDSctr + 1
-        EvalDataSets[[EDSctr]] = list(path=file[['path']],mimetype=file[['mimetype']],
+        EvalDataSetFiles[[EDSctr]] = list(path=file[['path']],mimetype=file[['mimetype']],
         	name=file[['name']])
     }else if( file[['type']] == "Benchmark") {
     	Bctr = Bctr + 1
-        Benchmarks[[Bctr]] = list(path=file[['path']],mimetype=file[['mimetype']],
+        BenchmarkFiles[[Bctr]] = list(path=file[['path']],mimetype=file[['mimetype']],
         	name=file[['name']],number=file[['number']]) # user rank of benchmark
     }
 }
-for(f in 1:MOctr){
-	print(paste("Model Output file: ",ModelOutputs[[f]][['path']]));
-}
-for(f in 1:EDSctr){
-	print(paste("Data Set file: ",EvalDataSets[[f]][['path']]));
-}
-if(Bctr !=0){
-	for(f in 1:Bctr){
-		print(paste("Bench file: ",Benchmarks[[f]][['path']]));
-	}
-}
+
+region = 'Australia'
 
 # Nominate variables to analyse here (use ALMA standard names) - fetches
 # alternate names, units, units transformations etc:
@@ -53,18 +44,16 @@ vars = GetVariableDetails(c('Qle'))
 genAnalysis = c('TimeMean','TimeSD','TimeRMSE','TimeCor') #'PDFall','PDF2D','Taylor')
 
 # Determine number of user-nominated benchmarks:
-nBench = NumberOfBenchmarks(Benchmarks,Bctr)
-
-cat('\nUser number of benchmarks:',nBench$number,'\n')
+nBench = BenchmarkInfo(BenchmarkFiles,Bctr)
 
 # Set up analysis data and analysis list so we can use lapply or parlapply:
 AnalysisList = list()
 
 # Load all variables from obs and model output
 for(v in 1:length(vars)){
-	obs = GetGLEAM_Aus(vars[[v]],EvalDataSets,force_interval='monthly')
-    model = GetModelOutput(vars[[v]],ModelOutputs)  
-    bench = GetBenchmarks(vars[[v]],Benchmarks,nBench)
+	obs = GetGLEAM_Aus(vars[[v]],EvalDataSetFiles,force_interval='monthly')
+    model = GetModelOutput(vars[[v]],ModelOutputFiles)  
+    bench = GetBenchmarks(vars[[v]],BenchmarkFiles,nBench)
     
 	# Add those analyses that are equally applicable to any variable to analysis list:
 	for(a in 1:length(genAnalysis)){
@@ -74,22 +63,14 @@ for(v in 1:length(vars)){
 }
 
 # Create cluster:
-#cl = makeCluster(getOption('cl.cores', detectCores()))
-#cl = makeCluster(getOption('cl.cores', 2))
+cl = makeCluster(getOption('cl.cores', detectCores()))
 
 	OutInfo = lapply(AnalysisList,DistributeGriddedAnalyses,vars=vars,
-		obs=obs,model=model,bench=bench)
+		obs=obs,model=model,bench=bench,region=region,cl)
 #	OutInfo = parLapply(cl=cl,AnalysisList,DistributeGriddedAnalyses,vars=vars,
-#		obs=obs,model=model,bench=bench)
-
-# Add multiple variable analysis to analysis list:
-# analysis_number = analysis_number + 1
-# AnalysisList[[analysis_number]] = list(vindex=0, type='EvapFrac')
-# analysis_number = analysis_number + 1
-# AnalysisList[[analysis_number]] = list(vindex=0, type='Conserve')
+#		obs=obs,model=model,bench=bench,region=region,NULL)
 
-# stop cluster
-#stopCluster(cl)
+stopCluster(cl) # Stop cluster
 
 # Write outinfo to output list for javascript:
 output = list(files=OutInfo);
diff --git a/scripts/SingleSiteExperiment.R b/scripts/SingleSiteExperiment.R
index da53f01..30bd2f4 100644
--- a/scripts/SingleSiteExperiment.R
+++ b/scripts/SingleSiteExperiment.R
@@ -9,10 +9,10 @@ files <- input[["files"]]
 
 # Retrieve model output, forcing and evaluation data set and benchmark location and 
 # meta data from javascript input list: 
-ModelOutputs = list()
-ForcingDataSets = list()
-EvalDataSets = list()
-Benchmarks = list()
+ModelOutputFiles = list()
+ForcingDataSetFiles = list()
+EvalDataSetFiles = list()
+BenchmarkFiles = list()
 MOctr = 0
 FDSctr = 0
 EDSctr = 0
@@ -21,31 +21,31 @@ for (i in 1:(length(files))  ) {
     file <- files[[i]]
     if( file[['type']] == "ModelOutput" ) {
     	MOctr = MOctr + 1
-        ModelOutputs[[MOctr]] = list(path=file[['path']],mimetype=file[['mimetype']],
+        ModelOutputFiles[[MOctr]] = list(path=file[['path']],mimetype=file[['mimetype']],
         	name=file[['name']])
     }else if( (file[['type']] == "DataSet")) {
     	EDSctr = EDSctr + 1
-        EvalDataSets[[EDSctr]] = list(path=file[['path']],mimetype=file[['mimetype']],
+        EvalDataSetFiles[[EDSctr]] = list(path=file[['path']],mimetype=file[['mimetype']],
         	name=file[['name']])
     }else if( file[['type']] == "Benchmark") {
     	Bctr = Bctr + 1
-        Benchmarks[[Bctr]] = list(path=file[['path']],mimetype=file[['mimetype']],
+        BenchmarkFiles[[Bctr]] = list(path=file[['path']],mimetype=file[['mimetype']],
         	name=file[['name']],number=file[['number']]) # user rank of benchmark
     }
 }
 
-print(paste("Model Output file: ",ModelOutputs[[1]][['path']]));
-print(paste("Data Set file: ",EvalDataSets[[1]][['path']]));
+#print(paste("Model Output file: ",ModelOutputFiles[[1]][['path']]));
+#print(paste("Data Set file: ",EvalDataSetFiles[[1]][['path']]));
 
 # Nominate variables to analyse here (use ALMA standard names) - fetches
 # alternate names, units, units transformations etc:
 vars = GetVariableDetails(c('Qle','Qh','NEE','Rnet','Qg','SWnet'))
 
 # Analyses that can apply to any variable:
-genAnalysis = c('AvWindow','Scatter','Timeseries','AnnualCycle','DiurnalCycle','PDF','Taylor')
+analyses = c('Timeseries','Taylor','AnnualCycle','DiurnalCycle','Scatter','PDF')
 
 # Determine number of user-nominated benchmarks:
-nBench = NumberOfBenchmarks(Benchmarks,Bctr)
+BenchInfo = BenchmarkInfo(BenchmarkFiles,Bctr)
 
 # Set up analysis data and analysis list so we can use lapply or parlapply:
 AnalysisData = list()
@@ -54,15 +54,15 @@ AnalysisList = list()
 # Load all variables from obs and model output
 analysis_number = 0
 for(v in 1:length(vars)){
-	obs = GetFluxnetVariable(vars[[v]],EvalDataSets[[1]])
-    model = GetModelOutput(vars[[v]],ModelOutputs)
-    bench = GetBenchmarks(vars[[v]],Benchmarks,nBench)
+	obs = GetFluxnetVariable(vars[[v]],EvalDataSetFiles[[1]])
+    model = GetModelOutput(vars[[v]],ModelOutputFiles)
+    bench = GetBenchmarks(vars[[v]],BenchmarkFiles,BenchInfo)
 	# Save model, obs, bench data for each variable:
 	AnalysisData[[v]] = list(obs=obs, model=model, bench = bench)
 	# Add those analyses that are equally applicable to any variable to analysis list:
-	for(a in 1:length(genAnalysis)){
-		analysis_number = (v-1)*length(genAnalysis) + a
-		AnalysisList[[analysis_number]] = list(vindex=v, type=genAnalysis[a])
+	for(a in 1:length(analyses)){
+		analysis_number = (v-1)*length(analyses) + a
+		AnalysisList[[analysis_number]] = list(vindex=v, type=analyses[a])
 	}
 }
 # Add multiple variable analysis to analysis list:
@@ -82,14 +82,15 @@ outinfo = parLapply(cl=cl,AnalysisList,DistributeSingleSiteAnalyses,data=Analysi
 # stop cluster
 stopCluster(cl)
 
+# Draw summary metric table here (adds to outinfo list):
+outinfo = MetricTableSingleSite(outinfo,BenchInfo)
+
 # Write outinfo to output list for javascript:
 output = list(files=outinfo);
 
-#check error propagation!
-
 for(i in 1: length(output[["files"]])){
 	cat('Output ',i,': \n')
-	cat('  type:',output[["files"]][[i]]$type,'\n')
+	cat('  type:',output[["files"]][[i]]$type)
 	cat('  filename:',output[["files"]][[i]]$filename,'\n')
 	cat('  bench error:',output[["files"]][[i]]$bencherror,'\n')
 	cat('  first metric for model - ',output[["files"]][[i]]$metrics[[1]]$name,':',
diff --git a/scripts/TestBench3_SingleSite.nc b/scripts/TestBench3_SingleSite.nc
new file mode 100755
index 0000000..b4a6a48
Binary files /dev/null and b/scripts/TestBench3_SingleSite.nc differ
diff --git a/scripts/TestGlobalExperiment.R b/scripts/TestGlobalExperiment.R
new file mode 100644
index 0000000..17bef9b
--- /dev/null
+++ b/scripts/TestGlobalExperiment.R
@@ -0,0 +1,9 @@
+library("RJSONIO")
+inputFile <- "TestInput_Global.json"
+input <- fromJSON(paste(readLines(inputFile), collapse=""));
+Rruntime = system.time(source("GlobalGSWP30.5Experiment.R"))
+print(paste('Time to run:',Rruntime[3]))
+output <- toJSON(output)
+fileConn<-file("output_Global.json")
+writeLines(output, fileConn)
+close(fileConn)
\ No newline at end of file
diff --git a/scripts/TestInput_Global.json b/scripts/TestInput_Global.json
new file mode 100644
index 0000000..9b4c9f3
--- /dev/null
+++ b/scripts/TestInput_Global.json
@@ -0,0 +1,57 @@
+{
+   "_id" : "mVJDHsjHSJWHK",
+   "files" : [
+       {
+           "type" : "ModelOutput",
+           "path" : "/Users/gab/data/global/CABLE_output/GSWP3/cable_output_2005.nc",
+           "mimetype" : "application/x-netcdf",
+	   "name" : "CABLEtestGlobal"
+       },
+       {
+           "type" : "ModelOutput",
+           "path" : "/Users/gab/data/global/CABLE_output/GSWP3/cable_output_2006.nc",
+           "mimetype" : "application/x-netcdf",
+	   "name" : "CABLEtestGlobal"
+       },
+       {   
+           "type" : "DataSet", 
+           "path" : "/Users/gab/data/global/GLEAM/GSWP0.5/GLEAM_ET_GSWP3_2005.nc",
+           "mimetype" : "application/x-netcdf",
+	       "name" : "GLEAM_GSWP3"
+       },
+       {   
+           "type" : "DataSet", 
+           "path" : "/Users/gab/data/global/GLEAM/GSWP0.5/GLEAM_ET_GSWP3_2006.nc",
+           "mimetype" : "application/x-netcdf",
+	       "name" : "GLEAM_GSWP3"
+       },
+       {   
+           "type" : "Benchmark", 
+           "path" : "/Users/gab/data/global/CABLE_output/GSWP3/cable_output_2000.nc",
+           "mimetype" : "application/x-netcdf",
+		   "number" : "1",
+	       "name" : "Cab0-1bench"
+       },
+       {   
+           "type" : "Benchmark", 
+           "path" : "/Users/gab/data/global/CABLE_output/GSWP3/cable_output_2001.nc",
+           "mimetype" : "application/x-netcdf",
+		   "number" : "1",
+	       "name" : "Cab0-1bench"
+       },
+       {   
+           "type" : "Benchmark", 
+           "path" : "/Users/gab/data/global/CABLE_output/GSWP3/cable_output_1993.nc",
+           "mimetype" : "application/x-netcdf",
+		   "number" : "2",
+	       "name" : "Cab93-4bench"
+       },
+       {   
+           "type" : "Benchmark", 
+           "path" : "/Users/gab/data/global/CABLE_output/GSWP3/cable_output_1994.nc",
+           "mimetype" : "application/x-netcdf",
+		   "number" : "2",
+	       "name" : "Cab93-4bench"
+       }
+   ]
+}
diff --git a/scripts/TestInput_SingleSite.json b/scripts/TestInput_SingleSite.json
index ee991f2..dcc4f37 100644
--- a/scripts/TestInput_SingleSite.json
+++ b/scripts/TestInput_SingleSite.json
@@ -32,6 +32,13 @@
            "mimetype" : "application/x-netcdf",
            "number" : "2",
 	       "name" : "PenmanMontElSaler2"
+       },
+       {
+	   "type" : "Benchmark", 
+           "path" : "TestBench3_SingleSite.nc",
+           "mimetype" : "application/x-netcdf",
+           "number" : "3",
+	       "name" : "CABLE_SLI_ElSaler2"
 	}
    ]
 }