;
; Create an Apr-Sep mean timeseries of gridded precip from monthly observed
; precip.  Only the northern hemisphere is kept.  Finally due to so much missing
; data I allow just 3 months of obs to be enough to get a Apr-Sep mean.
; Also plots the long-term s.d.
;
print,'Reading in Mike Hulmes gridded precip data'
ncid=ncdf_open('/cru/u2/f055/data/obs/grid/surface/precip_19001995.mon.nc')
fdmon=crunc_rddata(ncid,tst=[1922,0],tend=[1995,11],grid=g,info=i)
ncdf_close,ncid
;
; Compute Apr-Sep means from the data (requiring 3 or more months of data)
;
nmon=12
nyr=g.nt/nmon
fdmon=reform(fdmon,g.nx*g.ny,nmon,nyr)
fdseas=mkseason(fdmon,3,8,datathresh=3)
fdseas=reform(fdseas,g.nx,g.ny,nyr)
timey=findgen(nyr)+g.year(0)
;
; convert from mm/month to mm/day  (mean month length for DJFM is 30.25 days)
;
fdseas=fdseas/30.25
;
xlon=g.x
ylat=g.y
nx=g.nx
ny=g.ny
;
; Just keep the northern hemisphere stuff
;
ykeep=where(ylat ge 0.,ny)
ylat=ylat(ykeep)
fdseas=fdseas(*,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 30.,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
;
save,filename='obs_prec_as.idlsave',timey,fdseas,xlon,ylat,nyr,nx,ny,$
  fdltm,fdsd,missfrac
;
loadct,39
multi_plot,nrow=2
if !d.name eq 'X' then window,ysize=850
;
map=def_map(/npolar)  &  map.limit(0)=10.
coast=def_coast(/get_device)
labels=def_labels(/off)
;
def_1color,cr,cg,cb,50,color='lgrey'
;
levels=[findgen(5)*0.2,findgen(8)*0.5+1.,findgen(8)+6.,30.]
cols=indgen(21)+50
def_1color,cr,cg,cb,50,color='red'
def_1color,cr,cg,cb,53,color='orange'
def_1color,cr,cg,cb,56,color='lyellow'
def_1color,cr,cg,cb,59,color='lgreen'
def_1color,cr,cg,cb,62,color='dgreen'
def_1color,cr,cg,cb,66,color='blue'
def_1color,cr,cg,cb,70,color='purple'
def_smearcolor,cr,cg,cb,fromto=[50,53]
def_smearcolor,cr,cg,cb,fromto=[53,56]
def_smearcolor,cr,cg,cb,fromto=[56,59]
def_smearcolor,cr,cg,cb,fromto=[59,62]
def_smearcolor,cr,cg,cb,fromto=[62,66]
def_smearcolor,cr,cg,cb,fromto=[66,70]
inter_boxfd,fdltm,xlon,ylat,$
  /scale,$
  title='Long-term-mean Apr-Sep precipitation  (mm/day)',$
  map=map,coast=coast,labels=labels,$
  levels=levels,c_colors=cols
;
;loadct,39
levels=[findgen(21)*0.05,5.]
inter_boxfd,fdsd,xlon,ylat,$
  /scale,$
  title='Long-term-standard deviation Apr-Sep precipitation  (mm/day)',$
  map=map,coast=coast,labels=labels,$
  levels=levels,c_colors=cols
;
pause
plot,timey,missfrac,psym=10,/xstyle,xtitle='Year',$
  ytitle='Fraction of data missing per summer',yrange=[0.,0.5],/ystyle
;
end