#PART 1 of this script takes Jones data from CRU, changes -9999 to NA and saves R-tables v - with dimension 1769 x 2592 (2592=72*36)
#some of these steps take a few minutes since the dataset is large
#the downloading from CRU and unzipping needs to be done manually, since I don't know how to unzip within R although it looks possible

#PART 2 takes the saved table, truncates it to 1902-1995 and Mann's 1082 cells, scales to SD=1 and lat-weights

#PART 0 - MANUAL PART
#make directories "c:/climate", "c:/climate/data", "c:/climate/data/jones", "c:/climate/data/mann"
#Download Jones data from http://www.cru.uea.ac.uk/ftpdata/hadcrut2.dat.gz
#Unzip to "c:/climate/data/jones/hadcrut2.dat"
#open hadcrut2.dat in Notepad and see what year and last month are. As of July 14, 2003, it was (2003,5).
#the above is inefficient way of checking on final date, but it's only once	

#PART 1 - PUTTING JONES DATA INTO R- TABLE
loc<-"c:/climate/data/jones/hadcrut2.dat";
v<-scan(file=loc,n=36*(2006-1855)*12*75);#reads 4,776,300 items (July 2003)
n1<-length(v)/75;
v<-array(v,dim=c(75,n1));
v<-t(v);

##dim v (July 2003) is 63684x75 ## 63684=36 lat groups * 1769 months from 1856 Jan to 2003 May
##column 1 is year; column 2 is latitude group 1,,,36 north to south; column 3 is field in which 1=temperature; columns 4,..,75 are temperature anomalies for longitude groups 1,...72 from west to east (178.5W to 178.5E
##NA values indicated by -9999

#REPLACE -9999 WITH NA
marker<--9000;
temp<-(v<marker);
v[temp]<-NA;
rm(temp);

#MAKE ARRAY WITH COL1-MONTH, COL2-LONG GROUP, COL 3 - LAT GROUP
n<-n1/36;
v<-array(v,dim=c(36,n,75)) # this is a nice example of R: now indexed as (lat gp, month, long gp) and is very fast
#v now has dimension 36x 1769 (July 2003) x75
v<-aperm(v,perm=c(2,3,1))#dim = 1769x75 x 36
v<-v[,4:75,];##this truncates to temperature measurements only
save(v,file="c:/climate/data/jones/hadcruv.tab")##this saves 3-indexed array in new table dim=1769x72x36
q()

#this has been checked against Jones table and is a correct sorting#v[1769,,8] reproduces results from data-read
##v<-array(v,dim=c(n,36*72))##this is useful to have 2-D time series. you may prefer to save in this form instead of 3-D.
##this puts the time series for grid box (i,j) into column 36*(i-1)+(j-1)

#PART 2 - TRUNCATING, SCALING AND LAT-WEIGHTING JONES DATA PER MBH98 RECIPE 
#this requires prior downloading of Mann data from his FTP site. See "read.mann"

base.out<-"c:/climate/data/mann"

#LOAD MANN GRIDPOINTS
load("c:/climate/data/mann/gridpoints.tab"); #loads grid2 with grid2[,5] being 1082 Jones grid-box IDs
tempmann<-sort(grid2[,5])##loads Mann grid-boxes as sorted Jones ID numbers

#LOAD RAW JONES DATA AND MAKE 1902-95 SUBSET FOR MANN
load("c:/climate/data/jones/hadcruv.tab") ###month 553 is 1902 Jan;
v<-array(v,dim=c(dim(v)[1],36*72));####this puts the time series for grid box (i,j) into column 36*(i-1)+(j-1) and takes about 5 seconds

#truncate to Mann period of 1902-1995
n1<-1995-1901; ##this is period described in MBH98,786
period<-(1901-1855)*12+1:(n1*12) ##this selects month-index for 1902 Jan to 1995 Dec in JOnes v-matrix #length 1128
v<-v[period,]## this truncates v to 1902-1995 dim=1128x2592
dim(v); ##1128x2592

#truncate to Mann's 1082 cells
v<-v[,tempmann]
dim(v) #  1128 x 1082

#now scale to mean 0, sd 1 on each cell
mann.scale<-apply(v,2,sd,na.rm=TRUE)
mann.mean<-apply(v,2,mean,na.rm=TRUE)
save(mann.scale,file="c:/climate/data/mann/mann.scale.tab")
save(mann.mean,file="c:/climate/data/mann/mann.mean.tab")
v<-scale(v)

#CALCULATE COS(LATITUDE) TO CREATE LAT-BAND WEIGHTS FOR MANN EOF
mann.mu<-grid2[,4]
mann.mu<-pi*mann.mu/180
mann.mu<-cos(mann.mu) 

#CREATE LAT-WEIGHTED SCALED DATASET
mu1<-array(rep(mann.mu,1128),dim=c(1082,1128))
mu1<-t(mu1)
v<-v*mu1
rm(mu1)
save(v,file=file.path(base.out,"v.scaled.latweighted.tab")) #this saves 1128 x 1082 matrix v
q()