;
multi_plot,nrow=4,ncol=2,layout='large'
if !d.name eq 'X' then begin
  window,xsize=650,ysize=950
  !p.font=-1
endif else begin
  !p.font=0
  device,/helvetica,/bold,xsize=17
endelse
loadct,0
;
; Plots maps of chronology location
;
ncid=ncdf_open('tree_dens_nh.nc')
ncdf_diminq,ncid,'station',dummy,nstat
ncdf_varget,ncid,'country',country
ncdf_varget,ncid,'tree',tree
ncdf_varget,ncid,'year',xtime
ncdf_varget,ncid,'density',dens
ncdf_attget,ncid,'density','missing_value',valmiss
ncdf_varget,ncid,'latitude',ylat
ncdf_varget,ncid,'longitude',xlon
ncdf_close,ncid
misslist=where(dens eq valmiss,nmiss)
dens(misslist)=!values.f_nan
;
restore,filename='reglists.idlsave'
iall=where(regname eq 'ALL')
i=iall(0)
iwant=treelist(0:ntree(i)-1,i)
;
map=def_map(/npolar)  &  map.limit(0)=30.
coast=def_coast(/get_device)  &  coast.double=0
labels=def_labels(/off)
;
for i = 0 , 1 do begin
;
inter_boxfd,/nodata,map=map,coast=coast,labels=labels
;
if i eq 0 then yrwant=where(xtime eq 1700.) else yrwant=where(xtime eq 1900.)
dvals=reform(dens(yrwant,*))
dvals=dvals(iwant)
igot=iwant(where(finite(dvals)))
x=xlon(igot)
y=ylat(igot)
cpl_usersym,/circle,/fill
plots,x,y,psym=8,noclip=0,symsize=0.6
cpl_usersym,/circle
plots,x,y,psym=8,thick=1,noclip=0,color=!p.background,symsize=0.6
;
device,font_size=9
!p.ticklen=0.03
lon_polar,map=map,[-180,-150.,-120.,-90,-60.,-30.,0.,30.,60.,90,120.,150.],/dotick
device,font_size=14
;
endfor
;
;
; Reset to single column
;
!p.multi(1)=1
!p.multi(0)=3
;
; 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.
;
; Get adjusted regional timeseries
;
restore,filename='reglists.idlsave'
;
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 region
;
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
;
!y.margin=[0,3]
;
filter_cru,20,tsin=meanmxd,tslow=tslow,/nan 
cpl_barts,x,meanmxd,$
  /xstyle,xtickformat='nolabels',$
  ytitle='        Normalised mean density anomaly',$
  zeroline=tslow,yrange=[-7.5,3],ystyle=9
oplot,x,tslow,thick=2 
oplot,!x.crange,[0.,0.]
;
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)
device,font_size=12
for i = 0 , nlab-1 do xyouts,yrx(i),yry(i),/data,yrlab(i),size=0.5,align=0.5
device,font_size=14
;
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)'
;
!y.margin=[16,0]
;
xsk=[1400,1450,1500,1550,1600,1620,1640,1660,1680,1700,1720,1740,1760,1780,$
  1800,1820,1840,1860,1880,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]
plot,xsk,ysk,/xstyle,xtitle='Year',$
  /ystyle,yrange=[0.,0.9],ytitle='Correlation'
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]
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
;
end