; Computes SLP patterns to accompany PC1 of old SLP dataset (UKMO analyses,
; with considerable removal of data poor regions, rather than GMSLP data).
; Then it plots it.
;
; Computes the observed SLP, TEMP and PREC signals that go with the observed
; summer SLP modes.  The signals are the regression coefficients between SLP
; anomalies, TEMP anomalies and PREC % anomalies, and the normalised (over
; 1961-90) principal component timeseries of each mode.
; USE CORRELATIONS NOT REGRESSIONS!!!!
;
doold=1   ; 0=use EOFs from new SLP data, 1=use EOFs from old SLP data
;
; Get the summer SLP modes
;
print,'Restoring summer SLP modes'
if doold eq 0 then restore,filename='obs_summer_modes.idlsave' $
              else restore,filename='obs_summer_modes.idlsave.old2'
;  timey,nretain,useinfill,usecorr,userot,xlon,ylat,nyr,nx,ny,fillea,ev
pctime=timey
pcnyr=n_elements(pctime)
;
; Define variables to analyse
;
varname=['SLP']
nvar=n_elements(varname)
;
; Repeat for each variable separately
;
for ivar = 0 , nvar-1 do begin
  print,varname(ivar)
  ;
  ; First restore the appropriate data (assume they cover the same period!)
  ;
  print,'Restoring gridded data'
  case ivar of
    0: begin
      if doold eq 1 then restore,filename='obs_mslp_as_infill.idlsave' $
                    else restore,filename='obs_gmslp_as.idlsave'
      end
    1: begin
      restore,filename='obs_temp_as.idlsave'
      end
    2: begin
      restore,filename='obs_prec_as.idlsave'
      end
  endcase
  nyr=pcnyr
  kl=where((timey ge pctime(0)) and (timey le pctime(nyr-1)),ngot)
  timey=timey(kl)
  fdseas=fdseas(*,*,kl)
  if ngot lt nyr then begin
    timey=[timey,findgen(nyr-ngot)+timey(ngot-1)+1]
    fdseas2=fltarr(nx,ny,nyr)*!values.f_nan
    fdseas2(*,*,0:ngot-1)=fdseas(*,*,*)
    fdseas=fdseas2
  endif
  ;
  dummy=where((timey-pctime) ne 0,nerror)
  if nerror gt 0 then message,'Ooops!!!'
  fdseas=reform(fdseas,nx*ny,nyr)
  allylat=fltarr(nx,ny)
  for ix = 0 , nx-1 do allylat(ix,*)=ylat(*)
  allylat=reform(allylat,nx*ny)
  ;
  ; Define array to store results
  ;
  allarp=fltarr(nx,ny,nretain)
  ;
  ; Repeat for each mode
  ;
  for iretain = 0 , nretain-1 do begin
    print,iretain
    ;
    ; Extract timeseries and normalise it
    ;
    print,'Normalising timeseries'
    onets=reform(fillea(*,iretain))
    mknormal,onets,pctime,refperiod=[1961,1990]
    if ivar eq 0 then begin
      if iretain eq 0 then modets=fltarr(nyr,nretain)
      modets(*,iretain)=onets(*)
    endif
    ;
    print,'Computing signal pattern'
    patterndata,fdseas,onets,zcoeff=zcoeff,meanskill=meanskill,ylat=allylat,$
      zcorr=zcorr
    print,meanskill
    ;
    ; Convert to % anomalies for precip
    ;
;    case ivar of
;      0: allarp(*,*,iretain)=reform(zcorr,nx,ny)
;      1: allarp(*,*,iretain)=reform(zcoeff,nx,ny)
;      2: allarp(*,*,iretain)=100.*reform(zcoeff,nx,ny)/fdltm
;    endcase
     allarp(*,*,iretain)=reform(zcorr,nx,ny)
    ;
  endfor
  ;
  ; Keep signal patterns in appropriate variable names
  ;
  case ivar of
    0: begin
      slparp=allarp
      gslp={ xlon: xlon, ylat: ylat, nx: nx, ny: ny }
      end
    1: begin
      temparp=allarp
      gtemp={ xlon: xlon, ylat: ylat, nx: nx, ny: ny }
      end
    2: begin
      precarp=allarp
      gprec={ xlon: xlon, ylat: ylat, nx: nx, ny: ny }
      end
  endcase
  ;
endfor
;
; Plots the observed SLP, TEMP and PREC signals that go with the observed
; summer SLP modes.  The signals are the regression coefficients between SLP
; anomalies, TEMP anomalies and PREC % anomalies, and the normalised (over
; 1961-90) principal component timeseries of each mode.
; FOR SLP, IT IS THE CORRELATION NOT REGRESSION COEFF.
;
docol=0        ; 0=B&W 1=color
;
loadct,39
multi_plot,nrow=2,ncol=2,layout='large'
if !d.name eq 'X' then begin
  window,ysize=850
  !p.font=-1
endif else begin
  !p.font=0
  device,/helvetica,/bold,font_size=12
endelse
;
map=def_map(/npolar) & map.limit(0)=15.
labels=def_labels(/off)
sm=def_sm(/off) & sm.thresh=0.1
coast=def_coast(/get_device) & coast.double=0 & coast.fill=1
coast.fillcolor=10
if docol eq 0 then begin
  def_1color,10,color='palegrey'
  def_1color,20,color='mdgrey'
  def_1color,21,color='black'
  def_1color,22,color='black'
endif else begin
  def_1color,10,color='vlgreen'
  def_1color,20,color='mdgrey'
  def_1color,21,color='deepblue'
  def_1color,22,color='red'
endelse
;
nretain=1
;kretain=[1,2,6,8]
;nretain=nretain < 9
kretain=indgen(nretain)
kpve=[13,8,6,6,6,6,5,6,6]
for jretain = 0 , nretain-1 do begin
  iretain=kretain(jretain)
  titn=string(iretain+1,format='(I1)')
  tit1='  '+string(kpve(iretain),format='(I2)')+'%'
  ;
  ; Timeseries first
  ;
  pause
  yyy=reform(modets(*,iretain))
  plot,timey,yyy,$
    /xstyle,xtitle='Year',$
    ytitle='Normalised amplitude',ymargin=[10,10],$
    title='Summer mode #'+titn+tit1
  oplot,!x.crange,[0.,0.],linestyle=1
  filter_cru,10.,/nan,tsin=yyy,tslow=tslow
  oplot,timey,tslow,thick=5
  ;
  ; Now plot SLP pattern
  ;
  fd=reform(slparp(*,*,iretain))
  levs=findgen(8)*0.1+0.2
  levs=[-reverse(levs),levs]
  nlev=n_elements(levs)
  c_labels=replicate(1,nlev)
  c_thick=[replicate(2.,nlev/2),replicate(5.,nlev/2)]
  if docol eq 1 then c_thick=replicate(4.,nlev)
  c_colors=[replicate(21,nlev/2),replicate(22,nlev/2)]
;  sm.on=1
  ;
  inter_confd,fd,gslp.xlon,gslp.ylat,$
    coast=coast,map=map,labels=labels,sm=sm,$
    levels=levs,c_thick=c_thick,c_labels=c_labels,/follow,/hi_on,$
    miss_grey='white',c_colors=c_colors
  ;
  !p.ticklen=0.03
  lon_polar,map=map,[-150.,-120,-60,-30,30,60,120,150],/dotick
  ;
  normtop=convert_coord(map.limit(1),map.limit(0),/data,/to_normal)
  yr=!y.region
  yloc=normtop(1)+0.3*(yr(1)-normtop(1))
  xloc=0.5*total(!x.region)
  xyouts,xloc,yloc,/normal,align=0.5,$
    'Mode #'+titn+' vs. SLP anomalies'
;  xyouts,xloc,yr(0),/normal,align=0.5,$
;    'Correlations',charsize=!p.charsize*0.7
  ;
endfor
;
end