;
; Computes EOFs of observed Apr-Sep SLP, to test and compare different methods.
; (1) Compare infilled and uninfilled datasets
; (2) Compare non-rotated covariance & correlation, and rotated covariance
;     and correlation EOFs.
;
usegmslp=2      ; 0=UKMO analyses, 1=GMSLP, 2=GMSLP with lower polar resolution
useinfill=1     ; 0=un-infilled, 1=infilled (not necessary for GMSLP)
usecorr=1       ; 0=covariance,  1=correlation
userot=1        ; 0=unrotated,   1=rotated
nretain=8
compindex=1     ; 0=do not, 1=do compare with other circulation indices
;
; Restore NOAA mode summer timeseries
;
if compindex eq 1 then begin
  restore,filename='noaa_summer.idlsave'
  ;  nsummer,summerind,summerts,summertslow
  restore,filename='nao_summer.idlsave'
  ;  timey,seasts,tslow
  naosummer=seasts
  naoshort=naosummer(where(timey ge 1950))
  ; Select series to combine and overplot the PCs
  mnan=!values.f_nan
  iuse=intarr(nsummer+1,9)
  iuse(*,0)=[mnan,mnan,mnan,mnan,mnan,mnan,mnan,mnan,mnan,mnan,mnan]
  iuse(*,1)=[mnan,mnan,mnan,mnan,+0.5,-1.0,mnan,mnan,mnan,mnan,mnan]
  iuse(*,2)=[mnan,mnan,mnan,mnan,mnan,mnan,mnan,-1.0,mnan,mnan,+1.0]
  iuse(*,3)=[+1.0,mnan,mnan,mnan,mnan,mnan,-0.5,mnan,mnan,mnan,mnan]
  iuse(*,4)=[-0.5,-1.0,mnan,mnan,mnan,mnan,mnan,mnan,mnan,mnan,mnan]
  iuse(*,5)=[mnan,mnan,mnan,-1.0,mnan,mnan,mnan,mnan,mnan,mnan,mnan]
  iuse(*,6)=[+0.5,mnan,+1.0,mnan,mnan,mnan,-0.5,mnan,mnan,mnan,mnan]
  iuse(*,7)=[mnan,mnan,mnan,mnan,mnan,mnan,mnan,mnan,mnan,mnan,mnan]
  iuse(*,8)=[mnan,mnan,+1.0,mnan,mnan,mnan,mnan,mnan,mnan,mnan,+0.5]
endif
;
; Restore Apr-Sep MSLP gridded dataset
;
if usegmslp eq 0 then begin
  if useinfill eq 0 then restore,filename='obs_mslp_as.idlsave' $
                    else restore,filename='obs_mslp_as_infill.idlsave'
endif else begin
  restore,filename='obs_gmslp_as.idlsave'
  if usegmslp eq 2 then fdseas=fdseaslow
endelse
;  timey,fdseas,xlon,ylat,nyr,nx,ny,fdltm,fdsd,missfrac
;
; Select the years to analyse
;
ksub=where(missfrac le 0.19,nsub)
if (usegmslp eq 0) and (nsub ne 73) then message,'Ooops!'
print,'Analysing years:',nsub
;
fdsub=fdseas(*,*,ksub)
yrsub=timey(ksub)
;
; Compute them
;
print,'Computing EOFs'
myeof2d_rot,fdsub,ev,ea,lam,lampct,lamcum,$
  evrot,earot,lamrot,lampctrot,lamcumrot,$
  nretain=nretain,correlation=usecorr
;myeof2d,fdkeep,ev,ea,lam,lampct,lamcum,$
;  nretain=nretain,correlation=cormat
;
; Choose rotated or unrotated EOFs
;
if userot eq 1 then begin
  ev=evrot & ea=earot & lam=lamrot & lampct=lampctrot & lamcum=lamcumrot
  ; Rotated ones need to be sorted into order
  isort=reverse(sort(lampct))
  lam=lam(isort)
  lampct=lampct(isort)
  lamcum(0)=lampct(0)
  for iret = 1 , nretain-1 do lamcum(iret)=lamcum(iret-1)+lampct(iret)
  ev=ev(*,*,isort)
  ea=ea(*,isort)
endif
;
; Try to infill the amplitude timeseries for those years not used in the EOF
; analysis (indeed, apply it to all the data to see if it gives the correct
; values!).  Not necessary for GMSLP which has full data.
;
if usegmslp eq 0 then begin
; First remove long-term mean (not 61-90 mean!)
fdanom=reform(fdseas,nx*ny,nyr)
mkanomaly,fdanom
fdanom=reform(fdanom,nx,ny,nyr)
allea=fltarr(nyr,nretain)
for iretain = 0 , nretain-1 do begin
  for iyr = 0 , nyr-1 do begin
    allea(iyr,iretain)=total(fdanom(*,*,iyr)*ev(*,*,iretain),/nan)
  endfor
endfor
fillea=allea
fillea(ksub,*)=ea(*,*)
endif else begin
  fillea=ea
  allea=ea
endelse
;
; Now plot them
;
loadct,39
multi_plot,nrow=3,ncol=2,layout='large'
if !d.name eq 'X' then window,ysize=850
;
lampct=[!values.f_nan,lampct]
lamcum=[0,lamcum]
xt=timey
plot,lampct,/xstyle,xtitle='Mode',ytitle='% variance explained',$
  psym=10,xrange=[0,20],/ylog
plot,lampct,/xstyle,xtitle='Mode',ytitle='% variance explained',$
  psym=10,xrange=[0,20]
plot,lamcum,/xstyle,xtitle='Mode',ytitle='Cumulative % variance explained',$
  psym=10,xrange=[0,20]
plot,[0,1],/nodata,xstyle=5,ystyle=5
xyouts,0,0.9,'PCA of observed Apr-Sep SLP'
if usegmslp eq 0 then begin
  if useinfill eq 0 then ttt='Not infilled UKMO' else ttt='Infilled UKMO'
endif else begin
  ttt='GMSLP'
endelse
xyouts,0,0.6,ttt
if usecorr eq 0 then ttt='Covariance matrix' else ttt='Correlation matrix'
xyouts,0,0.3,ttt
if userot eq 0 then ttt='Not rotated' else ttt='Varimax rotated'
xyouts,0,0.,ttt
;
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.5         ; &  sm.method=3  & sm.npole=0
;
def_1color,cr,cg,cb,50,color='lgrey'
def_1color,cr,cg,cb,100,color='mdgrey'
;
levels=findgen(25)*0.02-0.24
c_thick=[replicate(2.,12),replicate(4.,13)]
th=10.
;
for i = 0 , nretain-1 do begin
  pause
  ;
  yt=fltarr(nyr)*!values.f_nan
  yt(ksub)=ea(*,i)
  filter_cru,th,tsin=yt,tslow=tslow,/nan
  plot,xt,yt,title='Mode'+string(i+1,format='(I2)'),$
    /xstyle,$
    xtitle='Year',ytitle='Amplitude'
  ;
  ml=where(finite(yt) eq 0,nmiss)
  if nmiss gt 0 then begin
    yt(ml)=allea(ml,i)
    oplot,xt,yt,linestyle=1
  endif
  filter_cru,th,tsin=yt,tslow=tslow,/nan
  if nmiss gt 0 then begin
    oplot,xt,tslow,thick=5,linestyle=1
    tslow(ml)=!values.f_nan
  endif
  oplot,xt,tslow,thick=5
  ;
  oplot,!x.crange,[0,0],linestyle=1
  ;
  ; Compute and print correlations with NOAA summer indices
  ;
if compindex eq 1 then begin
  kcorr=where((xt ge 1950) and (xt le 1997))
  pcshort=yt(kcorr)
  allr=fltarr(nsummer)
  for isummer = 0 , nsummer-1 do $
    allr(isummer)=correlate(yt(kcorr),summerts(isummer,*))
  print,i,allr,format='(I2,10F6.2)'
  ;
  ; ...and summer stations-based NAO
  ;
  print,i,correlate(yt,naosummer)
  ;
  ; Overlay combination of NOAA and NAO indices
  ;
  mknormal,yt,refmean=refmean,refsd=refsd
  ovx=xt(kcorr)
  nov=n_elements(kcorr)
  ovy=fltarr(nov)
  for iov = 0 , nov-1 do begin
    ovy(iov)=total(iuse(*,i)*[naoshort(iov),reform(summerts(*,iov))],/nan)
  endfor
  if (max(ovy,/nan)-min(ovy,/nan)) gt 0 then begin
    mknormal,ovy
    ovy=ovy*refsd(0)+refmean(0)
    filter_cru,th,tsin=ovy,tslow=tslow,/nan
    oplot,ovx,tslow,color=100,thick=6
  endif
endif
  ;
  fd=reform(ev(*,*,i))
  inter_confd,fd,xlon,ylat,$
    map=map,coast=coast,labels=labels,sm=sm,$
    /hi_on,/follow,levels=levels,c_thick=c_thick   ;,miss_grey='white'
endfor
;
save,filename='obs_summer_modes.idlsave',$
  timey,nretain,useinfill,usecorr,userot,xlon,ylat,nyr,nx,ny,fillea,ev,$
  lampct
;
end