;
; Reads in the MXD data set, grids it onto a 5 by 5 grid, adjusts each grid
; box for sample size variations (requires rbar for each box, which is also
; computed here), and then normalises over the period 1881-1960 (although
; the period can be altered).  When gridding, multiple chronologies are
; weighted by their "fraction of cores" values before averaging.
; Only those with complete data over 1800-1964 are kept.  Those with complete
; data over 1800-1976 are used to infill missing values between 1965 and 1976,
; and those complete over 1697-1964 are used to infill missing values between
; 1697 and 1799.  Any that are not complete over 1697-1976 after infilling
; are removed, and only the 1697-1976 period is kept and saved.
; *** NOW REMOVES MXD SERIES WITH POOR CORRELATIONS WITH LOCAL SUMMER TEMP ***
; *** NOW IT ALSO DOES LOTS OF INFILLING USING VARIOUS PERIODS ***
;
; *** NOW ALSO WORKS FOR TRW! ***
;
;----------------------------------------------------------------------
;
dosave=0    ; 0=do not, 1=do save the resulting data set
;
beglist=[1400,1453,1540,1583,1612,1660,1697,1743,1777]
nbeg=n_elements(beglist)
endlist=[1994,1988,1980]
nend=n_elements(endlist)
;
trv=0           ; selects tree-ring-variable: 0=MXD 1=TRW
cutr=0.22      ; a 'best' series is based on sites with local r >= cutr (0.22)
case trv of
  0: begin
    fnadd='mxd'
    iseas=18        ; use local r with Apr-Sep temp to select chronologies
    end
  1: begin
    fnadd='trw'
    iseas=20        ; use local r with Jun-Aug temp to select chronologies
    end
endcase
titadd=strupcase(fnadd)
;
; First read in the MXD data
;
print,'Reading '+titadd+' data'
restore,filename='../all'+fnadd+'.idlsave'
;  nchron,idno,idname,location,country,tree,yrstart,yrend,statlat,statlon,$
;  mxd,fraction,timey,nyr
;
; Next read in the correlations between MXD and local climate
;
print,'Reading '+titadd+' local correlations'
restore,filename='../datastore/'+fnadd+'_moncorr.idlsave'
;  allr,ncorr,nvar,nchron,varname,corrname,statlat,statlon,allp,moir,moir2
;
; Remove any that are almost entirely missing, or that
; aren't well correlated with their local summer temperature
;
nkeep=0
for i = 0 , nchron-1 do begin
  if allr(i,iseas,1) ge cutr then begin
    dummy=where(finite(mxd(*,i)),ngot)
    if ngot ge 20 then begin
      if nkeep eq 0 then kl=i else kl=[kl,i]
      nkeep=nkeep+1
    endif
  endif
endfor
print,'Had',nchron,' chronologies, but',nchron-nkeep,' are empty or poor'
nstat=nkeep
density=mxd(*,kl)
weight=fraction(*,kl)
statlon=statlon(kl)
statlat=statlat(kl)
statcty=country(kl)
;
;----------------------------------------------------------------------
;
; Get grid info (IPCC instrumental dataset grid), NH only
;
g=def_grid(/ipcc)
ky=where(g.y ge 0.,newny)
g = { nx: g.nx, ny: newny, x: g.x, y: g.y(ky) }
;
;----------------------------------------------------------------------
;
; Create a box by box list of which chronologies fall in each box
;
print,'Locating the chronologies into grid boxes'
chronlists=ingrid(g.x,g.y,statlon,statlat,nstat=chrondens)
maxlist=max(chrondens)
keepbox=where(chrondens gt 0.,maxbox)
keepi=keepbox mod g.nx
keepj=fix(keepbox/g.nx)
chronlists=reform(chronlists,maxlist,g.nx*g.ny)
chrondens=reform(chrondens,g.nx*g.ny)
;
;----------------------------------------------------------------------
;
; Define array to contain box by box density timeseries
;
boxdens=fltarr(maxbox,nyr)
cty=strarr(maxbox)
;
;----------------------------------------------------------------------
;
; For each box, extract chronologies, compute rbar if more than one,
; compute weighted mean chronology, adjust variance and normalise
;
print,'Computing variance adjusted, normalised, box-mean timeseries'
;print,replicate('-',maxbox),format='(80A1)'
refper=[1881,1960]
for ibox = 0 , maxbox-1 do begin
  ;print,'.',format='($,A1)'
  nlist=chrondens(keepbox(ibox))
  onelist=chronlists(0:nlist-1,keepbox(ibox))
  mxd=density(*,onelist)
  wt=weight(*,onelist)
  cty(ibox)=statcty(onelist(0))
  if nlist gt 1 then begin
    ; compute weighted mean of chronologies in the box
    totwtmxd=total(mxd*wt,2,/nan)
    totwt=total(float(finite(mxd))*wt,2)
    wtmxd=totwtmxd/float(totwt)
    ; compute rbar between chronologies
    inx=transpose(mxd)
    onerbar=mkrbar(inx)
    print,nlist,onerbar
    ; adjust variance for sample size
    nsample=total(mxd*0.+1.,2,/nan)
    adjmxd=mkeffective(wtmxd,nsample,rbar=onerbar)
  endif else begin
    ; if only one chronology in box, then no need to do anything to it
    print,nlist
    adjmxd=reform(mxd)
  endelse
  ; now normalise
  mknormal,adjmxd,timey,refperiod=refper
  boxdens(ibox,*)=adjmxd
endfor
print
;
; Now store the full gridded time series into a northern hemisphere array
;
mxdfd=fltarr(g.nx,g.ny,nyr)*!values.f_nan
for i = 0 , maxbox-1 do mxdfd(keepi(i),keepj(i),*)=boxdens(i,*)
mxdyear=timey
mxdnyr=nyr
mxdfd2=mxdfd          ; mxdfd2 will also contain the infilled data
;
;----------------------------------------------------------------------
;
; Remove those that are too short (assume that they are continuous, so just
; look for the first non-missing value)
;
maxstart=1810
print,'Removing boxes that dont precede',maxstart
yrstart=fltarr(maxbox)
yrend=fltarr(maxbox)
for ibox = 0 , maxbox-1 do begin
  adjmxd=boxdens(ibox,*)
  keeplist=where(finite(adjmxd),nkeep)
  yrstart(ibox)=timey(keeplist(0))
  yrend(ibox)=timey(keeplist(nkeep-1))
endfor
keeplist=where(yrstart le maxstart,nkeep)
print,'Reduce boxes from',maxbox,' to',nkeep
maxbox=nkeep
keepbox=keepbox(keeplist)
keepi=keepi(keeplist)
keepj=keepj(keeplist)
boxdens=boxdens(keeplist,*)
cty=cty(keeplist)
yrstart=yrstart(keeplist)
yrend=yrend(keeplist)
;
; Make a copy containing un-infilled data, and one to contain the final
; multi-infilled data
;
infilldens=boxdens
nofilldens=boxdens
;
;----------------------------------------------------------------------
;
; Locate those that are complete for 1800 to 1985 and use them to infill
; those that stop between 1965 and 1985
; *** NOW TRY IT FOR 1697 to 1976 (280 years) ***
;
for iend = 0 , nend-1 do begin
  lastyr=endlist(iend)
;  print,'Infilling to',lastyr
  baselist=where(yrend ge lastyr,nbase)
  if nbase eq 0 then message,'Cannot infill'
  filllist=where(yrend lt lastyr,nfill)
  print,'Out of',maxbox,' boxes,',nbase,' reach',lastyr,'. Try to fill in',nfill
  ncant1=0
  if nfill gt 0 then begin
    boxdens=nofilldens
    infill,boxdens,g.x(keepi),g.y(keepj),timey,boxname=cty,$
      baselist=baselist,filllist=filllist,$
      overlap=[1810,1970],fillperiod=[1971,lastyr],$
      search=[800.,1600.,2400.,3000.],$
      cantlist=cantlist1
    ncant1=n_elements(cantlist1)
    print,'Cannot fill',ncant1,', so can fill',nfill-ncant1
    newlist=where(finite(boxdens),nnew)
    if nnew gt 0 then infilldens(newlist)=boxdens(newlist)
  endif
  print
endfor
;
;----------------------------------------------------------------------
;
; Locate those that are complete for 1700 to 1964 and use them to infill
; those that start between 1700 and 1800
;
for ibeg = 0 , nbeg-1 do begin
  begyr=beglist(ibeg)
;  print,'Infilling to',begyr
  baselist=where(yrstart le begyr,nbase)
  if nbase eq 0 then message,'Cannot infill'
  filllist=where(yrstart gt begyr,nfill)
  print,'Out of',maxbox,' boxes,',nbase,' reach',begyr,'. Try to fill in',nfill
  if nfill gt 0 then begin
    boxdens=nofilldens
    infill,boxdens,g.x(keepi),g.y(keepj),timey,boxname=cty,$
      baselist=baselist,filllist=filllist,$
      overlap=[1810,1970],fillperiod=[begyr,1809],$
      search=[800.,1600.,2400.,3000.],$
      cantlist=cantlist2
    ncant2=n_elements(cantlist2)
    print,'Cannot fill',ncant2,', so can fill',nfill-ncant2
    newlist=where(finite(boxdens),nnew)
    if nnew gt 0 then infilldens(newlist)=boxdens(newlist)
  endif
  print
endfor
;
; Store the infilled ones in the full array
;
boxdens=infilldens
for i = 0 , maxbox-1 do mxdfd2(keepi(i),keepj(i),*)=boxdens(i,*)
;
;----------------------------------------------------------------------
;
; Get rid of any that couldn't be infilled, and keep only 1700 to 1985 only.
; Any that contain any missing data during that period are removed.
;
if ncant1 gt 0 then boxdens(cantlist1,*)=!values.f_nan
if ncant2 gt 0 then boxdens(cantlist2,*)=!values.f_nan
;
keeplist=where((timey ge begyr) and (timey le lastyr),nkeep)
timey=timey(keeplist)
boxdens=boxdens(*,keeplist)
ntime=nkeep
;
totdens=total(boxdens,2)
keeplist=where(finite(totdens),nkeep)
print,'Ended up with',nkeep,' boxes'
keepi=keepi(keeplist)
keepj=keepj(keeplist)
boxdens=boxdens(keeplist,*)
oldn=maxbox
maxbox=nkeep
;
; Also work out where the un-infilled (complete) boxes are in the new
; reduced dataset
;
ilist=fltarr(oldn)        ; define 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...
ilist(baselist)=1.        ; mark   0 0 1 1 0 1 0 0 1 1 1 0 0 1 1 0 1 0 0 ...
ilist=ilist(keeplist)     ; remove 0 0 1 1 0 0 1 1 1 0 1 1 0 0 ...
baselist=where(ilist eq 1,nbase) ; locate remaining 1's
;
;----------------------------------------------------------------------
;
if dosave ne 0 then begin
  save,filename=fnadd+'_infill.idlsave',ntime,maxbox,boxdens,g,keepi,keepj,timey,$
    baselist,nbase,mxdfd,mxdyear,mxdnyr,mxdfd2,beglist,endlist
endif
;
end