;
trv=0           ; selects tree-ring-variable: 0=MXD 1=TRW 2=MXD-TRW
case trv of
  0: fnadd='mxd'
  1: fnadd='trw'
  2: fnadd='mxd-trw'
endcase
titadd=strupcase(fnadd)
;
; Get chronology locations
;
print,'Reading '+titadd+' data'
if trv eq 2 then begin
  restore,filename='../alltrw.idlsave'
  trw=mxd
  restore,filename='../allmxd.idlsave'
  mxd=mxd-trw
endif else begin
  restore,filename='../all'+fnadd+'.idlsave'
  ;  nchron,idno,idname,location,country,tree,yrstart,yrend,statlat,statlon,$
  ;  mxd,fraction,timey,nyr
endelse
;
; Now read in the 1961-90 monthly means of precip and temperature
;
print,'Reading precipitation baseline'
restore,'/cru/u2/f055/data/obs/grid/surface/precip_new_19611990.idlsave'
;  g,nmon,ltm5
pre6190=ltm5
print,'Reading temperature baseline'
restore,'/cru/u2/f055/data/obs/grid/surface/lat_new_19611990.idlsave'
;  g,nmon,ltm5
lat6190=ltm5
;
; Now read in the land precipitation dataset.
; Although there is some missing data in Mark New's precip anomalies, all
; MXD boxes have sufficient data, so do not have to use surrounding boxes.
;
print,'Reading precipitation data'
ncid=ncdf_open('/cru/u2/f055/data/obs/grid/surface/precip_new_19011995.mon.nc')
pmon=crunc_rddata(ncid,tst=[1901,0],tend=[1994,11],grid=g)
ncdf_close,ncid
;
; Now read in the land air temperature dataset.
;
print,'Reading temperature data'
ncid=ncdf_open('/cru/u2/f055/data/obs/grid/surface/lat_jones_18511998.mon.nc')
tmon=crunc_rddata(ncid,tst=[1881,0],tend=[1994,11],grid=gtemp)
ncdf_close,ncid
;
; The precipitation series will be padded with missing values at its
; beginning, to make it the same length as the temperature.  Find pad length,
; and also find section length of MXD data to extract.
;
nmon=12
npad=(gtemp.nt-g.nt)/nmon
ppad=replicate(!values.f_nan,nmon*npad)
klmxd=where((timey ge 1881) and (timey le 1994),nyrm)
timey=timey(klmxd)
;
; ALTHOUGH THE FULL PERIOD IS ALWAYS READ (TO ENSURE CONSISTENT EVALUATION
; OF WHICH TEMPERATURE BOXES DO NOT HAVE 300 OR MORE MONTHS OF NON-MISSING
; 1881-1964 DATA) THE FOLLOWING LINE ALLOWS SELECTION OF A SMALLER SUB-PERIOD
; FOR ANALYSIS.
;
klperiod=where((timey ge 1881) and (timey le 1994),nyrper)
print,'LENGTH OF ANALYSIS PERIOD (YEARS) =',nyrper,'      *************'
;
; Define arrays
;
thalf=10.    ; filter for smoothed series, lower the better cos of sample size
ncorr=21     ; 16 months (JJASOND-JFMAMJJAS) & 5 seasons (O-S O-M A-S M-A JJA)
nvar=2       ; 0=MXD,precip, 1=MXD,temperature
nyrp=g.nt/nmon
nyrt=gtemp.nt/nmon
if (nyrp+npad) ne nyrm then message,'Ooops!!!'
if nyrt ne nyrm then message,'Ooops!!!'
nyr=nyrm
z=!values.f_nan
allr=fltarr(nchron,ncorr,nvar)*z      ; correlations with P or T
lowr=fltarr(nchron,ncorr,nvar)*z      ; low freq correlation with P or T
moir=fltarr(nchron)*z                 ; corr with moisture deficit
moir2=fltarr(nchron)*z                ; corr Apr-Sep T vs moist deficit
allp=fltarr(nchron,ncorr,nvar)*z      ; partial correlations with P or T
lowp=fltarr(nchron,ncorr,nvar)*z    ; low freq partial correlations with P or T
varname=['Prec','Temp']
corrname=['J','J','A','S','O','N','D','J','F','M','A','M','J','J','A','S',$
  'Oct-Sep','Oct-Mar','Apr-Sep','May-Aug','Jun-Aug']
;
; Analyse each chronology in turn
;
runlen=0.
runnum=0.
for i = 0 , nchron-1 do begin
  print,i,format='($,I4)'
  ;
  ; Extract MXD timeseries
  ;
  mxd1=reform(mxd(*,i))
  mxd1=mxd1(klmxd)
  ;
  ; Locate grid box containing site and extract monthly timeseries
  ;
  x=statlon(i)
  y=statlat(i)
  statval=[1.]
  fdout=gridit(g.nx,g.ny,g.x,g.y,x,y,statval)
  loclist=where(finite(fdout),nloc)
  if nloc ne 1 then message,'Ooops!!!'
  locx=loclist(0) mod g.nx
  locy=fix(loclist(0)/g.nx)
  ;
  ; Simple single box mean for precipitation, but pad with missing
  ;
  pre12=[ppad,reform(pmon(locx,locy,*))]
  pre12=reform(pre12,nmon,nyr)
  ;
  ; For temperature, if less than 300 months of 1881-1964 data use an
  ; average of the nine surrounding boxes (variance adjusted etc.)
  ;
  tem12=reform(tmon(locx,locy,*),nmon,nyr)
  kl=where(gtemp.year le 1964)
  dummy=where(finite(tem12(*,kl)),ngot)
  if ngot lt 300 then begin
    print
    ; extract timeseries from 9 surrounding boxes (NB. zonally cyclic)
    locy=[locy-1,locy,locy+1]
    locx=[locx-1,locx,locx+1] mod g.nx
    neglist=where(locx lt 0,nneg)
    if nneg gt 0 then locx(neglist)=locx(neglist)+g.nx
    temboxes=tmon(locx,*,*)
    temboxes=temboxes(*,locy,*)
    temboxes=reform(temboxes,9,nmon,nyr)
    ; average them and adjust variance on a month by month basis
    ; I've altered mkeffective to adjust variance to that of a single
    ; input timeseries
    for imon = 0 , nmon-1 do begin
      print,imon,format='($,I4)'
      tsin=reform(temboxes(*,imon,*))
      tsout=var_adjust(tsin,/no_anomaly,/td_rbar,/mkrbar,/mkcorr)
      tem12(imon,*)=tsout(*)
    endfor
    print
  endif
  ;
  save,filename='example_chron'+string(i,format='(I3.3)')+'.idlsave',$
    mxd1,pre12,tem12,nmon,nyr,timey,i,location,statlat,statlon
  ;
endfor
;
end