;
; 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, and
; compares them with composite regional and near-hemispheric instrumental
; temperature series which are also computed here (BUT NOT variance adjusted
; since for the instrumental series I do not want to go through a 2-stage
; averaging process, and the number of boxes involved in a direct average
; is too many to variance adjust - and the adjustment would be extremely small
; anyway!).
;
; The instrumental record can be based on either all available land
; observations in the region, or only those with an MXD site within their grid
; box (in the latter case, the region is defined by the FULL sample of MXD
; data, not by the limited sample of MXD sites actually being used).
; The instrumental data is based on Apr-Sep means, although this is adjustable.
; In addition, the near-hemispheric mean of all LAND north of 20N is used.
;
; 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.
;
; Comparison is made by computing correlations over the period 1881-1960,
; both filtered and unfiltered.
;
;-----------------------------------------------------------------
;
dotem=1         ; 0=use CRUTEM1, 1=use CRUTEM2v
if dotem eq 0 then begin
  fntem='lat_jones_18512000.mon.nc'
  savtem=''
endif else begin
  fntem='crutem2v_18512001.mon.nc'
  savtem='ct2v_'
endelse
;
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
    seasst=3
    seasen=8
    dthr=3
    end
  1: begin
    fnadd='trw'
    iseas=20        ; use local r with Jun-Aug temp to select chronologies
    seasst=5
    seasen=7
    dthr=2
    end
endcase
titadd=strupcase(fnadd)
;
; Specify options
;
fixedyr=1988              ; year of fixed grid
rbarper=[1881,1960]       ; period for computing correlation matrix
corrper=[1881,1960]    ; period over which comparison with temperature is made
doreg=0                ; 0=use co-located temperature only, 1=use full region
cutr=[0.6,0.5,0.4,0.3,0.22,0.1,0.0,-1.]    ; cutoffs for local correlations
ncut=n_elements(cutr)
;
; Identify analysis period:
; For speed, analyse the minimum possible. Must cover all of rbarper and
; corrper, plus a bit more either side of corrper for filtering.  Does not
; need to cover the fixedyr.
;
analper=[1860,2000]
anallen=analper(1)-analper(0)+1
;
; 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','NH']
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'
fnlocal='datastore/'+savtem+fnadd+'_moncorr.idlsave'
print,fnlocal
restore,filename=fnlocal
;  allr,ncorr,nvar,nchron,varname,corrname,statlat,statlon,allp,moir,moir2
;
; Read in the instrumental temperatures for the required period
;
print,'Reading temperatures'
print,fntem
ncid=ncdf_open('/cru/u2/f055/data/obs/grid/surface/'+fntem)
tsmon=crunc_rddata(ncid,tst=[analper(0),0],tend=[analper(1),11],grid=gtemp)
ncdf_close,ncid
nmon=12
ntemp=gtemp.nt
nyrtemp=ntemp/nmon
yrtemp=reform(gtemp.year,nmon,nyrtemp)
yrtemp=reform(yrtemp(0,*))
;
; Read in the regional breakdown of grid boxes
;
print,'Reading temperature regions'
restore,filename='reg_mxdboxes.idlsave'
;  nreg,boxlists,boxnumbs,regname,g,boxregs,boxprec,boxtemp
;
; Define arrays for storing all regional results
;
z=!values.f_nan
nsites=fltarr(nreg,ncut)       ; no. sites with data back to fixed yr & local r
compr=fltarr(ncomp+1,ncut,3)*z ; composite raw high low correls for cutoffs
compmxd=fltarr(anallen,ncomp,2)*z  ; MXD series: locr >= 0.22; all
comptemp=fltarr(anallen,ncomp+1)*z ; composite temperature timeseries & NH
;
; Now analyse each region in turn
;
multi_plot,nrow=4,ncol=2,layout='large',/landscape
loadct,39
def_1color,20,color='red'
def_1color,21,color='blue'
if !d.name eq 'X' then window,xsize=1000,ysize=850
for icomp = 0 , ncomp-1 do begin
  print,compname(icomp)
  ;
  ; Extract monthly instrumental data for current COMPOSITE region
  ;
  ; Choose mask (co-located data only, or all temperatures in region)
  ireg=avgreg(icomp,0:avgnum(icomp)-1)
  if doreg eq 0 then fdmask=reform(boxlists(*,*,ireg)) $
                else fdmask=reform(boxregs(*,*,ireg))
  fdmask=total(finite(fdmask),3) ; if any regional mask has data, comp has too
  ; Count boxes to use
  blist=where(fdmask gt 0,nbox)
  ; Define array to hold monthly grid box timeseries from region
  regts=fltarr(nbox,ntemp)
  ; Fill the array
  for imon = 0 , ntemp-1 do begin
     fd1mon=reform(tsmon(*,*,imon))
     regts(*,imon)=fd1mon(blist)
  endfor
  ; Separate months from years
  regts=reform(regts,nbox,nmon,nyrtemp)
  ; Compute seasonal average (April-September)
  regseas=mkseason(regts,seasst,seasen,datathresh=dthr)
  ; Define an array to hold latitude weighting
  locy=fix(blist/g.nx)
  regwt=cos(g.y(locy)*!dtor)
  wt=fltarr(nbox,nyrtemp)
  for iyr = 0 , nyrtemp-1 do wt(*,iyr)=regwt(*)
  ; Compute a regional average (latitude-weighted, but no variance adjustment)
  totval=total(regseas*wt,1,/nan)
  totnum=total(float(finite(regseas))*wt,1,/nan)
  meanlat=totval/totnum
  comptemp(*,icomp)=meanlat(*)
  mknormal,meanlat,yrtemp,refperiod=rbarper
  ;
  ; For the 'ALL' region, also compute the northern hemisphere >20N land series
  ;
  if icomp eq 2 then begin
    ; First extract >20N rows
    kl=where(g.y ge 20.)
    ylat=g.y(kl)
    tsnorth=tsmon(*,kl,*)
    ; Compute latitude-weighted mean
    nhmon=globalmean(tsnorth,ylat)
    ; Compute Apr-Sep mean
    nhmon=reform(nhmon,nmon,nyrtemp)
    nhseas=mkseason(nhmon,seasst,seasen,datathresh=dthr)
    comptemp(*,icomp+1)=nhseas(*)
    mknormal,nhseas,yrtemp,refperiod=rbarper
  endif
  ;
  ; Repeat regional mean for different subsets according to local skills
  ;
  !p.multi(0)=0
  pause
  for icut = 0 , ncut-1 do begin
    print,icut,cutr(icut)
    ;
    ; For the MXD, we need to do each sub-region individually, before composite
    ; averaging
    ;
    ndoreg=avgnum(icomp)
    indts=fltarr(ndoreg,anallen)*!values.f_nan
    indeps=fltarr(ndoreg,anallen)*!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(icut),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,icut)=nx
        if nx gt 0 then begin
          regmxd=regmxd(*,kl)
          regfra=regfra(*,kl)
          regcor=regcor(kl)
          xxx=xxx(kl)
          yyy=yyy(kl)
          ;
          ; Reduce time dimension to minimum required
          ;
          kl2=where((timey ge analper(0)) and (timey le analper(1)),nyr)
          print,nyr,' : THIS MUST BE NON-ZERO!'
          regmxd=regmxd(kl2,*)
          regfra=regfra(kl2,*)
          yrmxd=timey(kl2)
          ;
          ; Transpose to give space then time
          ;
          regmxd=transpose(regmxd)
          regfra=transpose(regfra)
          ;
          ; Compute correlation matrix and full rbar
          ;
          fullrbar=mkrbar(regmxd,grandr=corrmat)
          if nx eq 1 then fullrbar=1.
          ;
          ; Compute weighted regional mean with appropriate variance adjustment
          ;
          meanmxd=var_adjust(regmxd,yrmxd,weight=regfra,corrmat=corrmat,$
            refperiod=rbarper,/no_anomaly,/td_rbar)
          mknormal,meanmxd,yrmxd,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(anallen)*!values.f_nan
    endif else begin
      domiss=0
      ;
      ; Now compute an EPS-weighted, variance-adjusted composite regional mean
      ;
      meanmxd=var_adjust(indts,yrmxd,weight=indeps,$
        refperiod=rbarper,/no_anomaly,/td_rbar)
      mknormal,meanmxd,yrmxd,refperiod=rbarper
      ;
    endelse
    ;
    ; Store composite regional mean series if required
    ;
    if cutr(icut) eq 0.22 then begin
      compmxd(*,icomp,0)=meanmxd(*)
    endif
    if cutr(icut) eq -1. then begin
      compmxd(*,icomp,1)=meanmxd(*)
    endif
    ;
    ; Correlate MXD with temperature
    ;
    if domiss eq 0 then begin
      r3=mkcorrelation(meanmxd,meanlat,yrmxd,filter=10,refperiod=corrper,$
        datathresh=10)
      compr(icomp,icut,*)=r3(*)
      plot,yrmxd,meanmxd,/xstyle,title=string(cutr(icut))+'  '+string(r3(0))
      oplot,yrtemp,meanlat,color=20
    endif else begin
      compr(icomp,icut,*)=!values.f_nan
    endelse
    ;
    ; If region is 'ALL' then show and correlate with NH temperature too
    ;
    if icomp eq 2 then begin
      oplot,yrtemp,nhseas,color=21
      r3=mkcorrelation(meanmxd,nhseas,yrmxd,filter=10,refperiod=corrper,$
        datathresh=10)
      compr(icomp+1,icut,*)=r3(*)
    endif
    ;
    ; Repeat for next skill cutoff
    ;
  endfor
  ;
  ; Repeat for next composite region
  ;
endfor
;
; Now save results
;
timey=yrmxd
if doreg eq 0 then fnadd2='' else fnadd2='full'
save,filename='datastore/'+savtem+fnadd2+'compbest_'+fnadd+'_fixed'+$
  string(fixedyr,format='(I4)')+'.idlsave',$
  ncomp,compname,nreg,regname,nsites,compr,compmxd,comptemp,timey,$
  fixedyr,corrper,doreg,nchron,cutr,ncut,seasst,seasen
;
end