;
;*** COMPUTES PARTIAL CORRELATIONS BETWEEN TRW AND T & P, WITH MXD HELD
;*** CONSTANT.  WHAT CAN TRW GIVE US IN ADDITION TO MXD?
;
;*** CAN ALSO COMPUTE THE MULTIPLE CORRELATION WITH TRW & MXD, BUT THIS SHOULD
;*** BE USED WITH CARE BECAUSE IT WILL ALWAYS IMPROVE OVER JUST TRW OR JUST
;*** MXD, EVEN IF IT'S JUST NOISE - ALWAYS REFER BACK TO THE PARTIAL
;*** CORRELATIONS TRW|MXD TO CHECK FOR A CONSISTENT SIGNAL.
;
; THERE IS AN ERROR IN IDL's P_CORRELATE, SO I'VE HAD TO REMOVE OTHER INFLUENCES
; BY HAND.
;
; Analyses MXD data in terms of their correlation with local monthly
; temperature and precipitation.  All data up until the end of 1984
; is used (more recent data is not used, because of differing spatial
; temporal coverage.
; ALSO computes correlations after smoothing (to identify the decline).
;
;----------------------------------------------------
;
; Get chronologys
;
domult=0        ; 0=partial r (TRW|MXD), 1=multiple r (TRW&MXD)
restore,filename='alltrw.idlsave'
trw=mxd
restore,filename='allmxd.idlsave'
;  nchron,idno,idname,location,country,tree,yrstart,yrend,statlat,statlon,$
;  mxd,fraction,timey,nyr
;
; Now read in the land precipitation dataset.
; Although there is some missing data in Mark New's precip anomalies, all
; MXD boxes have sufficient data, so do not have to use surrounding boxes.
;
print,'Reading precipitation data'
ncid=ncdf_open('/cru/u2/f055/data/obs/grid/surface/precip_new_19011995.mon.nc')
pmon=crunc_rddata(ncid,tst=[1901,0],tend=[1984,11],grid=g)
ncdf_close,ncid
;
; Now read in the land air temperature dataset.
;
print,'Reading temperature data'
ncid=ncdf_open('/cru/u2/f055/data/obs/grid/surface/crutem2v_18512001.mon.nc')
tmon=crunc_rddata(ncid,tst=[1881,0],tend=[1984,11],grid=gtemp)
ncdf_close,ncid
;
; The precipitation series will be padded with missing values at its
; beginning, to make it the same length as the temperature.  Find pad length,
; and also find section length of MXD data to extract.
;
nmon=12
npad=(gtemp.nt-g.nt)/nmon
ppad=replicate(!values.f_nan,nmon*npad)
klmxd=where((timey ge 1881) and (timey le 1984),nyrm)
timey=timey(klmxd)
;
; ALTHOUGH THE FULL PERIOD IS ALWAYS READ (TO ENSURE CONSISTENT EVALUATION
; OF WHICH TEMPERATURE BOXES DO NOT HAVE 300 OR MORE MONTHS OF NON-MISSING
; 1881-1964 DATA) THE FOLLOWING LINE ALLOWS SELECTION OF A SMALLER SUB-PERIOD
; FOR ANALYSIS.
;
klperiod=where((timey ge 1881) and (timey le 1984),nyrper)
print,'LENGTH OF ANALYSIS PERIOD (YEARS) =',nyrper,'      *************'
;
; Define arrays
;
thalf=10.    ; filter for smoothed series, lower the better cos of sample size
ncorr=21     ; 16 months (JJASOND-JFMAMJJAS) & 5 seasons (O-S O-M A-S M-A JJA)
nvar=2       ; 0=MXD,precip, 1=MXD,temperature
nyrp=g.nt/nmon
nyrt=gtemp.nt/nmon
if (nyrp+npad) ne nyrm then message,'Ooops!!!'
if nyrt ne nyrm then message,'Ooops!!!'
nyr=nyrm
z=!values.f_nan
allr=fltarr(nchron,ncorr,nvar)*z      ; correlations with P or T
lowr=fltarr(nchron,ncorr,nvar)*z      ; low freq correlation with P or T
allp=fltarr(nchron,ncorr,nvar)*z      ; partial correlations with P or T
lowp=fltarr(nchron,ncorr,nvar)*z    ; low freq partial correlations with P or T
varname=['Prec','Temp']
corrname=['J','J','A','S','O','N','D','J','F','M','A','M','J','J','A','S',$
  'Oct-Sep','Oct-Mar','Apr-Sep','May-Aug','Jun-Aug']
;
; Analyse each chronology in turn
;
runlen=0.
runnum=0.
for i = 0 , nchron-1 do begin
  print,i,format='($,I4)'
  ;
  ; Extract MXD timeseries
  ;
  trw1=reform(trw(*,i))
  trw1=trw1(klmxd)
  mxd1=reform(mxd(*,i))
  mxd1=mxd1(klmxd)
  ;
  ; Locate grid box containing site and extract monthly timeseries
  ;
  x=statlon(i)
  y=statlat(i)
  statval=[1.]
  fdout=gridit(g.nx,g.ny,g.x,g.y,x,y,statval)
  loclist=where(finite(fdout),nloc)
  if nloc ne 1 then message,'Ooops!!!'
  locx=loclist(0) mod g.nx
  locy=fix(loclist(0)/g.nx)
  ;
  ; Simple single box mean for precipitation, but pad with missing
  ;
  pre12=[ppad,reform(pmon(locx,locy,*))]
  pre12=reform(pre12,nmon,nyr)
  ;
  ; For temperature, if less than 300 months of 1881-1964 data use an
  ; average of the nine surrounding boxes (variance adjusted etc.)
  ;
  tem12=reform(tmon(locx,locy,*),nmon,nyr)
  kl=where(gtemp.year le 1964)
  dummy=where(finite(tem12(*,kl)),ngot)
  if ngot lt 300 then begin
    print
    ; extract timeseries from 9 surrounding boxes (NB. zonally cyclic)
    locy=[locy-1,locy,locy+1]
    locx=[locx-1,locx,locx+1] mod g.nx
    neglist=where(locx lt 0,nneg)
    if nneg gt 0 then locx(neglist)=locx(neglist)+g.nx
    temboxes=tmon(locx,*,*)
    temboxes=temboxes(*,locy,*)
    temboxes=reform(temboxes,9,nmon,nyr)
    ; average them and adjust variance on a month by month basis
    ; I've altered mkeffective to adjust variance to that of a single
    ; input timeseries
    for imon = 0 , nmon-1 do begin
      print,imon,format='($,I4)'
      tsin=reform(temboxes(*,imon,*))
      tsout=var_adjust(tsin,/no_anomaly,/td_rbar,/mkrbar,/mkcorr)
      tem12(imon,*)=tsout(*)
    endfor
    print
  endif
  ;
  ; Extract period for analysis
  ;
  trw1=trw1(klperiod)
  mxd1=mxd1(klperiod)
  pre12=pre12(*,klperiod)
  tem12=tem12(*,klperiod)
  ;
  ; Repeat for precipitation then for temperature
  ;
  for ivar = 0 , nvar-1 do begin
    if ivar eq 0 then begin
      var12=pre12             ; variable to analyse
      con12=tem12             ; variable to hold constant
    endif else begin
      var12=tem12
      con12=pre12
    endelse
    ;
    ; Now correlate the MXD with monthly data and seasonal data
    ;
    for icorr = 0 , ncorr-1 do begin
      ; Select month (possibly lagged) and season of instrumental data
      if icorr le 6 then begin
        imon=icorr+5
        ; Full correlation
        pre1=reform(var12(imon,*))
        pre1=shift(pre1,1)
        pre1(0)=!values.f_nan
        ; Partial correlation
        con1=reform(con12(imon,*))
        con1=shift(con1,1)
        con1(0)=!values.f_nan
      endif
      if (icorr ge 7) and (icorr le 15) then begin
        imon=icorr-7
        ; Full correlation
        pre1=reform(var12(imon,*))
        ; Partial correlation
        con1=reform(con12(imon,*))
      endif
      if icorr ge 16 then begin
        stseas=[9,9,3,4,5]
        enseas=[8,2,8,7,7]
        dtseas=[10,5,5,3,2]
        iseas=icorr-16
       pre1=mkseason(var12,stseas[iseas],enseas[iseas],datathresh=dtseas[iseas])
       con1=mkseason(con12,stseas[iseas],enseas[iseas],datathresh=dtseas[iseas])
      endif
      ; Now do the correlations if there's enough data
      kl=where(finite(mxd1+pre1+trw1),nkeep)
      if nkeep gt 5 then begin
       if domult eq 0 then begin
        ; predict pre1 and trw1 from the mxd1 data
        dummy=linfit(mxd1[kl],pre1[kl],yfit=pfit)
        dummy=linfit(mxd1[kl],trw1[kl],yfit=tfit)
        ; correlate trw1 against the residual after subtracting the prediction
        allr(i,icorr,ivar)=correlate( trw1[kl]-tfit , pre1[kl]-pfit )
        ; old approach was to use partial correlations - should be the same
        ; and indeed is so when only one influence is removed - FAILS FOR MORE
        ; THAN ONE!!!
;allr(i,icorr,ivar)=p_correlate( trw1(kl) , pre1(kl) , reform(mxd1(kl),1,nkeep) )
       endif else begin
        mpred=transpose(reform([mxd1[kl],trw1[kl]],nkeep,2))
        allr(i,icorr,ivar)=m_correlate( mpred , pre1[kl] )
       endelse
        filter_cru,thalf,/nan,tsin=mxd1,tslow=tl1
        filter_cru,thalf,/nan,tsin=pre1,tslow=tl2
        filter_cru,thalf,/nan,tsin=trw1,tslow=tl3
       if domult eq 0 then begin
        ; predict pre1 from the mxd1 data
        dummy=linfit(tl1[kl],tl2[kl],yfit=yfit2)
        dummy=linfit(tl1[kl],tl3[kl],yfit=yfit3)
        ; correlate trw1 against the residual after subtracting the prediction
        lowr(i,icorr,ivar)=correlate( tl3[kl]-yfit3 , tl2[kl]-yfit2 )
        ; old approach was to use partial correlations - should be the same
        ; and indeed is so when only one influence is removed - FAILS FOR MORE
        ; THAN ONE!!!
;lowr(i,icorr,ivar)=p_correlate( tl3(kl) , tl2(kl) , reform(tl1(kl),1,nkeep) )
       endif else begin
        mpred=transpose(reform([tl1[kl],tl3[kl]],nkeep,2))
        lowr(i,icorr,ivar)=m_correlate( mpred , tl2[kl] )
       endelse
      endif
      ; Now do the partial correlations if there's enough data
      kl=where(finite(mxd1+pre1+con1+trw1),nkeep)
      if nkeep gt 5 then begin
       if domult eq 0 then begin
        ; OLD APPROACH FAILED BECAUSE OF BUG IN P_CORRELATE WHEN >1 INFLUENCE
        ; TO REMOVE.  HENCE REMOVE MULTIPLE INFLUENCES USING REGRESS, FROM
        ; BOTH VARIABLES BEING USED!
        reminf=transpose(reform([con1[kl],mxd1[kl]],nkeep,2))
        dummy=regress(reminf,pre1[kl],yfit=pfit)
        dummy=regress(reminf,trw1[kl],yfit=tfit)
        allp(i,icorr,ivar)=correlate( trw1[kl]-tfit , pre1[kl]-pfit )
;        allp(i,icorr,ivar)=p_correlate( trw1(kl) , pre1(kl) , $
;          reform([con1(kl),mxd1(kl)],2,nkeep) )
        filter_cru,thalf,/nan,tsin=con1,tslow=tl4
        reminf=transpose(reform([tl4[kl],tl1[kl]],nkeep,2))
        dummy=regress(reminf,tl2[kl],yfit=yfit2)
        dummy=regress(reminf,tl3[kl],yfit=yfit3)
        lowp(i,icorr,ivar)=correlate( tl3[kl]-yfit3 , tl2[kl]-yfit2 )
;        lowp(i,icorr,ivar)=p_correlate( tl3(kl) , tl2(kl) , $
;          reform([tl4(kl),tl1(kl)],2,nkeep) )
       endif else begin
        ; Not done if computing multiple correlations
        allp(i,icorr,ivar)=!values.f_nan
        lowp(i,icorr,ivar)=!values.f_nan
       endelse
      endif
      ;
    endfor
    ;
    ; Repeat for next variable
    ;
  endfor
  ;
  ; Repeat for next site
  ;
endfor
print
;
runlen=runlen/runnum
;
; Save all the data
;
if domult eq 0 then fnout='trwmxd_moncorr.idlsave' $
               else fnout='trwmxdmult_moncorr.idlsave'
;
save,filename=fnout,$
  allr,ncorr,nvar,nchron,varname,corrname,statlat,statlon,$
  allp,lowr,lowp,thalf
;
end