;
; Creates composite regional timeseries based on a limited sample of MXD
; series, with regional series variance adjusted with a time-dependent rbar,
; and then averaged with EPS-weighting and themselves variance adjusted.
;
; The MXD sites are limited first by specifying a year for which each site
; used must have had data for (e.g., specifying 1700 would use only sites
; with data back to at least 1700).  After this selection, all sites regardless
; of which sub-region they are in, that have correlations (with local Apr-Sep
; temperatures) >= CUTOFF are used, where CUTOFF is 0.6, 0.5, 0.4,
; 0.3, 0.22, 0.1, 0.0, -1.  *0.22 is used*
;
;-----------------------------------------------------------------
;
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                 ; cutoffs for local correlations
;
; 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
;
; Now define composite regions
;
ncomp=3
compname=['HILAT','LOLAT','ALL']
avgnum=[5,4,9]
avgreg=intarr(ncomp,max(avgnum))
avgreg(0,0:avgnum(0)-1)=[0,2,3,7,8]
avgreg(1,0:avgnum(1)-1)=[1,4,5,6]
avgreg(2,0:avgnum(2)-1)=indgen(nreg)
;
; 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
nsites=fltarr(nreg)       ; no. sites with data back to fixed yr & local r
cutx=fltarr(nchron,ncomp)*z    ; location of sites with local r >= 0.22
cuty=fltarr(nchron,ncomp)*z
compmxd=fltarr(nyr,ncomp)*z  ; MXD series, locr >= 0.22
;
; Now analyse each region in turn
;
for icomp = 0 , ncomp-1 do begin
  print,compname(icomp)
  ;
  ; For the MXD, we need to do each sub-region individually, before composite
  ; averaging
  ;
  ndoreg=avgnum(icomp)
  indts=fltarr(ndoreg,nyr)*!values.f_nan
  indeps=fltarr(ndoreg,nyr)*!values.f_nan
  for idoreg = 0 , ndoreg-1 do begin
    ireg=avgreg(icomp,idoreg)
    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))        ; 18/20=Apr-Sep/Jun-Aug 1=temperature
    xxx=statlon(tl)
    yyy=statlat(tl)
    ;
    ; Select those sites with local skill >= current cutoff
    ;
    kl=where(regcor ge cutr,nkeep)
    print,'Skill-based subset',nkeep
    if nkeep gt 0 then begin
      regmxd=regmxd(*,kl)
      regfra=regfra(*,kl)
      regcor=regcor(kl)
      xxx=xxx(kl)
      yyy=yyy(kl)
      ;
      ; 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
      nsites(ireg)=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 and full rbar
        ;
        fullrbar=mkrbar(regmxd,grandr=corrmat)
        ;
        ; Compute weighted regional mean with appropriate variance adjustment
        ;
        meanmxd=var_adjust(regmxd,timey,weight=regfra,corrmat=corrmat,$
          refperiod=rbarper,/no_anomaly,/td_rbar)
        mknormal,meanmxd,timey,refperiod=rbarper
        ;
        ; Store series for later compositing
        ;
        indts(idoreg,*)=meanmxd(*)
        ;
        ; Compute time-dependent EPS for later weighting (need time-dependent
        ; sample size and full rbar)
        ;
        nsamp=float(total(finite(regmxd),1))
        indeps(idoreg,*)=fullrbar / (fullrbar + (1.-fullrbar)/nsamp(*) )
        ;
      endif
    endif
    ;
    ; Repeat for next region
    ;
  endfor
  ;
  ; Check in case we have no series at all (all failed the local r cutoff)
  ;
  if total(finite(indts)) eq 0 then begin
    domiss=1
    meanmxd=fltarr(nyr)*!values.f_nan
  endif else begin
    domiss=0
    ;
    ; Now compute an EPS-weighted, variance-adjusted composite regional mean
    ;
    meanmxd=var_adjust(indts,timey,weight=indeps,$
      refperiod=rbarper,/no_anomaly,/td_rbar)
    mknormal,meanmxd,timey,refperiod=rbarper
    ;
  endelse
  ;
  ; Store composite regional mean series if required
  ;
  cutx(0:nx-1,icomp)=xxx(*)
  cuty(0:nx-1,icomp)=yyy(*)
  compmxd(*,icomp)=meanmxd(*)
  ;
  ; Repeat for next composite region
  ;
endfor
;
; Now save results
;
save,filename='compts_'+fnadd+'_fixed'+string(fixedyr,format='(I4)')+'.idlsave',$
  ncomp,compname,nsites,compmxd,timey,$
  fixedyr,nchron,cutx,cuty
;
end