;
; Reads in density data, puts a low-frequency filter through each retained
; grid-box, and plots anomaly (normalised w.r.t. 1901-40) map
;
; Get display ready
;
yrwant=1980
yrwantst=string(yrwant,format='(I4)')
thalf=10.
;
loadct,39
multi_plot,nrow=2,layout='large'
if !d.name eq 'X' then begin
  window,ysize=900,xsize=650
  fac=1.
endif else begin
  fac=2.
endelse
;
; N. Hemi. map down to 30N, no labels
;
map=def_map(/npolar)
map.limit(0)=30.
labels=def_labels(/off)
labels.gridon=1  &  labels.dlon=90.
coast=def_coast(/get_device)
coast.double=0
;
; 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)
;
fns='instr'+yrwantst+'.idlsave'
dummy=findfile(fns,count=nfile)
if nfile gt 0 then begin
  restore,filename=fns
endif else begin
;
; For each available box, normalise w.r.t. 1901-40, and filter with thalf
;
statvals=fltarr(nstat)
iwant=where(x eq yrwant)
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)
  dummy=where((totn gt 0) and (x ge 1901.) and (x le 1940.),nkkkk)
  print,i,nkkkk
  if nkkkk gt 10 then begin
    keeplist=where(totn gt 0,nkeep)
    ts2=totval*!values.f_nan
    ts2(keeplist)=totval(keeplist)/float(totn(keeplist))
    mknormal,ts2,x,refperiod=[1901,1940]
    filter_cru,thalf,tsin=ts2,tslow=ts3,/nan
    statvals(i)=ts3(iwant)
  endif else begin
    statvals(i)=!values.f_nan
  endelse
endfor
save,filename=fns,statvals
endelse
;
; Remove missing data and data south of 30N
;
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
;
; 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
;
fns='dens'+yrwantst+'.idlsave'
dummy=findfile(fns,count=nfile)
if nfile gt 0 then begin
  restore,filename=fns
endif else begin
;
; For each available box, normalise w.r.t. 1901-40, and filter with thalf
;
statdens=fltarr(nstat)
iwant=where(x eq yrwant)
fracdens=frac(iwant,*)
print,'nstat=',nstat
for i = 0 , nstat-1 do begin
  ts1=reform(density(*,i))
  dummy=where(finite(ts1) and (x ge 1901.) and (x le 1940.),nkkkk)
  print,i,nkkkk
  if nkkkk gt 10 then begin
    mknormal,ts1,x,refperiod=[1901,1940]
    filter_cru,thalf,tsin=ts1,tslow=ts3,/nan
    statdens(i)=ts3(iwant)
  endif else begin
    statdens(i)=!values.f_nan
  endelse
endfor
save,filename=fns,statdens,fracdens
endelse
;
; 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)
;
levels=findgen(16)*0.4-3.
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,g,b,10,color='purple'
def_1color,r,g,b,16,color='lgreen'
def_smearcolor,r,g,b,fromto=[10,16]
def_1color,r,g,b,17,color='lgrey'
def_1color,r,g,b,18,color='lyellow'
def_1color,r,g,b,24,color='red'
def_smearcolor,r,g,b,fromto=[18,24]
levcol=indgen(15)+10
;
scale_vert,levels=levels,c_colors=levcol
;
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=yrwantst+' (20 yr) normalised density - temperature difference'
;
contour,fdsm.fd,fdsm.x,fdsm.y,levels=levels,/overplot
;
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
;
fns='rwid'+yrwantst+'.idlsave'
dummy=findfile(fns,count=nfile)
if nfile gt 0 then begin
  restore,filename=fns
endif else begin
;
; For each available box, normalise w.r.t. 1901-40, and filter with thalf
;
statrwid=fltarr(nstat)
iwant=where(x eq yrwant)
fracrwid=frac(iwant,*)
print,'nstat=',nstat
for i = 0 , nstat-1 do begin
  ts1=reform(density(*,i))
  dummy=where(finite(ts1) and (x ge 1901.) and (x le 1940.),nkkkk)
  print,i,nkkkk
  if nkkkk gt 10 then begin
    mknormal,ts1,x,refperiod=[1901,1940]
    filter_cru,thalf,tsin=ts1,tslow=ts3,/nan
    statrwid(i)=ts3(iwant)
  endif else begin
    statrwid(i)=!values.f_nan
  endelse
endfor
save,filename=fns,statrwid,fracrwid
endelse
;
; 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=yrwantst+' (20 yr) normalised width - temperature difference'
;
contour,fdsm.fd,fdsm.x,fdsm.y,levels=levels,/overplot
;
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
;
end