plotpc <- function(x,i=1,n=11,bw=FALSE,cont=FALSE,...){

# Plot EOFs and Principal Components (PCs, time series) that are output of the
# eof function 
#
# Description:
#
#      Plots the principal component time series and the map of loadings 
#      of a given mode
#
# Usage:
#
#      plotpc(x,i,n=11,bw=FALSE,cont=FALSE,...)
#
# Input:
#
#      x: the output of the function eofs
#
#      i: the number of the mode to be plotted
#
#
#   cont: Logical. If TRUE adds contours to the image plot. Default is FALSE.
#
#      n: An odd number of colours/intervals for the colourbar. Default is 11.
#
#     bw: Logical. If TRUE produces black, gray and white plot. Default is FALSE.  
#
#    ...: Additional arguments passed to `contour'.
#
# Output:
#
#
# Authors:
#
#      Dag Johan Steinskog <dag.johan.steinskog@nersc.no> 10 June 2005
#      Caio Coelho <c.a.d.s.coelho@reading.ac.uk> 
#
# Example:
#
#      lon <- seq(0,360,5)
#      lat <- seq(-90,90,5)
#      t <- 365
#      x <- array(rnorm(length(lon)*length(lat)*t,sd=10),dim=c(length(lon),length(lat),t))
#      sss<-eof(lon,lat,x)
#      plotpc(sss,2,cont=T)  # plots the second mode 
#      plotpc(sss,10,bw=T,cont=T) # plots the tenth mode 


rg<-range(x$eof[,,i]*100, na.rm=T)

lowerlim <- rg[1]
upperlim <- rg[2]

maximum <- max(abs(c(lowerlim,upperlim)))

if(lowerlim<0) breaks <- seq(-maximum,maximum,(maximum-(-maximum))/n)
else breaks <- seq(lowerlim,upperlim,(upperlim-(lowerlim))/n)

if(!bw){

# check if range includes negative values to use appropriate colour scale
if (rg[1] <0) {
colours = bluered(seq(0,n-1,by=1), "linear", yellow =TRUE)
}
else {
colours = bluered(seq(0,n-1,by=1), "linear", white=0, invert=T)
}

}

else{
if (rg[1] <0) {
colours <- grey(seq(0, 1, length = length(breaks)-1))
colours <- c(colours[1:((n-1)/2)],rev(colours[(((n-1)/2)+1):n]))
}
else {
colours <- grey(seq(0, 1, length = length(breaks)-1))
}


}



layout(matrix(1:3, ncol = 1, nrow=3), heights = c(5,8,1))
par(mar = c(6, 4, 4, 0.5))

# Plot PC time series
plot(x$pcs[,i],type='l',xlab='Time',ylab=paste('PC',as.character(i)))
title(paste('Principal component',as.character(i),sep=" "))

# Plot loaddings
par(mar = c(6, 4, 2, 0.5))
image(x$lon,x$lat,x$eof[,,i]*100,col=colours,xlab = '', ylab = '',breaks=breaks)

if(cont) {contour(x$lon,x$lat,x$eof[,,i]*100,levels=round(breaks,1),labels=round(breaks,1),add=T,...)}

title(paste('EOF',as.character(i),' ',as.character(x$accountedvar[i]),'%'))

# check if lon is from 0 to 360 or -180 to 180 to use appropriate world map
if (min(x$lon)<0){
map('world',interior = FALSE,add = T, lwd=1) # Low resolution world map (lon -180 to 180)
}
else{
map('world2',interior = FALSE,add = T, lwd=1) # Low resolution world map (lon 0 360)
}
box()
map.axes()
#title(maintit)


# adding colorbar
par(mar = c(2.5, 2, 0, 0.5), mgp = c(1.5, 0.3, 0), las = 1)
image(c(1:n),1, t(t(c(1:n))), axes = F, col = colours,xlab = '', ylab = '')

box()
par(cex = 1.1)



## add the tick marks to the plot without plotting labels
axis(1, at = seq(1.5, length(breaks)-1.5), labels = F)

if(maximum>1){
## add the labels to the plot, without plotting tick marks
axis(1, at = seq(1.5, length(breaks)-1.5), tick = F,
     labels = round(breaks[2:(length(breaks)-1)],1))
}
else{
axis(1, at = seq(1.5, length(breaks)-1.5), tick = F,
     labels = round(breaks[2:(length(breaks)-1)],2))
}
## redifine font size
par(cex = 1.2)
## add the title
title(xlab = 'loadings')

}