#reviewed March 11, 2010 
	#new data set ftp://srvx6.img.univie.ac.at/pub/rich_gridded_2009.nc
 # see  http://www.univie.ac.at/theoret-met/research/raobcore/ Feb 2011 

source("http://www.climateaudit.info/scripts/utilities/getx_var_ncdf.txt") 

get.rich.trp=function(loc= "ftp://srvx7.img.univie.ac.at/pub/rich_gridded_2010.nc", 
            method="trp") {
	#loc="ftp://srvx6.img.univie.ac.at/pub/rich_gridded.nc"  #formerly 2008 
	# 2011/01/17 downloaded Feb 2011
	#loc="ftp://srvx6.img.univie.ac.at/pub/rich_gridded_2009.nc"  #formerly 2008 
		### loc="ftp://srvx6.img.univie.ac.at/pub/rich_gridded_clim_1979-1998.nc" : climatology 
	download.file(loc,"temp.nc",mode="wb")
	v<-open.ncdf("temp.nc")
	  length(v$dim)
	  names(v$dim) # 
	  	#[1] "lon"      "lat"      "pressure" "time"  
	  sapply(v$dim,function(A) A$len)
	  	# lon      lat pressure     time 
      	  	# 36       18       12      612 
	  v$dim$lon$vals # 10 degree Dateline to Dateline
	  v$dim$lat$vals # 10 degree N to S
	  v$dim$pressure$vals # 850 700 500 400 300 250 200 150 100  70  50  30
	  K=length(v$dim$pressure$vals) #12
	  v$dim$time$vals # "days since 1958-01-01 00:00:00"   	#monthly
	  year=c(time(ts( v$dim$time$vals,start=1958,freq=12) ))
	    # 1958.000 2008.917   12.000
	  names(v$var)  	  #anomalies
	  sapply(v$var$anomalies$dim,function(A) A$len)
		#  36  18  12 612
	instr <- getx.var.ncdf( v, v$var[[1]]) # 1850 2006
	   dim(instr) #36  18  12 588
	temp= abs(v$dim$lat$vals)<20;sum(temp) #4 
	instr=instr[,temp,,];dim(instr) # 36   4  12 612 
	instr=aperm(instr,c(4,1,2,3)); dim(instr) # 588  36   4  12
	dim0=dim(instr)
	instr=array(instr,dim=c(dim0[1],144,12) );dim(instr) # 624 144  12
	 mask=make.ten_degree.mask() # this gives tropical mask with grid arranged N to S and then dateline to dateline
	 mask= mask$mask 
	 X=rep( list(NA),3); names(X)=c("land","ocean","total")
	for (k in 1:3) { X[[k]]= ts( array (NA,dim=c(dim0[1],12) ),start= 1958,freq=12);  #tsp(rich[[k]]) # 1958.000 2006.917 
	  	 dimnames(X[[k]])[[2]]= v$dim$pressure$vals }# 850 700 500 400 300 250 200 150 100  70  50  30
	 for(i in 1:12)  X[["total"]][,i]= round(apply(instr[,,i],1,mean,na.rm=T) ,4)
	 for(i in 1:12)  X[["land"]][,i]= round(apply(instr[,,i],1,function(x) weighted.mean(x,1-mask,na.rm=T) ) ,4)
	 for(i in 1:12)  X[["ocean"]][,i]= round(apply(instr[,,i],1,function(x) weighted.mean(x,mask,na.rm=T) ) ,4)
	return(X)
}


make.ten_degree.mask=function() {
	loc="ftp://ftp.cdc.noaa.gov/Datasets/noaamergedtemp/lsmask.nc"	
	download.file(loc,"temp.dat",mode="wb")
	download.file(loc,"d:/climate/data/gridcell/noaa/lsmask.nc",mode="wb")

	v<-open.ncdf("temp.dat")
	names(v$var) # [1] lsmask
	mask <- get.var.ncdf( v, v$var[[1]]) # 1880 2008
	dim(mask) # 72 36
	names(v$dim) #[1] "time"  "lat"   "lon"   "nbnds"
	v$dim$lat$vals # 87.5N to 87.5S
	v$dim$lon$vals # 177.5W to 177.5E
	  #range(mask)# 0 to 1

	grid= data.frame(lat=as.numeric(as.character(gl(36,72,labels=v$dim$lat$vals))),lon=rep( v$dim$lon$vals,36),mask=c(mask))
	temp= abs(grid$lat)<20;sum(temp)
	grid=grid[temp,]
 	grid$latclass= 5+ 10*floor(grid$lat/10);  length(unique(grid$latclass))
	grid$longclass= 5+ 10*floor(grid$lon/10);  length(unique(grid$longclass))
	grid$class= 100*grid$latclass+grid$longclass; length(unique(grid$class)) #144
	Grid=data.frame(lat=as.numeric( as.character(gl(4,36,labels=seq(15,-15,-10)))),lon=rep( seq(-175,175,10),4 ) )
	Grid$class= 100*Grid$lat+Grid$lon
	x=tapply(grid$mask,grid$class,mean) 
	test=match(Grid$class,as.numeric(names(x))); sum(!is.na(test)) #144
	Grid$mask=x[test]
	Grid=Grid[,c(1,2,4)]
	return(Grid)
}

	#mask=make.ten_degree.mask()