;
; Computes regional mean timeseries of instrumental data on a monthly basis
;
;----------------------------------------------------
;
; Get region info
;
restore,filename='timevary.idlsave'
g=def_grid(/ipcc)
nboxkp=maxbox
print,'No. of boxes to keep=',nboxkp
;
; Define arrays
;
nyrmax=1992-1600+1
nmon=12
regts=fltarr(nyrmax*nmon,nreg)+!values.f_nan
regcover=fltarr(nyrmax*nmon,nreg)
;
print,'Opening netCDF file and reading 1881-1992 temperature data'
ntime=(1992-1881+1)*nmon
misstime=(1880-1856+1)*nmon
elest=(1881-1600)*nmon
eleen=(1992-1600)*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.,misstime],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=',regname(ireg)
  for i = 0 , maxtime-1 do begin      ; maxtime given in timevary.idlsave
    fdmask=fltarr(g.nx,g.ny)*!values.f_nan
    n=boxcount(i,ireg)
    if n gt 0 then begin
      keeplist=reform(boxlists(0:n-1,i,ireg))
      fdmask(keeplist)=1.
      mon12=indgen(nmon)+elest+(i*nmon)
      mon12temp=indgen(nmon)+(i*nmon)
  regts(mon12,ireg)=globalmean(fdall(*,*,mon12temp),g.y,mask=fdmask,cover=cover)
      regcover(mon12,ireg)=cover
      print,i,n,total(cover(3:8))/6.,format='($,I4,2I3)'
    endif
  endfor
endfor
;
; Now save them
;
x=findgen(nyrmax)+1600.
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='instrvary_'+regname(ireg)+'.idlsave',x,winstr,wcover
endfor
;
end