;
; Computes EOFs of observed (partly infilled) monthly surface temperature.
; (NH only).
;
usecoloc=0      ; 0=full data, 1=co-located with MXD sites only -1=fullregions
usefilt=0       ; 0=unfiltered,  1=pre-filtered
thalf=40.       ; cutoff for high-pass filter (this is in months not years!)
nretain=15
;
; Restore Apr-Sep land air temperature gridded dataset
;
print,'Reading in temperature data'
restore,filename='obs_temp_mon.idlsave'
;  timey,fdmon,xlon,ylat,nyr,nx,ny,fdltm,fdsd,missfrac
;
; Restore regional breakdown of grid boxes
;
print,'Reading tree and regional locations'
restore,filename='../reg_mxdboxes.idlsave'
;  nreg,boxlists,boxnumbs,regname,g,boxregs,boxprec,boxtemp
;
; Compute required mask to apply to data
;
if usecoloc ne 0 then begin
  print,'Masking out unused data'
  if usecoloc gt 0 then fdmask=total(finite(boxlists),3) $
                   else fdmask=total(finite(boxregs),3)
  kl=where(fdmask gt 0,nbox)
  print,'Retaining maximum of boxes = ',nbox
  fdmask=fltarr(nx,ny)*!values.f_nan
  fdmask(kl)=1.
  for iyr = 0 , nyr-1 do begin
    for imon = 0 , nmon-1 do begin
      fdmon(*,*,imon,iyr)=fdmon(*,*,imon,iyr)*fdmask(*,*)
    endfor
  endfor
endif
;
; Optionally pre-filter (high-pass) each grid-box time series
;
if usefilt eq 1 then begin
  print,'Filtering data: rows=',nx
  for i = 0 , nx-1 do begin 
    print,i,format='($,I4)'
    for j = 0 , ny-1 do begin
      onets=reform(fdmon(i,j,*,*),nmon*nyr)
      dummy=where(finite(onets),ngot)
      if ngot gt 5. then begin
        filter_cru,thalf,/nan,tsin=onets,tshigh=onehi
        fdmon(i,j,*,*)=reform(onehi(*),nmon,nyr)
      endif
    endfor
  endfor
  print
endif
;
; Select the years to analyse
;
ksub1=where((timey ge 1906) and (timey le 1913),nsub1)
ksub2=where((timey ge 1917) and (timey le 1993),nsub2)
ksub=[ksub1,ksub2]
nsub=nsub1+nsub2
if nsub ne 85 then message,'Ooops!'
;
fdsub=reform(fdmon(*,*,*,ksub),nx*ny,nmon*nsub)
mknormal,fdsub,refmean=pcapermn,refsd=pcapersd
fdsub=reform(fdsub,nx,ny,nmon*nsub)
pcapermn=reform(pcapermn,nx,ny)
pcapersd=reform(pcapersd,nx,ny)
yrsub=timey(ksub)
;
; Compute them
;
print,'Computing EOFs'
myeof2d,fdsub,ev,ea,lam,lampct,lamcum,$
  nretain=nretain,correlation=0    ; 0 cos we have already normalised!
ea=reform(ea,nmon,nsub,nretain)
;
; 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!).
;
; First remove long-term mean (not 61-90 mean!) of the PCA analysis period
print,'Infilling time series gaps'
print,'Making normalised anomalies'
fdanom=fdmon
for iyr = 0 , nyr-1 do begin
  print,timey(iyr),format='($,I4)'
  for imon = 0 , nmon-1 do begin
    fdanom(*,*,imon,iyr)=(fdanom(*,*,imon,iyr)-pcapermn(*,*))/pcapersd(*,*)
  endfor
endfor
print
print,'Projecting data onto EOFs',nretain
fdanom=reform(fdanom,nx,ny,nmon*nyr)
allea=fltarr(nmon*nyr,nretain)
for iretain = 0 , nretain-1 do begin
  print,iretain,format='($,I4)'
  for i = 0 , (nmon*nyr)-1 do begin
    allea(i,iretain)=total(fdanom(*,*,i)*ev(*,*,iretain),/nan)
  endfor
endfor
print
print,'Combining PCs with infilled PCs'
allea=reform(allea,nmon,nyr,nretain)
fillea=allea
fillea(*,ksub,*)=ea(*,*,*)
;
; 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 monthly temperature'
case usecoloc of
  -1: ttt='Full MXD regions'
  0: ttt='Northern Hemisphere'
  1: ttt='Co-located with MXD'
endcase
xyouts,0,0.6,ttt
if usefilt eq 1 then ttt='Pre-filtered ( < '+string(thalf,format='(I3)')+' months )' $
                else ttt='Unfiltered'
xyouts,0,0.3,ttt
xyouts,0,0.,'Correlation matrix, unrotated'
;
map=def_map(/npolar)  &  map.limit(0)=0.
coast=def_coast(/get_device)
labels=def_labels(/off)
;
levels=findgen(21)*0.02-0.2
levels(0)=-0.3
levels(20)=0.3
c_colors=indgen(20)+50
def_1color,cr,cg,cb,50,color='dblue'
def_1color,cr,cg,cb,59,color='lgreen'
def_1color,cr,cg,cb,60,color='lyellow'
def_1color,cr,cg,cb,69,color='dred'
def_smearcolor,fromto=[50,59]
def_smearcolor,fromto=[60,69]
;
th=10.
;
; Compute annual means of the monthly PCs
;
ea=total(ea,1)/12.
allea=total(allea,1)/12.
fillea=total(fillea,1)/12.
;
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)
  yt(ml)=allea(ml,i)
  oplot,xt,yt,linestyle=1
  filter_cru,th,tsin=yt,tslow=tslow,/nan
  oplot,xt,tslow,thick=5,linestyle=1
  tslow(ml)=!values.f_nan
  oplot,xt,tslow,thick=5
  ;
  oplot,!x.crange,[0,0],linestyle=1
  ;
  fd=reform(ev(*,*,i))
  inter_boxfd,fd,xlon,ylat,$
    map=map,coast=coast,labels=labels,$
    levels=levels,c_colors=c_colors
endfor
;
save,filename='obst_mon_modes.idlsave',$
  timey,nretain,usefilt,usecoloc,thalf,xlon,ylat,nyr,nx,ny,fillea,ev
;
end