;
; Reads in density data-instrumental, computes two decadal means and
; differences them. Plots difference (normalised w.r.t. 1881-1975) map
;
; Get display ready
;
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]
;
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)
;
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)
;
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=titst+' density - temperature difference (normalised)'
;
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
;
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=titst+' width - temperature difference (normalised)'
;
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