;
; Create an Apr-Sep mean timeseries of gridded IST from monthly observed 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 I allow just 3 months of obs to be enough to get a AMJJAS mean, and then
; 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, year-by-year timeseries are considered.  Any
; missing values that have non-missing values before and after are infilled
; by a mean of those values.
; Also plots the long-term s.d.
;
; CAN OPTIONALLY COMPUTE THE JJA TEMP INSTEAD OF AMJJAS TEMP
;
dojja=0        ; 0=Apr-Sep 1=Jun-Aug
dosave=0       ; 0=do not, 1=do save the resulting data set
;
; Read in the mask of where the MXD gridded data are, so I can check what
; data is available over that region only.
;
print,'Getting MXD coverage'
restore,'mxd_infill.idlsave'
numval=total(finite(mxdfd),3)
mxdlist=where(numval gt 0,nmxd)
;
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=[1998,11],grid=g,info=i)
ncdf_close,ncid
;
; Compute Apr-Sep means from the data (requiring 3 or more months of data)
; or Jun-Aug means from the data (requiring 2 or more months of data)
;
print,'Computing seasonal means'
nmon=12
nyr=g.nt/nmon
fdmon=reform(fdmon,g.nx*g.ny,nmon,nyr)
if dojja eq 0 then begin
  print,'April-September'
  fdseas=mkseason(fdmon,3,8,datathresh=3)
  fullseas=mkseason(fdmon,3,8,datathresh=6)
endif else begin
  print,'June-August'
  fdseas=mkseason(fdmon,5,7,datathresh=2)
  fullseas=mkseason(fdmon,5,7,datathresh=3)
endelse
timey=findgen(nyr)+g.year(0)
;
; Check how much extra data we have over the MXD region by allowing
; incomplete seasonal means
;
testper=where((timey ge 1911) and (timey le 1990),ntest)
maskfd=fdseas[mxdlist,*]
maskfd=maskfd[*,testper]
dummy=where(finite(maskfd),nvalfd)
maskfd=fullseas[mxdlist,*]
maskfd=maskfd[*,testper]
dummy=where(finite(maskfd),nvalfull)
print
print,'Extra MXD region 1911-1990 values obtained by allowing incomplete seasonal data'
print,nvalfull,nvalfd,ntest*nmxd
print,float([nvalfull,nvalfd])/(ntest*nmxd)
print
;
fdseas=reform(fdseas,g.nx,g.ny,nyr)
xlon=g.x
ylat=g.y
nx=g.nx
ny=g.ny
;
; Now do some infilling if possible from surrounding values
;
print,'Spatial infilling',nyr
for i = 0 , nyr-1 do begin
  print,i,format='($,I4)'
  onefd=reform(fdseas(*,*,i))
  fdinf=tim_infill(onefd,3.,/cyclic,thresh=0.42,ylat=ylat)
  fdseas(*,*,i)=fdinf(*,*)
endfor
print
;
; Check how much extra data we have over the MXD region by allowing
; spatial infilling
;
maskfd=(reform(fdseas,g.nx*g.ny,nyr))[mxdlist,*]
maskfd=maskfd[*,testper]
dummy=where(finite(maskfd),nvalfd)
print
print,'Extra MXD region 1911-1990 values obtained by spatial infilling'
print,nvalfull,nvalfd,ntest*nmxd
print,float([nvalfull,nvalfd])/(ntest*nmxd)
print
;
; Now do some infilling if possible from previous/subsequent values
;
print,'Temporal infilling'
newfd=fdseas
for ix = 0 , nx-1 do begin
  print,ix
  for iy = 0 , ny-1 do begin
    print,iy,format='($,I3)'
    for iyr = 1 , nyr-2 do begin
      if finite(newfd(ix,iy,iyr)) eq 0 then begin
        newfd(ix,iy,iyr)=0.5*(fdseas(ix,iy,iyr-1)+fdseas(ix,iy,iyr+1))
      endif
    endfor
  endfor
  print
endfor
fdseas=newfd
;
; Check how much extra data we have over the MXD region by allowing
; temporal infilling
;
maskfd=(reform(fdseas,g.nx*g.ny,nyr))[mxdlist,*]
maskfd=maskfd[*,testper]
dummy=where(finite(maskfd),nvalfd)
print
print,'Extra MXD region 1911-1990 values obtained by temporal infilling'
print,nvalfull,nvalfd,ntest*nmxd
print,float([nvalfull,nvalfd])/(ntest*nmxd)
print
;
; 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)
fdseas=fdseas(*,ykeep,*)
landonly=landonly(*,ykeep)
seaonly=seaonly(*,ykeep)
mainlysea=mainlysea(*,ykeep)
;
; Now compute long-term mean and s.d.
;
fdd=double(fdseas)
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 year-by-year timeseries of the fraction of missing data
;
nall=total(finite(fdltm))
nmiss=float(total(finite(fdseas),1))
nmiss=1.-(total(nmiss,1)/nall)
missfrac=nmiss
;
if dojja eq 0 then fn='obs_temp_as.idlsave' else fn='obs_temp_jja.idlsave'
if dosave ne 0 then begin
  save,filename=fn,timey,fdseas,xlon,ylat,nyr,nx,ny,$
    fdltm,fdsd,landonly,seaonly,mainlysea,missfrac
endif
;
loadct,39
multi_plot,nrow=2
if !d.name eq 'X' then window,ysize=850
;
map=def_map(/npolar)  &  map.limit(0)=15.
coast=def_coast(/get_device)
labels=def_labels(/off)
;
if dojja eq 0 then tit='Apr-Sep' else tit='Jun-Aug'
levels=findgen(17)*0.1
inter_boxfd,fdsd,xlon,ylat,$
  /scale,$
  title='Long-term-standard deviation '+tit+' temperature',$
  map=map,coast=coast,labels=labels,$
  levels=levels
;
plot,timey,missfrac,psym=10,/xstyle,xtitle='Year',$
  ytitle='Fraction of data missing per summer',yrange=[0.,0.5],/ystyle
;
end