;
; Creates regional timeseries of MXD series, variance adjusted with a
; time-dependent rbar.
;
; The MXD sites are limited by requiring a minimum correlation against
; their local temperature.
;
;-----------------------------------------------------------------
;
trv=0           ; selects tree-ring-variable: 0=MXD 1=TRW
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)
;
; Specify options
;
fixedyr=1950              ; year of fixed grid
rbarper=[1881,1960]       ; period for computing correlation matrix
cutr=0.22       ; a 'best' series is based on sites with local r >= cutr
;
; 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 regional breakdown
;
print,'Reading MXD regions'
restore,filename='reg_mxdlists.idlsave'
;  ntree,treelist,nreg,regname
;
; Next read in the correlations between MXD and local climate
;
print,'Reading '+titadd+' local correlations'
restore,filename=fnadd+'_moncorr.idlsave'
;  allr,ncorr,nvar,nchron,varname,corrname,statlat,statlon,allp,moir,moir2
;
; Define arrays for storing all regional results
;
z=!values.f_nan
maxntree=max(ntree)
bestmxd=fltarr(nyr,nreg)*z ; MXD series
sampmxd=intarr(nyr,nreg)   ; MXD series
cutsites=fltarr(nreg)      ; no. of sites with data to fixed yr & with local r
cutx=fltarr(maxntree,nreg)*z
cuty=fltarr(maxntree,nreg)*z
;
; Now analyse each region in turn
;
for ireg = 0 , nreg-1 do begin
  print,regname(ireg)
  ;
  ; First select sites in region
  ;
  nx=ntree(ireg)
  print,'Maximum number of sites',nx
  tl=treelist(0:nx-1,ireg)
  ;
  regmxd=mxd(*,tl)
  regfra=fraction(*,tl)
  regcor=reform(allr(tl,iseas,1)) ; iseas=18/20=Apr-Sep/Jun-Aug 1=temperature
  xxx=statlon(tl)
  yyy=statlat(tl)
  ;
  ; Select those sites with data back to required year
  ;
  iyr=where(timey eq fixedyr)
  kl=where(finite(regmxd(iyr,*)),nx)
  print,'Fixed grid number of sites',nx
  if nx gt 0 then begin
    regmxd=regmxd(*,kl)
    regfra=regfra(*,kl)
    regcor=regcor(kl)
    xxx=xxx(kl)
    yyy=yyy(kl)
    ;
    ; Transpose to give space then time
    ;
    regmxd=transpose(regmxd)
    regfra=transpose(regfra)
    ;
    ; Compute correlation matrix of all retained sites
    ;
    fullrbar=mkrbar(regmxd,grandr=corrmat)
    ;
    ; Identify decreasing order of local skill
    ;
    slist=reverse(sort(regcor))
    ;
    ; Identify those with significant local correlations (approx.)
    ;
    dummy=where(regcor ge cutr,n2cut)
    cutsites(ireg)=n2cut
    ;
    ; Compute regional series for subset (based on local r)
    ;
    print,'Number of sites being used:',n2cut
    ;
    if n2cut gt 0 then begin
      iuse=n2cut-1
      ;
      ; Extract skillful subset and subset of the correlation matrix
      ;
      submxd=regmxd(slist(0:iuse),*)
      subfra=regfra(slist(0:iuse),*)
      subcorr=corrmat(slist(0:iuse),*)
      subcorr=subcorr(*,slist(0:iuse))
      subcorr=reform(subcorr,iuse+1,iuse+1)     ; make sure it is 2D
      ;
      ; Compute weighted regional mean with appropriate variance adjustment
      ;
      meanmxd=var_adjust(submxd,timey,weight=subfra,corrmat=subcorr,$
        refperiod=rbarper,/no_anomaly,/td_rbar,actualsample=sampsize)
      mknormal,meanmxd,timey,refperiod=rbarper
      ;
      bestmxd(*,ireg)=meanmxd(*)
      sampmxd(*,ireg)=sampsize(*)
      cutx(0:iuse,ireg)=xxx(slist(0:iuse))
      cuty(0:iuse,ireg)=yyy(slist(0:iuse))
      ;
    endif
    ;
  endif
  ;
endfor
;
; Now save results
;
save,filename='regts_'+fnadd+'_fixed'+string(fixedyr,format='(I4)')+'.idlsave',$
  nreg,regname,bestmxd,timey,fixedyr,maxntree,cutx,cuty,cutsites,sampsize
;
end