;
; First stage in computing timeseries of instrumental data for Harry's regions
;
; Get region info
;
restore,filename='harryboxes.idlsave'
g=def_grid(/ipcc)
;
; Define arrays
;
nyrmax=1992-1856+1
nmon=12
regts=fltarr(nyrmax*nmon,nreg)+!values.f_nan
regcover=fltarr(nyrmax*nmon,nreg)
;
; Now also want to keep groups of grid-box timeseries for the boxes that go
; into each regional mean
;
for i = 0 , nreg-1 do begin
  dummy=where(finite(boxlists(*,*,i)),dumn)
  if i eq 0 then begin
    nboxkp=dumn
    boxkp=dummy
  endif else begin
    nboxkp=nboxkp+dumn
    boxkp=[boxkp,dummy]
  endelse
  print,i,dumn,nboxkp
endfor
a=boxkp
b=a(uniq(a,sort(a)))          ; get rid of any duplicates
boxkp=b
nboxkp=n_elements(b)
print,'No. of boxes to keep=',nboxkp
boxts=fltarr(nboxkp,nyrmax*nmon)*!values.f_nan
;
print,'Opening netCDF file and reading 1856-1992 temperature data'
ntime=(1992-1856+1)*nmon
elest=(1856-1856)*nmon
eleen=(1992-1856)*nmon+nmon-1
ncid=ncdf_open('/cru/dave3/f055/detect/obsdat/tair_monthly.nc')
ncdf_attget,ncid,'temperature','missing_value',missv
ncdf_varget,ncid,'temperature',fdall,$
  offset=[0.,0.,0.],count=[g.nx,g.ny,ntime]
print,'Altering the missing values to NaN'
misslist=where(fdall eq missv,nmiss)
if nmiss gt 0 then fdall(misslist)=!values.f_nan
print,'Computing regional means'
for ireg = 0 , nreg-1 do begin
  print,'region=',ireg
  print,'Making mask'
  fdmask=boxlists(*,*,ireg)
  print,'Finding mean'
  regts(elest:eleen,ireg)=globalmean(fdall,g.y,mask=fdmask,cover=cover)
  regcover(elest:eleen,ireg)=cover
endfor
print,'Now extracting individual timeseries'
fd=reform(fdall,g.nx*g.ny,ntime)
boxts(*,elest:eleen)=fd(boxkp,*)
ncdf_close,ncid
;
; Now save them
;
x=findgen(nyrmax)+1856.
for ireg = 0 , nreg-1 do begin
  ; extract one region and separate into months, years and then transpose into
  ; year,months
  winstr=reform(regts(*,ireg),nmon,nyrmax)
  winstr=transpose(winstr)
  wcover=reform(regcover(*,ireg),nmon,nyrmax)
  wcover=transpose(wcover)
  save,filename='harry.instrts_'+regname(ireg)+'.idlsave',x,winstr,wcover
endfor
allx=fltarr(g.nx,g.ny)
for i = 0 , g.ny-1 do allx(*,i)=g.x
ally=fltarr(g.nx,g.ny)
for i = 0 , g.nx-1 do ally(i,*)=g.y
boxlon=allx(boxkp)
boxlat=ally(boxkp)
boxts=reform(boxts,nboxkp,nmon,nyrmax)
save,filename='harry.instrboxes.idlsave',x,boxts,boxlon,boxlat
;
end