;
; Create an Apr-Sep mean timeseries of gridded MSLP from monthly observed MSLP.
; Also plots the long-term mean and s.d.
; USES THE GMSLP DATASET!!!
;
dolowres=0    ; 0=keep with 5 by 5, 1=degrade to 10 by 5
dolowpole=1   ; 0=keep all boxes, 1=subsample at higher latitudes
; This lower resolution data set is stored in addition to the full data
;
ncid=ncdf_open('/cru/u2/f055/data/obs/grid/surface/gmslp3.2_18711998.mon.nc')
fdmon=crunc_rddata(ncid,tst=[1922,0],tend=[1998,11],grid=g,info=i)
ncdf_close,ncid
;
; Compute April-September means from the data (requiring full data by default)
;
nmon=12
nyr=g.nt/nmon
fdmon=reform(fdmon,g.nx*g.ny,nmon,nyr)
fdseas=mkseason(fdmon,3,8)
fdseas=reform(fdseas,g.nx,g.ny,nyr)
timey=findgen(nyr)+g.year(0)
;
xlon=g.x
ylat=g.y
nx=g.nx
ny=g.ny
;
if dolowres eq 1 then begin
  ;
  ; Degrade the longitudinal resolution by sampling every other box
  ;
  nx=nx/2
  kl=indgen(nx)*2
  xlon=xlon(kl)
  fdseas=fdseas(kl,*,*)
  ;
endif
;
if dolowpole eq 1 then begin
  fdseasfull=fdseas
  ; set all polar values except one to be missing, set the one to the mean
  kl=where(ylat eq 90.,ngot)
  if ngot eq 1 then begin
    mv=total(fdseas(*,kl,*),1)/float(nx)
    fdseas(*,kl,*)=!values.f_nan
    fdseas(nx/2,kl,*)=mv(*)
  endif
  ; keep 4 values at 85N
  kl=where(ylat eq 85.,ngot)
  if ngot eq 1 then begin
    sl=indgen(4)*18
    fd2d=reform(fdseas(*,kl,*))
    fd2dm=fd2d*!values.f_nan
    fd2dm(sl,*)=fd2d(sl,*)
    fdseas(*,kl,*)=fd2dm(*,*)
  endif
  ; keep 8 values at 80N
  kl=where(ylat eq 80.,ngot)
  if ngot eq 1 then begin
    sl=indgen(8)*9
    fd2d=reform(fdseas(*,kl,*))
    fd2dm=fd2d*!values.f_nan
    fd2dm(sl,*)=fd2d(sl,*)
    fdseas(*,kl,*)=fd2dm(*,*)
  endif
  ; keep 18 values at 75N
  kl=where(ylat eq 75.,ngot)
  if ngot eq 1 then begin
    sl=indgen(18)*4
    fd2d=reform(fdseas(*,kl,*))
    fd2dm=fd2d*!values.f_nan
    fd2dm(sl,*)=fd2d(sl,*)
    fdseas(*,kl,*)=fd2dm(*,*)
  endif
  fdseaslow=fdseas
  fdseas=fdseasfull
endif
;
if 0 eq 1 then begin
; Following discussion with Phil, certain regions are removed.
; We've already removed everything prior to 1922.
;
; (1) Remove 80N and 85N (due to Arctic High problems etc.)
;
kl=where(ylat lt 80.,ny)
ylat=ylat(kl)
fdseas=fdseas(*,kl,*)
;
; (2) Remove 15N
;
kl=where(ylat gt 15.,ny)
ylat=ylat(kl)
fdseas=fdseas(*,kl,*)
;
; (3) Set to missing the pre-1950 data at 70N & 75N from 100E to 60W
; (40W if HARSH option is selected)
;
if doharsh eq 0 then xlmax=300. else xlmax=320.
kx=where( (xlon gt 100.) and (xlon le xlmax) , nxmiss)
ky=where( ylat ge 70. , nymiss)
kt=where( timey lt 1950. , ntmiss)
for ix = 0 , nxmiss-1 do begin
  for iy = 0 , nymiss-1 do begin
    for it = 0 , ntmiss-1 do begin
      fdseas(kx(ix),ky(iy),kt(it))=!values.f_nan
    endfor
  endfor
endfor
;
; Extra data removal if HARSH option is selected
;
if doharsh eq 1 then begin
  ;
  ; (4) Set to missing the pre-1950 data at 20N from 150E to 70W (290E) and
  ; south of 30N from 150E to 130W (230E)
  ;
  kx=where( (xlon ge 150.) and (xlon le 290.) , nxmiss)
  ky=where( ylat eq 20. , nymiss)
  kt=where( timey lt 1950. , ntmiss)
  for ix = 0 , nxmiss-1 do begin
    for iy = 0 , nymiss-1 do begin
      for it = 0 , ntmiss-1 do begin
        fdseas(kx(ix),ky(iy),kt(it))=!values.f_nan
      endfor
    endfor
  endfor
  kx=where( (xlon ge 150.) and (xlon le 230.) , nxmiss)
  ky=where( ylat le 30. , nymiss)
  for ix = 0 , nxmiss-1 do begin
    for iy = 0 , nymiss-1 do begin
      for it = 0 , ntmiss-1 do begin
        fdseas(kx(ix),ky(iy),kt(it))=!values.f_nan
      endfor
    endfor
  endfor
endif
endif
  ;
  ; (5) Set to missing all data south of 25N from 10W (350E) to 110E and
  ; south of 40N from 40E to 110E
  ;
;  kx=where( (xlon ge 350.) or (xlon le 110.) , nxmiss)
  kx=where( (xlon ge -10.) and (xlon le 110.) , nxmiss)
  ky=where( ylat le 25. , nymiss)
  for ix = 0 , nxmiss-1 do begin
    for iy = 0 , nymiss-1 do begin
      fdseaslow(kx(ix),ky(iy),*)=!values.f_nan
    endfor
  endfor
  kx=where( (xlon ge 40.) and (xlon le 110.) , nxmiss)
  ky=where( ylat le 40. , nymiss)
  for ix = 0 , nxmiss-1 do begin
    for iy = 0 , nymiss-1 do begin
      fdseaslow(kx(ix),ky(iy),*)=!values.f_nan
    endfor
  endfor
  ;
;
; Now compute 1961-90 mean and s.d. (requiring 25 yrs of data as minimum)
;
kt=where( (timey ge 1961) and (timey le 1990) , nref)
fdd=double(fdseas(*,*,kt))
fdtot=total(fdd,3,/nan)
fdnum=float(total(finite(fdd),3))
fdltm=fdtot/fdnum
ml=where(fdnum lt 25.,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)
;
; Convert all values to anomalies
;
for i = 0 , nyr-1 do fdseas(*,*,i)=fdseas(*,*,i)-fdltm(*,*)
if dolowpole eq 1 then begin
  for i = 0 , nyr-1 do fdseaslow(*,*,i)=fdseaslow(*,*,i)-fdltm(*,*)
endif
;
; Also produce 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 dolowpole eq 1 then begin
  save,filename='obs_gmslp_as.idlsave',timey,fdseas,xlon,ylat,nyr,nx,ny,$
    fdltm,fdsd,missfrac,fdseaslow
endif else begin
  save,filename='obs_gmslp_as.idlsave',timey,fdseas,xlon,ylat,nyr,nx,ny,$
    fdltm,fdsd,missfrac
endelse
;
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)  &  coast.double=0  &  coast.fill=1
coast.fillcolor=50
labels=def_labels(/off)
sm=def_sm(/off)  &  sm.thresh=0.1
;
def_1color,cr,cg,cb,50,color='vlgrey'
;
levels=findgen(21)*2.+1000.
c_thick=fltarr(21)+3.
;
inter_confd,fdltm,xlon,ylat,$
  title='1961-90 mean Apr-Sep MSLP',$
  map=map,coast=coast,labels=labels,sm=sm,$
  /hi_on,/follow,levels=levels,c_thick=c_thick     ;,miss_grey='white'
;
levels=findgen(21)*0.5
c_thick=fltarr(21)+3.
;
inter_confd,fdsd,xlon,ylat,$
  title='1961-90 standard deviation Apr-Sep MSLP',$
  map=map,coast=coast,labels=labels,sm=sm,$
  /hi_on,/follow,levels=levels,c_thick=c_thick     ;,miss_grey='white'
;
end