library(ncdf)
library(maps)
library(fields)


# rm(list=ls()) #  Commented out, Remove pound sign to activate. ACTIVATION REMOVES ALL OBJECTS IN R WORKSPACE -- BE CAREFUL TO INSURE YOU WISH TO DO THIS.
 

#####Definitions to get file of reconstructed grid and to save graphics output file


Filepath<-("//Dinobast/paleo_fileshare/2000-reconstruction/MLOST/")  #REPLACE WITH APPROPRIATE FILEPATH FOR LOCAL CIRCUMSTANCE.

Filename<-("mlost-ann.nc")

Outputfilecase<-("MLOST_example_timeseries_plus_spatialplots.jpg")


#####Open output file of MLOST grid, as above


z<-open.ncdf(paste(Filepath,Filename,sep=""))

print (z)


time<-get.var.ncdf(z,"time")
dim(time)

#YearAD<-get.var.ncdf(z,"YearAD")
#dim(YearAD)


lat<-get.var.ncdf(z,"lat_bnds")
dim(lat)

lon<-get.var.ncdf(z,"lon_bnds")
dim(lon)

zlev<-get.var.ncdf(z,"zlev")
dim(zlev)


temp<-get.var.ncdf(z,"surf_temp_anom")  #[,,1:(dim(time))]
dim(temp)


lon2<-colMeans(lon)
lat2<-colMeans(lat)

lon<-lon2
lat<-lat2


dimnames(temp)<-list(lon,lat,time)


#####Checking overall mean of MLOST field, and setting up index and spatial mean values over Period of Record


mean(temp, na.rm=TRUE)  


meantemp_time<-colMeans(temp,dims=2, na.rm=TRUE)


temp_latlon<-array(dim=c(length(time),length(lat),length(lon)),dimnames=list(time,lat,lon))

 for (l in 1:length(time)) {
   for (j in 1:length(lat)) {
     for (k in 1: length(lon)) {
        temp_latlon[l,j,k]<-temp[k,j,l] }}}

dim(temp_latlon)    

mean(temp_latlon, na.rm=TRUE)  


meantemp_latlon<-colMeans(temp_latlon, dims=1, na.rm=TRUE)

dim(meantemp_latlon)


#####Plotting time series and mean spatial values -- transformations required to use R's "image.plot" function in correct orientation and addition of ocean boundaries


latinv<-(-lat)

meantemp_latlon2<-meantemp_latlon
temp_latlon2<-temp_latlon

  for (i in 1:length(lat)) {
    meantemp_latlon2[i,]<-meantemp_latlon[lat==latinv[i],]
    temp_latlon2[,i,]<-temp_latlon[,lat==latinv[i],] }       # these two lines reverse the T's in terms of lat, making them run from S to N, which conforms to the input requirements of R's "image.plot" function. 

LON<-lon
LAT<-latinv


timemeanmax<-max(meantemp_latlon2, na.rm=TRUE)  # These lines are used to check ranges of T values for full data set and averages over period of record
timemeanmax

timemeanmin<-min(meantemp_latlon2, na.rm=TRUE)
timemeanmin

overallmax<-max(temp_latlon2, na.rm=TRUE)
overallmax

overallmin<-min(temp_latlon2, na.rm=TRUE)
overallmin


    ## Plotting functions for different types of outputs saved in specified file.  Examples given for global average and spatial plots for entire period of record mean, bi-decadal and tri-decadal means, and specific year.


jpeg(filename=paste(Filepath,Outputfilecase,sep=""), width=1200, height=1200, pointsize=18, quality=100, bg="white")


par(mfrow=c(3,2), oma = c(2,2,2,2))


plot(time, meantemp_time, type = "l", cex.main = 2, cex.lab = 2, xlab = paste("Year"), ylab = paste("global mean temp  deg C"))
abline(h=0)

image.plot(LON,LAT,t(meantemp_latlon2),zlim=c(-2,2), cex.main = 2, cex.lab = 2, main=paste("1880-2010"), xlab = paste("lon"), ylab = paste ("lat"), legend.width = 1.5)
map("world",interior=FALSE,add=TRUE)


image.plot(LON,LAT,t(colMeans(temp_latlon2[time<=1920 & time>=1901,,], dims=1, na.rm=TRUE)),zlim=c(-2,2), cex.main = 2, cex.lab = 2, main=paste("1901-1920"), xlab = paste("lon"), ylab = paste ("lat"), legend.width = 1.5)
map("world",interior=FALSE,add=TRUE)


image.plot(LON,LAT,t(colMeans(temp_latlon2[time<=1950 & time>=1931,,], dims=1, na.rm=TRUE)),zlim=c(-2,2), cex.main = 2, cex.lab = 2, main=paste("1931-1950"), xlab = paste("lon"), ylab = paste ("lat"), legend.width = 1.5)
map("world",interior=FALSE,add=TRUE)


image.plot(LON,LAT,t(colMeans(temp_latlon2[time<=2000 & time>=1971,,], dims=1, na.rm=TRUE)),zlim=c(-2,2), cex.main = 2, cex.lab = 2, main=paste("1971-2000:  ref climatology"), xlab = paste("lon"), ylab = paste ("lat"), legend.width = 1.5)
map("world",interior=FALSE,add=TRUE)


image.plot(LON,LAT,t(temp_latlon2[time==1998,,]),zlim=c(-3,3), cex.main = 2, cex.lab = 2, main=paste("1998"), xlab = paste("lon"), ylab = paste ("lat"), legend.width = 1.5)
map("world",interior=FALSE,add=TRUE)


dev.off()



#####********************************* END ***************************************