;
; Create a monthly mean timeseries of gridded IST.
; Only the northern hemisphere is kept.  Also, a land-sea mask is read
; in and the northerm hemisphere extracted.  Finally due to so much missing
; data, if any point is missing at a particular time but 3 or more surrounding
; points have values then it is infilled with the mean of those points.
; After this spatial infilling, temporal gaps are considered.  Any
; missing values that have non-missing values before and after are infilled
; by a mean of those values (in fact I now allow missing runs of up to 5
; values to be infilled by linear interpolation).
; Also plots the long-term s.d.
;
loadct,39
multi_plot,nrow=2
if !d.name eq 'X' then window,ysize=850
def_1color,20,color='red'
;
print,'Reading in IPCC global IST gridded data'
ncid=ncdf_open('/cru/u2/f055/data/obs/grid/surface/ist_ipcc_18561999.mon.nc')
fdmon=crunc_rddata(ncid,tst=[1856,0],tend=[1999,11],grid=g,info=i)
ncdf_close,ncid
;
xlon=g.x
ylat=g.y
nx=g.nx
ny=g.ny
;
nmon=12
nyr=g.nt/nmon
fdmon=reform(fdmon,nx,ny,nmon,nyr)
timey=findgen(nyr)+g.year(0)
;
; Now do some infilling if possible from surrounding values
;
print,'Smoothing',nyr
for i = 0 , nyr-1 do begin
  print,i
  for imon = 0 , nmon-1 do begin
    onefd=reform(fdmon(*,*,imon,i))
    fdinf=tim_infill(onefd,3.,/cyclic,thresh=0.42,ylat=ylat)
    fdmon(*,*,imon,i)=fdinf(*,*)
  endfor
endfor
;
; Now do some infilling if possible from previous/subsequent values
;
fdmon=reform(fdmon,nx,ny,nmon*nyr)
newfd=fdmon
for ix = 0 , nx-1 do begin
  for iy = 0 , ny-1 do begin
    ts1=reform(fdmon(ix,iy,*))
    kl=where(finite(ts1),nval)
    if nval gt 1 then begin
      kout=kl(0)
      km1=kl(0)     ; always remember the last index stored in output list
      dofill=0
      for i = 1 , nval-1 do begin
        k=kl(i)
        kgap=k-km1-1
        if (kgap gt 0) and (kgap lt 11) then begin
          dofill=1
          kout=[kout,indgen(kgap)+km1+1]
        endif
        kout=[kout,k]
        km1=k
      endfor
      print,ix,iy,dofill
      if dofill ne 0 then begin
        ts2=interpol(ts1(kl),kl,kout)
        newfd(ix,iy,kout)=ts2(*)
      endif
      if (ix eq 36) and (iy eq 10) then begin
        print,ts1,format='(10F8.2)'
        print,kl,format='(20I4)'
        print,kout,format='(20I4)'
        print,ts2,format='(10F8.2)'
        ots1=reform(fdmon(ix,iy,*))
        nts1=reform(newfd(ix,iy,*))
        plot,nts1(0:500),/xstyle
        oplot,ots1(0:500),color=20
        oele=kl & nele=kout
        plot,kl
        oplot,kout,color=20
      endif
    endif
  endfor
endfor
fdmon=reform(newfd,nx,ny,nmon,nyr)
;
; Now read in the land/sea fraction (1=land 0=ocean)
;
lsm=fltarr(nx,ny)
openr,1,'/cru/u2/f055/detect/obsdat/landsea.dat'
readf,1,format='(72I6)',lsm
close,1
lsm=lsm/100.
landonly=lsm*!values.f_nan
landonly(where(lsm eq 1.))=1.
seaonly=lsm*!values.f_nan
seaonly(where(lsm eq 0.))=1.
mainlysea=lsm*!values.f_nan
mainlysea(where(lsm lt 0.25))=1.
;
; Just keep the northern hemisphere stuff
;
ykeep=where(ylat ge 0.,ny)
ylat=ylat(ykeep)
fdmon=fdmon(*,ykeep,*,*)
landonly=landonly(*,ykeep)
seaonly=seaonly(*,ykeep)
mainlysea=mainlysea(*,ykeep)
;
; Now compute long-term mean and s.d.
;
fdd=double(reform(fdmon,nx,ny,nmon*nyr))
fdtot=total(fdd,3,/nan)
fdnum=float(total(finite(fdd),3))
fdltm=fdtot/fdnum
ml=where(fdnum lt 20.,nmiss)
if nmiss gt 0 then fdltm(ml)=!values.f_nan
;
fdsqu=total(fdd^2,3,/nan)
fdsqu=fdsqu/fdnum
fdsd=sqrt(fdsqu-fdltm^2)
;
; Now compute a monthly timeseries of the fraction of missing data
;
nall=total(finite(fdltm))
nmiss=float(total(finite(fdmon),1))
nmiss=1.-(total(nmiss,1)/nall)
missfrac=nmiss
;
fn='obs_temp_mon.idlsave'
save,filename=fn,timey,fdmon,xlon,ylat,nyr,nx,ny,$
  fdltm,fdsd,landonly,seaonly,mainlysea,missfrac,nmon
;
map=def_map(/npolar)  &  map.limit(0)=0.
coast=def_coast(/get_device)
labels=def_labels(/off)
;
tit='monthly'
levels=findgen(17)*0.4
inter_boxfd,fdsd,xlon,ylat,$
  /scale,$
  title='Long-term-standard deviation '+tit+' temperature',$
  map=map,coast=coast,labels=labels,$
  levels=levels
;
xxx=findgen(nmon*nyr)/float(nmon)+timey(0)
plot,xxx,missfrac,psym=10,/xstyle,xtitle='Year',$
  ytitle='Fraction of data missing per month',yrange=[0.,1],/ystyle
;
end