;
; Creates MEAN timeseries from adjusted regional timeseries and adjusts the
; MEAN for changing sample size using a time-dependent rbar.
; Timeseries are weighted by their EPS when computing the mean.
;
doland=0         ; 0=portrait, 1=landscape
;
ploton,1,/color8
loadct,39
def_1color,50,color='red'
def_1color,51,color='deepblue'
def_1color,52,color='green'
def_1color,53,color='purple'
def_1color,54,color='orange'
if doland eq 0 then multi_plot,nrow=3,layout='large' $
               else multi_plot,nrow=1,/landscape
if !d.name eq 'X' then begin
  window,ysize=900
  !p.font=-1
  initx
endif else begin
  !p.font=0
  if doland eq 0 then device,/helvetica,/bold,font_size=12,xsize=17.9 $
                 else device,/helvetica,/bold,font_size=12,ysize=12.,$
                             xoffset=6.
endelse
;
; Get adjusted regional timeseries
;
restore,filename='reglists.idlsave'
nreg=nreg-2          ; not ARCTIC or EXTRA-ARCTIC
;
for i = 3 , nreg-1 do begin       ; not 0, 1 & 2 (ALL, NORTH & SOUTH)
  ;
  ; Get adjusted regional timseries
  ;
  restore,filename='densadj_'+regname(i)+'.idlsave'
  ;
  if i eq 3 then begin
    nyr=n_elements(x)
    ndata=nreg-3           ; not 0, 1 & 2
    allmxd=fltarr(nyr,ndata)
    alleps=fltarr(nyr,ndata)
  endif
  allmxd(*,i-3)=densadj(*)
  alleps(*,i-3)=denseps(*)
  ;
endfor
;
; Get grand correlation matrix between regions
;
restore,filename='rbartd_dens_MEAN.idlsave'
;
; Make a time-dependent rbar next
;
print,'Time-dependent rbar for ',nyr,'  years'
rbar_td=fltarr(nyr)*!values.f_nan
for it = 0 , nyr-1 do begin
  print,it,format='($,I4)'
  mask = where( finite(allmxd(it,*)) , nsamp)
  if nsamp gt 0 then begin
    if nsamp eq 1 then begin
      rbar_td(it)=1.               ; if n=1 then n'=1 regardless of rbar
    endif else begin
      ; get reduced grandr
      redr=grandr(mask,*)
      redr=redr(*,mask)
      ; Total it and average it, but only do i<>j.
      totr=0.
      totn=0.
      for j = 1 , nsamp-1 do begin
        for k = 0 , j-1 do begin
          totr=totr+redr(j,k)        ; should have no missing values I hope
          totn=totn+1.
          if redr(j,k) eq 1. then message,'Ooops!'
        endfor
      endfor
      rbar_td(it)=totr/totn
    endelse
  endif
endfor
print
;
; Now compute a weighted mean of the available regional timeseries
;
n=float(total(finite(allmxd),2))   ; compute timeseries of number of regions
totmxd=total(allmxd*alleps,2,/nan)
toteps=total(alleps,2,/nan)
meanmxd=totmxd/toteps
;
; Now need to adjust for varying sample size
;
xadj=mkeffective(meanmxd,n,rbar=rbar_td,neff=neff)
meanmxd=xadj
;
; Now normalise w.r.t. 1881-1960
;
mknormal,meanmxd,x,refperiod=[1881,1960]
;
; Now plot them
;
filter_cru,20,tsin=meanmxd,tslow=tslow,/nan 
cpl_barts,x,meanmxd,$
  ;title='Normalised adjusted MEAN of all regional MXD timeseries',$ 
  xtitle='Year',/xstyle,$ 
  ytitle='        Normalised mean density anomaly  (standard deviation units)',$
  yminor=4,yticks=4,ytickv=findgen(5)*2.-6.,$
  zeroline=tslow,yrange=[-9,3],ystyle=9,bar_color=[50,51]
oplot,x,tslow,thick=2
oplot,!x.crange,[0.,0.]
;
axis,xaxis=1,/xstyle           ;,xtitle='Year'
axis,xaxis=0,/xstyle,xtickformat='nolabels'
;
;xyouts,1420.,1.9,/data,'NHD1',size=1.5,color=54
;
yrlab=['1448','1453','1495','1544','1587','1601','1641/2/3','1666',$
  '1675','1695','1698/9','1740','1783','1816/7/8','1836/7','1884','1912']
yry  =[-3.1  ,-4.85 ,-3.03 ,-2.63 ,-3.32 ,-6.0  ,-5.       ,-3.42 ,$
  -3.81 ,-4.15 ,-3.52   ,-3.35 ,-3.05 ,-5.03     ,-3.5    ,-3.6  ,-4.1]
yrx  =[1440.5,1453  ,1495  ,1544  ,1587  ,1589  ,1642      ,1658  ,$
  1675  ,1695  ,1711    ,1740  ,1783  ,1817      ,1837    ,1884  ,1912]
nlab=n_elements(yrlab)
for i = 0 , nlab-1 do xyouts,yrx(i),yry(i),/data,yrlab(i),size=0.6,align=0.5
;
keepcol=!p.color
!p.color=53
yrvei=[1450,1452,1471,1477,1480,1482,1580,1586,1593,1600,1640,1641,1660,$
  1663,1667,1673,1680,1707,1739,1800,1815,1835,1853,1854,1883,1886,1902,$
  1907,1912,1932,1956,1980,1982,1991]
yrsiz=[5,6,5,5,5,5,6,5,5,6,5,6,6,5,5,5,5,5,5,5,7,5,5,5,6,5,6,5,6,5,5,5,5,6]
nvei=n_elements(yrvei)
ym=!y.crange(0)
for i = 0 , nvei-1 do begin
  z=yrvei(i)
  y=yrsiz(i)
  case y of
    5: plots,[z,z,z-1,z,z+1,z],[ym,ym+0.4,ym+0.4,ym+0.7,ym+0.4,ym+0.4]
    6: begin
      plots,[z,z],[ym,ym+1.1],thick=3
      polyfill,[z-2.,z,z+2.,z-2.],[ym+1.1,ym+1.6,ym+1.1,ym+1.1]
      end
    7: begin
      plots,[z,z],[ym,ym+2.2],thick=5
      polyfill,[z-3.,z,z+3.,z-3.],[ym+2.2,ym+2.7,ym+2.2,ym+2.2]
      end
  endcase
endfor
;
xyouts,1449.,ym+0.3,align=1.,'?',size=0.6
xyouts,1579.,ym+0.3,align=1.,'?',size=0.6
xyouts,1659.,ym+0.3,align=1.,'?',size=0.6
;
xyouts,1420.,ym+0.3,orient=90.,'VEI',size=0.9
plots,[1432.,1432.],[ym,ym+2.7],thick=3
plots,[1432.,1438.],[ym+2.7,ym+2.7],thick=3
plots,[1432.,!x.crange(1)],[ym+2.7,ym+2.7],linestyle=1
plots,[1432.,1438.],[ym+1.6,ym+1.6],thick=3
plots,[1432.,!x.crange(1)],[ym+1.6,ym+1.6],linestyle=1
plots,[1432.,1438.],[ym+0.7,ym+0.7],thick=3
plots,[1432.,!x.crange(1)],[ym+0.7,ym+0.7],linestyle=1
xyouts,1431.,ym+2.5,'7',size=0.6,align=1.
xyouts,1431.,ym+1.4,'6',size=0.6,align=1.
xyouts,1431.,ym+0.5,'5',size=0.6,align=1.
!p.color=keepcol
;
axis,yaxis=1,ystyle=1,yrange=!y.crange*0.116626,$         ;-0.105023,$
  yminor=4,yticks=5,ytickv=findgen(6)*0.2-0.8,$
  ytitle='        Estimated Northern Hemisphere temperature anomaly  (!Uo!NC)'
;
;plot,x,n,ytitle='Sample size',/xstyle
;
;plot,x,rbar_td,/xstyle,ytitle='rbar'
;
;plot,x,neff,/xstyle,ytitle='Effective sample size'
;
; Now save the weighted mean timeseries, for further analysis
;
densadj=meanmxd
save,filename='densadj_MEAN.idlsave',x,densadj,n,neff
;
; Also plot the fixed grid correlations below
;
xsk=[1400,1450,1500,1550,1600,1620,1640,1660,1680,1700,1720,1740,1760,1780,$
  1800,1820,1840,1860,1880,1901,1970,1980,1990,1994]
ysk=[0.32,0.37,0.33,0.34,0.35,0.38,0.39,0.40,0.46,0.46,0.48,0.50,0.52,0.54,$
  0.55,0.55,0.55,0.55,0.55,0.57,0.57,0.51,0.39,0.29]
;plot,xsk,ysk,/xstyle,xtitle='Year of fixed grid',$
;  /ystyle,yrange=[0.,1.],ytitle='Fixed-grid correlation over 1881-1960'
;ysk=[0.51,0.65,0.60,0.62,0.61,0.61,0.60,0.62,0.68,0.67,0.69,0.71,0.73,0.74,$
;  0.76,0.77,0.77,0.77,0.78,0.80,0.80,0.74,0.58,0.40]
;oplot,xsk,ysk,thick=5
;oplot,!x.crange,[0.2,0.2],linestyle=1
;oplot,!x.crange,[0.4,0.4],linestyle=1
;oplot,!x.crange,[0.6,0.6],linestyle=1
;oplot,!x.crange,[0.8,0.8],linestyle=1
;
plotoff,/gv
;
end