;
; 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=0   ; 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','Temperature','Precipiation']
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
;      restore,filename='obs_mslp_as_infill.idlsave'
      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
;
save,filename='obs_stp_modes.idlsave',$
  slparp,temparp,precarp,modets,timey,nyr,nvar,varname,$
  gslp,gtemp,gprec,nretain
;
end