;
g=def_grid(/ipcc)
;
multi_plot,nrow=2,ncol=1,layout='large'
if !d.name eq 'X' then begin
  window,ysize=850
  fac=1.
endif else begin
  fac=2.
endelse
!x.omargin=[0,10]
!y.omargin=[0,0]
!x.margin=[0,1]
!y.margin=[0,0]
;
loadct,39
def_1color,r,gr,b,100,color='lgrey'
;
map=def_map(/npolar)  &  map.limit(0)=30.
map.xmargin=[0,1]
map.ymargin=[0,0]
labels=def_labels(/off)  &  labels.gridon=1  &  labels.dlon=30.
labels.dlat=10.
;
; Reads in density data-instrumental, computes two decadal means and
; differences them. Plots difference (normalised w.r.t. 1881-1975) map
;
navg=[10,10]
yrwant=[1935,1975]
aaa=string(yrwant,format='(I4)')
bbb=string(yrwant+navg-1,format='(I4)')
titst=aaa(1)+'-'+bbb(1)+' minus '+aaa(0)+'-'+bbb(0)
normper=[1881,1975]
;
; Small amount of smoothing, with filling in around the edges to extend to
; region that is plotted
;
sm=def_sm()
sm.method=2
sm.thresh=1.0
;
; Read in data
;
restore,filename='nc.instrboxes.idlsave'
nstat=n_elements(boxlon)
;
for iii = 0 , 1 do begin
  ;
  fns='instr'+aaa(iii)+'.idlsave'
  print,fns
  dummy=findfile(fns,count=nfile)
  if nfile gt 0 then begin
    restore,filename=fns
  endif else begin
    ;
    ; For each available box, normalise and compute decadal mean
    ;
    statvals=fltarr(nstat)
    iwant=where( (x ge yrwant(iii)) and (x lt yrwant(iii)+navg(iii)) )
    print,'iwant=',iwant
    print,'nstat=',nstat
    for i = 0 , nstat-1 do begin
      ts1=reform(boxts(i,*,*))              ; to give a mon , yr array
      totval=total(ts1(3:8,*),1,/nan)
      totn=total(finite(ts1(3:8,*)),1)
      keeplist=where(totn gt 0,nkeep)
      ts2=totval*!values.f_nan
      if nkeep gt 0 then $
        ts2(keeplist)=totval(keeplist)/float(totn(keeplist))
      dummy=where(finite(ts2) and (x ge normper(0)) and (x le normper(1)),nkk)
      print,i,nkk,format='($,I4,I3)'
      if nkk gt 10 then begin
        mknormal,ts2,x,refperiod=normper
        ts3=ts2(iwant)
        statvals(i)=total(ts3,/nan)/total(finite(ts3))
      endif else begin
        statvals(i)=!values.f_nan
      endelse
    endfor
    save,filename=fns,statvals
  endelse
  if iii eq 0 then keepinstr=statvals $
  else statvals=statvals-keepinstr
endfor
;
; Remove missing data and data south of 30N
;
blon1=boxlon
blat1=boxlat
keeplist=where(finite(statvals) and (boxlat ge 30.),nkeep)
if nkeep gt 0 then begin
   statvals=statvals(keeplist)
   boxlat=boxlat(keeplist)
   boxlon=boxlon(keeplist)
endif
print,'No. of boxes with data:',nkeep
blon2=boxlon
blat2=boxlat
;
; Read in data
;
ncid=ncdf_open('tree_dens_nh.nc')
ncdf_diminq,ncid,'time',dummy,ntime
ncdf_diminq,ncid,'station',dummy,nstat
ncdf_varget,ncid,'year',x
ncdf_varget,ncid,'density',density
ncdf_attget,ncid,'density','missing_value',valmiss
ncdf_varget,ncid,'fraction',frac
ncdf_varget,ncid,'longitude',statlon
ncdf_varget,ncid,'latitude',statlat
ncdf_close,ncid
misslist=where(density eq valmiss,nmiss)
density(misslist)=!values.f_nan
;
for iii = 0 , 1 do begin
  ;
  fns='dens'+aaa(iii)+'.idlsave'
  print,fns
  dummy=findfile(fns,count=nfile)
  if nfile gt 0 then begin
    restore,filename=fns
  endif else begin
    ;
    ; For each available site, normalise and compute decadal mean
    ;
    statdens=fltarr(nstat)
    iwant=where( (x ge yrwant(iii)) and (x lt yrwant(iii)+navg(iii)) )
    f2=frac(iwant,*)
    fracdens=total(f2,1,/nan)/total(finite(f2),1)
    print,'iwant=',iwant
    print,'nstat=',nstat
    for i = 0 , nstat-1 do begin
      ts1=reform(density(*,i))
      dummy=where(finite(ts1) and (x ge normper(0)) and (x le normper(1)),nkk)
      print,i,nkk,format='($,I4,I3)'
      if nkk gt 10 then begin
        mknormal,ts1,x,refperiod=normper
        ts3=ts1(iwant)
        statdens(i)=total(ts3,/nan)/total(finite(ts3))
      endif else begin
        statdens(i)=!values.f_nan
      endelse
    endfor
    save,filename=fns,statdens,fracdens
  endelse
  if iii eq 0 then keepdens=statdens $
  else statdens=statdens-keepdens
endfor
;
; Remove missing data and data south of 30N
;
keeplist=where(finite(statdens) and (statlat ge 30.),nkeep)
if nkeep gt 0 then begin
   statdens=statdens(keeplist)
   statlat=statlat(keeplist)
   statlon=statlon(keeplist)
   fracdens=fracdens(keeplist)
endif
print,'No. of boxes with data:',nkeep
;
; First of all, store each station value into its 0.5 by 0.5 grid box.
; When more than one falls in a box, average them - this is where the
; weighting by the fraction of cores available comes into it!  It also
; prevents duplicate points going forward, which can upset spherical
; triangulation.
;
dx=5.  &  dy=5.
gnx=360./dx
gny=(90.-30.)/dy
gx=findgen(gnx)*dx-177.5
gy=87.5-findgen(gny)*dy
griddens=gridit(gnx,gny,gx,gy,statlon,statlat,statdens,fracdens,nstat=chron)
;
; Convert boxes back to a list of stations
;
gx2d=gx # (intarr(gny)+1.)
gy2d=transpose(gy # (intarr(gnx)+1.))
keeplist=where(finite(griddens),nkeep)
griddens=griddens(keeplist)
gx=gx2d(keeplist)
gy=gy2d(keeplist)
;
; Now difference the fields
;
griddiff1=griddens*!values.f_nan
print,'No. of density boxes',nkeep
for i = 0 , nkeep-1 do begin
  dval=griddens(i)
  dxxx=gx(i)
  dyyy=gy(i)
  instrlist=where((boxlon eq dxxx) and (boxlat eq dyyy),ngot)
  if ngot gt 1 then message,'Strange - two boxes the same!'
  if ngot eq 1 then griddiff1(i)=dval-statvals(instrlist)
endfor
keeplist=where(finite(griddiff1),nkeep)
print,'No. of that have both',nkeep
if nkeep eq 0 then message,'Nothing to plot'
griddiff1=griddiff1(keeplist)
gx=gx(keeplist)
gy=gy(keeplist)
;
lll=[0.2 , 0.4 , 0.6 , 0.8 , 1.2 , 2.0]
levels=[-reverse(lll),lll]
nlev=n_elements(levels)
zeroloc=where(levels eq 0,nzero)
c_labels=intarr(nlev)+1  &  if nzero gt 0 then c_labels(zeroloc)=0
c_charsize=0.6
;c_thick=(levels lt 0.)*fac+2.  &  if nzero gt 0 then c_thick(zeroloc)=1.
c_thick=1.5*fac  ; &  if nzero gt 0 then c_thick(zeroloc)=1.
;print,levels
;print,c_thick
def_1color,r,gr,b,10,color='purple'
def_1color,r,gr,b,11,color='dblue'
def_1color,r,gr,b,12,color='blue'
def_1color,r,gr,b,14,color='lgreen'
def_smearcolor,r,gr,b,fromto=[12,14]
def_1color,r,gr,b,15,color='lgrey'
def_1color,r,gr,b,16,color='lyellow'
def_1color,r,gr,b,20,color='red'
def_smearcolor,r,gr,b,fromto=[16,20]
levcol=indgen(11)+10
;
;scale_vert,levels=levels,c_colors=levcol
;
labels.gridon=0
;
inter_const,griddiff1,gx,gy,map=map,$
  maxdist=3000.,$
  gs=[2.,2.],$
  sm=sm,labels=labels,$
;  shade=1,sh_thresh=0.,sh_grey=[235,180],miss_grey='white',$
;  levels=levels,c_labels=c_labels,c_charsize=c_charsize,c_thick=c_thick,$
;  /follow,$
  miss_grey='white',$
  levels=levels,c_colors=levcol,$
  /cell_fill,fdsm=fdsm
;  title='!5Ring density decline!3'
;
;xyouts,/normal,!x.region(0)+0.05,!y.region(1)-0.1,'!5a!3',$
;  charsize=!p.charsize*2.5
;
contour,fdsm.fd,fdsm.x,fdsm.y,levels=levels,/overplot,$
  c_colors=replicate(!d.n_colors-1,nlev)
;
xcir=findgen(361)-180.
ycir=replicate(30.,361)
map_plots,xcir,ycir
;
x1=[180.,131.,75.,38.,-40.,-99.]
y1=[30.,30.,30.,30.,30.,30.]
x2=[180.,131.,75.,38.,-40.,-99.]
y2=[90.,90.,90.,90.,90.,90.]
n=n_elements(x1)
for i = 0 , n-1 do begin
  map_plots,[x1(i),x2(i)],[y1(i),y2(i)],thick=5
endfor
;
xs=findgen(79)-40.
ys=replicate(50.,79)
map_plots,xs,ys,thick=5
;
xs=findgen(82)-180.
ys=replicate(56.,82)
map_plots,xs,ys,thick=5
;
; 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,'latitude',ylat
ncdf_varget,ncid,'longitude',xlon
ncdf_close,ncid
;
restore,filename='reglists.idlsave'
iall=where(regname eq 'ALL')
i=iall(0)
;
x=xlon(treelist(0:ntree(i)-1,i))
y=ylat(treelist(0:ntree(i)-1,i))
cpl_usersym,/circle,/fill
plots,x,y,psym=8,color=!p.background,symsize=1.05,noclip=0
cpl_usersym,/circle
plots,x,y,psym=8,thick=1,symsize=1.05,noclip=0
;
regname(0:2)=['ESIB','CSIB','WSIB']
rnx=[170,95,50,-35,-33,-48,-128,-160]
rny=[40,45,42,53,33,33,75,36]
for ir = 0 , 7 do $
  xyouts,rnx(ir),rny(ir),regname(ir)
;
!x.omargin(1)=0
!y.margin=[1,1]
!p.multi(0)=2
scale_vert,levels=levels,c_colors=levcol,ytickformat='(F4.1)',$
  noextremes=['<','>']
scale_vert,/off
;
labels.gridon=1
!x.omargin(1)=10
!y.margin=[0,0]
!p.multi(0)=0
pause
;
;for i = 0 , nkeep-1 do begin
;  hx=0.5*dx
;  hy=0.5*dy
;  plots,[-hx,hx,hx,-hx,-hx]+gx(i),[-hy,-hy,hy,hy,-hy]+gy(i),color=!d.n_colors-2
;endfor
;
;scale_vert,/off
;
; Read in data
;
ncid=ncdf_open('tree_rwid_nh.nc')
ncdf_diminq,ncid,'time',dummy,ntime
ncdf_diminq,ncid,'station',dummy,nstat
ncdf_varget,ncid,'year',x
ncdf_varget,ncid,'width',density
ncdf_attget,ncid,'width','missing_value',valmiss
ncdf_varget,ncid,'fraction',frac
ncdf_varget,ncid,'longitude',statlon
ncdf_varget,ncid,'latitude',statlat
ncdf_close,ncid
misslist=where(density eq valmiss,nmiss)
density(misslist)=!values.f_nan
;
for iii = 0 , 1 do begin
  ;
  fns='rwid'+aaa(iii)+'.idlsave'
  print,fns
  dummy=findfile(fns,count=nfile)
  if nfile gt 0 then begin
    restore,filename=fns
  endif else begin
    ;
    ; For each available site, normalise and compute decadal mean
    ;
    statrwid=fltarr(nstat)
    iwant=where( (x ge yrwant(iii)) and (x lt yrwant(iii)+navg(iii)) )
    f2=frac(iwant,*)
    fracrwid=total(f2,1,/nan)/total(finite(f2),1)
    print,'iwant=',iwant
    print,'nstat=',nstat
    for i = 0 , nstat-1 do begin
      ts1=reform(density(*,i))
      dummy=where(finite(ts1) and (x ge normper(0)) and (x le normper(1)),nkk)
      print,i,nkk,format='($,I4,I3)'
      if nkk gt 10 then begin
        mknormal,ts1,x,refperiod=normper
        ts3=ts1(iwant)
        statrwid(i)=total(ts3,/nan)/total(finite(ts3))
      endif else begin
        statrwid(i)=!values.f_nan
      endelse
    endfor
    save,filename=fns,statrwid,fracrwid
  endelse
  if iii eq 0 then keeprwid=statrwid $
  else statrwid=statrwid-keeprwid
endfor
;
; Remove missing data and data south of 30N
;
keeplist=where(finite(statrwid) and (statlat ge 30.),nkeep)
if nkeep gt 0 then begin
   statrwid=statrwid(keeplist)
   statlat=statlat(keeplist)
   statlon=statlon(keeplist)
   fracrwid=fracrwid(keeplist)
endif
print,'No. of boxes with data:',nkeep
;
; First of all, store each station value into its 0.5 by 0.5 grid box.
; When more than one falls in a box, average them - this is where the
; weighting by the fraction of cores available comes into it!  It also
; prevents duplicate points going forward, which can upset spherical
; triangulation.
;
dx=5.  &  dy=5.
gnx=360./dx
gny=(90.-30.)/dy
gx=findgen(gnx)*dx-177.5
gy=87.5-findgen(gny)*dy
gridrwid=gridit(gnx,gny,gx,gy,statlon,statlat,statrwid,fracrwid,nstat=chron)
;
; Convert boxes back to a list of stations
;
gx2d=gx # (intarr(gny)+1.)
gy2d=transpose(gy # (intarr(gnx)+1.))
keeplist=where(finite(gridrwid),nkeep)
gridrwid=gridrwid(keeplist)
gx=gx2d(keeplist)
gy=gy2d(keeplist)
;
; Now difference the fields
;
griddiff2=gridrwid*!values.f_nan
print,'No. of width boxes',nkeep
for i = 0 , nkeep-1 do begin
  dval=gridrwid(i)
  dxxx=gx(i)
  dyyy=gy(i)
  instrlist=where((boxlon eq dxxx) and (boxlat eq dyyy),ngot)
  if ngot gt 1 then message,'Strange - two boxes the same!'
  if ngot eq 1 then griddiff2(i)=dval-statvals(instrlist)
endfor
keeplist=where(finite(griddiff2),nkeep)
print,'No. of that have both',nkeep
if nkeep eq 0 then message,'Nothing to plot'
griddiff2=griddiff2(keeplist)
gx=gx(keeplist)
gy=gy(keeplist)
;
;!p.multi(0)=1
;scale_vert,levels=levels,c_colors=levcol
;
inter_const,griddiff2,gx,gy,map=map,$
  maxdist=3000.,$
  gs=[2.,2.],$
  sm=sm,labels=labels,$
;  shade=1,sh_thresh=0.,sh_grey=[235,180],miss_grey='white',$
;  levels=levels,c_labels=c_labels,c_charsize=c_charsize,c_thick=c_thick,$
;  /follow,$
  miss_grey='white',$
  levels=levels,c_colors=levcol,$
  /cell_fill,fdsm=fdsm
;  title='!5Ring width decline!3'
;
;xyouts,/normal,!x.region(0)+0.05,!y.region(1)-0.1,'!5b!3',$
;  charsize=!p.charsize*2.5
;
contour,fdsm.fd,fdsm.x,fdsm.y,levels=levels,/overplot,$
  c_colors=replicate(!d.n_colors-1,nlev)
;
xcir=findgen(361)-180.
ycir=replicate(30.,361)
map_plots,xcir,ycir
;
;scale_vert,/off
;
restore,filename='regboxes.idlsave'
iall=where(regname eq 'ALL')
i=iall(0)
;
fd1=boxlists(*,*,i)
dummy=where(finite(fd1),nboxxx)
fdline=fd1*0.+4.
fdline(dummy)=2.
boxplot,fd1,g.x,g.y,highlight=fdline,/overmap,/overplot,thick=4
;
latname=['30','60','120','150','-30','-60','-120','-150']
nlat=n_elements(latname)
rnx=float(latname)
rny=[25,25.5,26.5,27,25,25,25.5,26.5]
for ir = 0 , nlat-1 do $
  xyouts,rnx(ir),rny(ir),latname(ir),align=0.5
;
!x.omargin(1)=0
!y.margin=[1,1]
!p.multi(0)=2
scale_vert,levels=levels,c_colors=levcol,ytickformat='(F4.1)',$
  noextremes=['<','>']
scale_vert,/off
;
end