close,/all
	print,'data? 0-NCEP, 1-ERA40, 2-20C'
	read,idata

	if(idata eq 0) then begin
	  cdir='/Volumes/datasets/ncep.reanalysis.derived/pressure/'
	  cvar0=['air','hgt','uwnd','vwnd','omega']
	  yr0=1948.0
	  nyr=64
	  mul=fltarr(5,6)+1.0
	  mul(0,*)=[-1,-1,-1,-1,-1,-1]
	endif
	if(idata eq 1) then begin
	  cdir='/Volumes/datasets/ecmwf.reanalysis-40.derived/pressure/'
	  cvar0=['air','geopot','uwnd','vwnd','omega']
	  yr0=1958.0
	  nyr=43
	  mul=fltarr(5,6)+1.0
	endif
	if(idata eq 2) then begin
	  cdir='/Volumes/datasets2/20C/pressure/'
	  cvar0=['air','hgt','uwnd','vwnd','omega']
	  yr0=1870.0
	  nyr=140
	  mul=fltarr(5,6)+1.0
	endif
	print,cvar0
	read,ivar
	cvar=cvar0(ivar)

	print,'remove seasonal? 0-n, 1-y'
	read,iseason

	file=cdir+cvar+".mon.mean.nc"
	file_id=ncdf_open(file)
	file_info=ncdf_inquire(file_id)
     	nvars=file_info.nvars      
     	print, "No. of variables", nvars
     	for var_id=0,nvars-1 do begin
          varinfo=ncdf_varinq(file_id, var_id)
          print, var_id, "  ", varinfo.name        
     	endfor

	ncdf_varget,file_id, "lon", lon
	ncdf_varget,file_id, "lat", lat
	ncdf_varget,file_id, "time", time
	ncdf_varget,file_id, "level", p0
	ncdf_varget,file_id, cvar, d0
	ncdf_attget, file_id, cvar, "missing_value", missing_value
	ncdf_attget, file_id, cvar, "scale_factor", scale_factor
	ncdf_attget, file_id, cvar, "add_offset", add_offset
	print,"missing_value", missing_value
	print, "scale_factor", scale_factor
	print, "add_offset", add_offset

        nlon=n_elements(lon)
	d=total(d0,1)/nlon
	d=d*scale_factor+add_offset

	ncdf_close,file_id

        nlat=n_elements(lat)
        nt=n_elements(time)
	nt=nyr*12
	t=yr0+findgen(nt)/12

	d=d(*,0:11,0:nt-1)
	p=p0(0:11)
	np=n_elements(p)

	np2=10
	p2=1000-findgen(np2)*100.0
	d2=fltarr(nlat,np2,nt)
	for j=0,nlat-1 do begin
	  for l=0,nt-1 do begin
	    temp=d(j,*,l)
	    d2(j,*,l)=interpol(temp,p,p2)
	  endfor
	endfor
	np=np2
	p=p2
	d=d2

	nlat2=18
	lat2=-85.0+findgen(nlat2)*10.0
	d2=fltarr(nlat2,np,nt)
	lat3=-88.75+findgen(72)*2.5
	for k=0,np-1 do begin
	  for l=0,nt-1 do begin
	    temp=reform(d(*,k,l))
	    temp2=interpol(temp,lat,lat3)
	    for j=0,nlat2-1 do begin
	      j1=j*4
	      j2=j1+3
	      d2(j,k,l)=total(temp2(j1:j2))/4.0
	    endfor
	  endfor
	endfor
	nlat=nlat2
	lat=lat2
	d=d2

	dc=fltarr(12)
	for j=0,nlat-1 do begin
	  for k=0,np-1 do begin
	    temp=reform(d(j,k,*))
	    if(iseason eq 0) then goto,jump9
	    dc(*)=0.0
	    for iyr=0,nyr-1 do begin
	      l1=iyr*12
	      l2=l1+11
	      dc=dc+temp(l1:l2)
	    endfor
	    dc=dc/nyr
	    for iyr=0,nyr-1 do begin
	      l1=iyr*12
	      l2=l1+11
	      temp(l1:l2)=temp(l1:l2)-dc
	    endfor
	    temp=smooth(temp,13,/edge_truncate)
jump9:
	    temp=temp-mean(temp)
	    temp=temp/stdev(temp)
	    d(j,k,*)=temp
	  endfor
	endfor

        calc_eof,d,u,var,c0

        !p.multi=[0,3,6] 
	loadct,39

        fs_base,nt,nf,f,cosu,sinu
        e=fltarr(nlat,np)
        for imode=0,5 do begin
          e(*,*)=u(imode,*,*)*mul(ivar,imode)

	  em=max(abs(e))
	  nlev=6
	  lev=(findgen(nlev*2+1)-nlev)/nlev*em
	  lin=[intarr(nlev)+1,intarr(nlev+1)]
	  thi=[intarr(nlev)+1,4,intarr(nlev)+1]

          contour,e,lat,p,xrange=[-90,90],xticks=6,yrange=[1000,100],xstyle=1,ystyle=1, $
                  charsize=2.0, $
                  xtitle="Latitude (deg)",ytitle="P (mb)", $
		  levels=lev,c_linestyle=lin,c_thick=thi, $
                  title="EOF"+string(imode+1)+" Var="+string(var(imode))+"%"
;          surface,e,lon,lat,shades=bytscl(e,top=255)

	  c=reform(c0(imode,*))*mul(ivar,imode)
	  c=c-mean(c)
          fs,nt,cosu,sinu,c,a,b,psd
	  psd2=smooth(psd,3,/edge_truncate)

	  plot,t,c,thick=4.0,xrange=[t(0),t(0)+nyr],xstyle=1,charsize=2.0
	  oplot,t,t*0.0
	  period=[50,20,10,3,1,1.0/12]
	  xti=n_elements(period)-1
	  plot,f(1:*),psd(1:*),/xlog,xticks=xti,xtickv=1.0/12/period,charsize=2.0
	  oplot,f,psd2,thick=4
        endfor

        end