####
## make GISS TRP averages from zonal averages archived at GISS
###
### giss.trp - a list of 4 time series 
###
### see comparison to averages from gridded data in notes

	#eff_month=paste( substr(Sys.Date(),1,4),"_",numbertochar( as.numeric(substr(Sys.Date(),6,7) )-1),sep="")


source("http://www.climateaudit.info/scripts/gridcell/read.binary.txt")
     

######################################################################
## ZON VERSIONS
######################################################################
numbertochar =function(x) if (x<10) paste ("0",x,sep="") else paste(x)


#this constructs zonal mean from gridded

make.giss.trp=function(dset="loti_1200"){ 
  if(dset=="loti_250")  loc="ftp://ftp.cdc.noaa.gov/Datasets/gistemp/combined/250km/air.2x2.250.mon.anom.comb.nc"
  if(dset=="loti_1200") loc="ftp://ftp.cdc.noaa.gov/Datasets/gistemp/combined/1200km/air.2x2.1200.mon.anom.comb.nc"
  if(dset=="land_250")  loc="ftp://ftp.cdc.noaa.gov/Datasets/gistemp/landonly/250km/air.2x2.250.mon.anom.land.nc"
  if(dset=="land_1200") loc="ftp://ftp.cdc.noaa.gov/Datasets/gistemp/landonly/1200km/air.2x2.1200.mon.anom.land.nc"	
	download.file(loc,"temp.dat",mode="wb") #190 MB
	nc<-open.ncdf("temp.dat")

    #only one dim here
	v1      <- nc$var[[1]]
	varsize <- v1$varsize #180 90 1595
	ndims   <- v1$ndims
	nt      <- varsize[ndims];nt   # Remember timelike dim is always the LAST dimension!
	giss.trp=ts(rep(NA,varsize[3]),start=c(1880,1),freq=12)	 #to Feb 2008
	temp=(v1$dim[[2]]$vals<=25)&(v1$dim[[2]]$vals> -25)

    for( i in 1:nt ) {
        # Initialize start and count to read one timestep of the variable.
        start <- rep(1,ndims)   # begin with start=(1,1,1,...,1)
        start[ndims] <- i       # change to start=(1,1,1,...,i) to read timestep i
        count <- varsize        # begin w/count=(nx,ny,nz,...,nt), reads entire var
        count[ndims] <- 1       # change to count=(nx,ny,nz,...,1) to read 1 tstep
        data3 <- get.var.ncdf( nc, v1, start=start, count=count ) #360 180
	giss.trp[i]= mean(data3[,temp],na.rm=TRUE)
        }
   	close.ncdf(nc)
	rm(nc)
	rm(data3)
	return(giss.trp)
} # 

	#giss[["loti_1200"]]=make.giss.trp(dset="loti_250")
get_zonal=function() {
         #was "ftp://data.giss.nasa.gov/pub/gistemp/GISS_Obs_analysis/ARCHIVE/03/ZON/LOTI.zon.web"
	 url=file.path("ftp://data.giss.nasa.gov/pub/gistemp/GISS_Obs_analysis/zonal_means","LOTI.zon.web")
	   #ftp://data.giss.nasa.gov/pub/gistemp/GISS_Obs_analysis/zonal_means/  
	   #now  ftp://data.giss.nasa.gov/pub/gistemp/GISS_Obs_analysis/zonal_means/LOTI.zon.web
	loti=read.zon.web(url);  tsp(loti)
		#[1] 1880.00 2009.25   12.00
	  #dim loti 1573 x 46: the columns are -90 to +90 by 2
	lat=as.numeric(dimnames(loti)[[2]]); temp= abs(lat)<=20;sum(temp)
	giss.trp[["loti"]]=ts(apply(loti[,temp],1,mean,na.rm=TRUE),start=1880,freq=12)

  #sst
	url=file.path("ftp://data.giss.nasa.gov/pub/gistemp/GISS_Obs_analysis/zonal_means","SST.zon.web")
  	x=read.zon.web(url)
	tsp(x)
		#[1] 1880.00 2009.25   12.00
	lat=as.numeric(dimnames(x)[[2]]); temp= abs(lat)<=20
	giss.trp[["sst"]]=ts(apply(x[,temp],1,mean,na.rm=TRUE),start=1880,freq=12)


   #land1200
	url=file.path("ftp://data.giss.nasa.gov/pub/gistemp/GISS_Obs_analysis/zonal_means","Ts.zon.web")
  	x=read.zon.web(url) ;	tsp(x)
		#[1] 1880.00 2009.25   12.00
	lat=as.numeric(dimnames(x)[[2]]); temp= abs(lat)<=20
	giss.trp[["land1200"]]=ts(apply(x[,temp],1,mean,na.rm=TRUE),start=1880,freq=12)

	
   #land250
	url=file.path("ftp://data.giss.nasa.gov/pub/gistemp/GISS_Obs_analysis/zonal_means","Ts.zon250.web")
  	x=read.zon.web(url) ;	tsp(x)
		#[1] 1880.00 2009.25   12.00
	lat=as.numeric(dimnames(x)[[2]]); temp= abs(lat)<=20
	giss.trp[["land250"]]=ts(apply(x[,temp],1,mean,na.rm=TRUE),start=1880,freq=12)

return(giss.trp)
}


 #CHECK ftp://data.giss.nasa.gov/pub/gistemp/GISS_Obs_analysis/ARCHIVE	
	
	

#gives giss.trp: REQUIRES VERY LARGE DOWNLOOAD
#has to be done step by step since matrix large

#Users interested in the entire gridded temperature anomaly data may download 
#the three basic binary files from our ftp site. Also available there are various FORTRAN programs 
#and instructions to create (time series of) regular gridded anomaly maps from these files. 
#This should make the maintenance of the files mentioned below unnecessary. 



###################################################
### three gridded versions: SST and two land extrapolated
###
### ftp://data.giss.nasa.gov/pub/gistemp/SBBX.HadR2 SST
### #ftp://data.giss.nasa.gov/pub/gistemp/SBBX1880.Ts.GHCN.CL.PA.250  18577
### #ftp://data.giss.nasa.gov/pub/gistemp/SBBX1880.Ts.GHCN.CL.PA.1200  42809
###
###
### ftp://data.giss.nasa.gov/pub/gistemp/GISS_Obs_analysis/lplots/LOW.Ts.GHCN.CL.PA.lpl
## ftp://data.giss.nasa.gov/pub/gistemp/GISS_Obs_analysis/GISTEMP_sources/
### zonal 
###
### read.Zon.Ts - function to read zonal binary files at GISS 
###
### lplots/LOW.Ts.GHCN.CL.PA.lpl 
###    LOW LATITUDES: 23.6 S TO  23.6 N CR 1200KM     GHCN V2 TemperatureANOM (C)                                        
###
### ftp://data.giss.nasa.gov/pub/gistemp/GISS_Obs_analysis/lplots/LOW.Tsho2.GHCN.CL.PA.lpl


#http://data.giss.nasa.gov/work/modelEt/time_series/work/JTmap.txt
#ftp://data.giss.nasa.gov/pub/gistemp/download/README.txt
#ftp://data.giss.nasa.gov/pub/gistemp/GISS_Obs_analysis/ARCHIVE/ds/STEP6/results/
	#Ethese are lpl plots
#http://www.climateaudit.org/?p=1228#comment-92539