;
; Plots regional and large-scale averages of correlations between MXD and
; monthly & seasonal precip. and temperature.
;
; Restore all correlations
;
restore,filename='mxd_moncorr.idlsave'
;  allr,ncorr,nvar,nchron,varname,corrname,statlat,statlon
;
; Restore regional breakdown
;
restore,filename='reg_mxdlists.idlsave'
;  ntree,treelist,nreg,regname
;
; For each region, average the site-by-site correlations together (and also
; find the range, by using the mean +- 2sd, unless < 10 sites in which case
; max and min are used)
;
regr=fltarr(nreg,ncorr,nvar)
regmin=fltarr(nreg,ncorr,nvar)
regmax=fltarr(nreg,ncorr,nvar)
for i = 0 , nreg-1 do begin
  n=ntree(i)
  oner=allr(treelist(0:n-1,i),*,*)
  totr=total(oner,1)/float(n)
  regr(i,*,*)=totr(*,*)
  for j = 0 , ncorr-1 do begin
    for k = 0 , nvar-1 do begin
      rrr=reform(oner(*,j,k))
      kl=where(finite(rrr),nrrr)
      rrr=rrr(kl)
      if nrrr lt 10 then begin
        regmin(i,j,k)=min(rrr)
        regmax(i,j,k)=max(rrr)
      endif else begin
        dummy=moment(rrr,sdev=sdev)
        regmin(i,j,k)=(regr(i,j,k)-2.*sdev) > (-1.)
        regmax(i,j,k)=(regr(i,j,k)+2.*sdev) < (1.)
      endelse
    endfor
  endfor
endfor
;
; Now generate three extra composite regions
;
newr=fltarr(nreg+3,ncorr,nvar)
newmin=fltarr(nreg+3,ncorr,nvar)
newmax=fltarr(nreg+3,ncorr,nvar)
newr(0:nreg-1,*,*)=regr(*,*,*)
newmin(0:nreg-1,*,*)=regmin(*,*,*)
newmax(0:nreg-1,*,*)=regmax(*,*,*)
regname=[regname,'HILAT','LOLAT','NH']
;
rlist=[0,2,3,7,8]
newr(nreg,*,*)=total(regr(rlist,*,*),1)/float(n_elements(rlist))
newmin(nreg,*,*)=total(regmin(rlist,*,*),1)/float(n_elements(rlist))
newmax(nreg,*,*)=total(regmax(rlist,*,*),1)/float(n_elements(rlist))
;
rlist=[1,4,5,6]
newr(nreg+1,*,*)=total(regr(rlist,*,*),1)/float(n_elements(rlist))
newmin(nreg+1,*,*)=total(regmin(rlist,*,*),1)/float(n_elements(rlist))
newmax(nreg+1,*,*)=total(regmax(rlist,*,*),1)/float(n_elements(rlist))
;
newr(nreg+2,*,*)=total(regr,1)/float(nreg)
newmin(nreg+2,*,*)=total(regmin,1)/float(nreg)
newmax(nreg+2,*,*)=total(regmax,1)/float(nreg)
;
regr=newr
regmin=newmin
regmax=newmax
nreg=nreg+3
;
; Now prepare for plotting
;
loadct,39
multi_plot,nrow=3,ncol=2,layout='caption'
if !d.name eq 'X' then begin
  window,ysize=800
  !p.font=-1
endif else begin
  !p.font=0
  device,/helvetica,/bold,font_size=14
endelse
def_1color,20,color='vlgrey'
def_1color,21,color='mgrey'
;
for i = 0 , nreg-1 do begin
  for ivar = 0 , nvar-1 do begin
    pause
    x=findgen(28)
    yy=reform(regr(i,*,ivar))
    yymin=reform(regmin(i,*,ivar))
    yymax=reform(regmax(i,*,ivar))
    y=fltarr(28)*!values.f_nan
    ymin=fltarr(28)*!values.f_nan
    ymax=fltarr(28)*!values.f_nan
    yloc=[indgen(16)+1,indgen(5)*2+18]
    y(yloc)=yy(*)
    ymin(yloc)=yymin(*)
    ymax(yloc)=yymax(*)
    plot,x,y,/nodata,$
      /ystyle,yrange=[-1.,1.],ytitle='Correlation',$
      xstyle=5,title=regname(i)+'  r(MXD,'+varname(ivar)+')'
    for j = -10,10,2 do oplot,!x.crange,replicate(float(j)*0.1,2),linestyle=1
    oplot,!x.crange,[0.,0.]
    oplot,[7.5,7.5],!y.crange,thick=4
    cpl_barts,x,y,outline=-1,/overplot,bar_color=[21,20],thick=2
    if !d.name eq 'PS' then device,font_size=8
    for j = 0 , ncorr-1 do begin
      x1=replicate(x(yloc(j)),2)
      y1min=ymin(yloc(j))
      y1max=ymax(yloc(j))
      plots,x1,y1min,psym=def_sym(12)
      plots,x1,y1max,psym=def_sym(11)
      x1=x(yloc(j))+0.28
      y1=!y.crange(0)-0.05
      if !p.multi(0) le 1 then xyouts,x1,y1,corrname(j),align=1,orient=90.
    endfor
    xyouts,7.,-0.85,align=1.,'Previous!Cyear'
    xyouts,8.,-0.85,align=0.,'Same!Cyear'
    if !d.name eq 'PS' then device,font_size=14
  endfor
endfor
;
end