;
; Computes the observed SLP, TEMP and PREC signals that go with the MXD
; 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!!!!
;
; Get the MXD modes
;
print,'Restoring MXD modes'
;
restore,filename='mxdmodes_16971976unrcoradjraw.idlsave'
;  mxdyear,nretain,usefilt,thalf,g,mxdnyr,ea,ev,lampct,usecorr,$
;  usedecline,useper,userot
pctime=mxdyear
;
; Define variables to analyse
;
varname=['SLP','Temperature','Precipiation']
nvar=n_elements(varname)
;
; Repeat for each variable separately
;
nvar=2
instrts=fltarr(nvar,nretain,mxdnyr)
instrts(*,*,*)=!values.f_nan
for ivar = 0 , nvar-1 do begin
  print,varname(ivar)
  ;
  ; First restore the appropriate data
  ;
  print,'Restoring gridded data'
  case ivar of
    0: begin
      restore,filename='obs_mslp_as_infill.idlsave'
      dthresh=0.36     ; 0.36*(1976-1992+1)=20 yrs needed
      end
    1: begin
      restore,filename='obs_temp_as.idlsave'
      dthresh=0.17     ; 0.17*(1976-1856+1)=20 yrs needed
      end
    2: begin
      restore,filename='obs_prec_as.idlsave'
      end
  endcase
  ; extract overlap period
  minyr=timey(0)
  maxyr=pctime(mxdnyr-1)
  kmxd=where(pctime ge minyr,n1)
  kins=where(timey le maxyr,n2)
  if n1 ne n2 then message,'oops!'
  dummy=where((timey(kins)-pctime(kmxd)) ne 0,nerror)
  if nerror gt 0 then message,'Ooops!!!'
  ;
  fdseas=reform(fdseas,nx*ny,nyr)
  fdseas=fdseas(*,kins)
  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(ea(*,iretain))
;***TEMPORARY REPLACEMENT OF TIME SERIES BY RANDOM NOISE!
;   nele=n_elements(onets)
;   onets=randomn(seed,nele)
;   for iele = 1 , nele-1 do onets(iele)=onets(iele)+0.35*onets(iele-1)
;***END
    mknormal,onets,pctime,refperiod=[1881,1960]
    if ivar eq 0 then begin
      if iretain eq 0 then modets=fltarr(mxdnyr,nretain)
      modets(*,iretain)=onets(*)
    endif
    ;
    ; Leading mode is contaminated by decline, so pre-filter it (but not
    ; the gridded datasets!)
    ;
;    if iretain eq 0 then begin
;      filter_cru,40.,/nan,tsin=onets,tshigh=tshi
;      onets=tshi
;    endif
    ;
    print,'Computing signal pattern'
    patterndata,fdseas,onets(kmxd),zcoeff=zcoeff,meanskill=meanskill,$
      ylat=allylat,zcorr=zcorr,thresh=dthresh,zeof=zeof
    print,meanskill
    ;
    ; Project the instrumental data onto the instrumental pattern to estimate
    ; the time series associated with the pattern.  Cannot use projectdata.pro
    ; since we have missing data.
    ;
    ; First remove long-term mean (not 61-90 mean!)
    fd1=fdseas
    mkanomaly,fd1
    allea=fltarr(n1)
    for iyr = 0 , n1-1 do begin
      allea(iyr)=total(fd1(*,iyr)*zeof(*),/nan)
    endfor
    instrts(ivar,iretain,kmxd)=allea(*)
    ;
    ; 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
;
timey=pctime
nyr=mxdnyr
save,filename='mxd_stp_modes.idlsave',$
  slparp,temparp,modets,timey,nyr,nvar,varname,$
  gslp,gtemp,nretain,lampct,instrts       ; precarp,gprec
;
end