;
; On a site-by-site basis, computes MXD timeseries from 1902-1976, and
; computes Apr-Sep temperature for same period, using surrounding boxes
; if necessary.  Normalises them over as common a period as possible, then
; takes 5-yr means of each (fairly generous allowance for
; missing data), then takes the difference.
; Results are then saved for briffa_sep98_decline2.pro to perform rotated PCA
; on, to obtain the 'decline' signal!
;
;----------------------------------------------------
;
; Get chronology locations and MXD series
;
print,'Reading MXD data'
restore,filename='allmxd.idlsave'
;  nchron,idno,idname,location,country,tree,yrstart,yrend,statlat,statlon,$
;  mxd,fraction,timey,nyr
;
; 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=[1984,11],grid=g)
ncdf_close,ncid
;
nmon=12
klmxd=where((timey ge 1881) and (timey le 1984),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 1984),nyrper)
print,'LENGTH OF ANALYSIS PERIOD (YEARS) =',nyrper,'      *************'
;
; Define arrays
;
avelen=5.    ; length of averaging blocks to use
nyrt=g.nt/nmon
if nyrt ne nyrm then message,'Ooops!!!'
nyr=nyrm
z=!values.f_nan
allmxd=fltarr(nchron,nyr)*z
alltemp=fltarr(nchron,nyr)*z
;
; Analyse each chronology in turn
;
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)
  ;
  ; 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(g.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
  ;
  ; Extract period for analysis
  ;
  mxd1=mxd1(klperiod)
  tem12=tem12(*,klperiod)
  ;
  ; Now compute Apr-Sep seasonal mean
  ;
  tem1=mkseason(tem12,3,8,datathresh=3)          ; weak threshold!
  ;
  ; Now need to normalise them.  Could normalise over an early common period,
  ; but some instrumental records are not long enough.  Still, those will
  ; drop out of the PCA anyway!
  ;
  mknormal,mxd1,timey,refperiod=[1901,1960]
  mknormal,tem1,timey,refperiod=[1901,1960]
  allmxd(i,*)=mxd1(*)
  alltemp(i,*)=tem1(*)
  ;
  ; Repeat for next site
  ;
endfor
print
;
; Now extract the 1902-1976 period
;
kl=where((timey ge 1902) and (timey le 1976),nyr2)
mxd5yr=allmxd(*,kl)
temp5yr=alltemp(*,kl)
;
; Now average into 5-yr periods, requiring at least 2 years data in each block
;
mxd5yr=reform(mxd5yr,nchron,5,nyr2/5)
totval=total(mxd5yr,2,/nan)
totnum=total(finite(mxd5yr),2)
ml=where(totnum lt 2,nmiss)
mxd5yr=totval/totnum
if nmiss gt 0 then mxd5yr(ml)=!values.f_nan
;
temp5yr=reform(temp5yr,nchron,5,nyr2/5)
totval=total(temp5yr,2,/nan)
totnum=total(finite(temp5yr),2)
ml=where(totnum lt 2,nmiss)
temp5yr=totval/totnum
if nmiss gt 0 then temp5yr(ml)=!values.f_nan
;
; Now save the data
;
nblock=nyr2/5
blockyr=findgen(nblock)*5.+1904.
save,filename='briffa_sep98_decline1.idlsave',$
  allmxd,alltemp,mxd5yr,temp5yr,nchron,nyr,nblock,statlat,statlon,$
  blockyr,timey
;
end