;
; 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 partial correlations, and correlations
; between a moisture deficit value and MXD.
; ALSO computes correlations after smoothing (to identify the decline).
; OPTIONALLY can replace precipitation with cloud cover!
;
;----------------------------------------------------
;
dotem=0         ; 0=use CRUTEM1, 1=use CRUTEM2v
if dotem eq 0 then begin
  fntem='lat_jones_18512000.mon.nc'
  savtem=''
endif else begin
  fntem='crutem2v_18512001.mon.nc'
  savtem='ct2v_'
endelse
;
docld=0         ; 0=T & precip, 1=T & cloud
trv=1           ; selects tree-ring-variable: 0=MXD 1=TRW 2=MXD-TRW
case trv of
  0: fnadd='mxd'
  1: fnadd='trw'
  2: fnadd='mxd-trw'
endcase
titadd=strupcase(fnadd)
;
; Get chronology locations
;
print,'Reading '+titadd+' data'
if trv eq 2 then begin
  restore,filename='alltrw.idlsave'
  trw=mxd
  restore,filename='allmxd.idlsave'
  mxd=mxd-trw
endif else begin
  restore,filename='all'+fnadd+'.idlsave'
  ;  nchron,idno,idname,location,country,tree,yrstart,yrend,statlat,statlon,$
  ;  mxd,fraction,timey,nyr
endelse
;
; Now read in the 1961-90 monthly means of precip and temperature
;
print,'Reading precipitation baseline'
if docld then cldfn='cld_new_19311990' else cldfn='precip_new_19611990'
restore,'/cru/u2/f055/data/obs/grid/surface/'+cldfn+'.idlsave'
;  g,nmon,ltm5
pre6190=ltm5
print,'Reading temperature baseline'
restore,'/cru/u2/f055/data/obs/grid/surface/lat_new_19611990.idlsave'
;  g,nmon,ltm5
lat6190=ltm5
;
; 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'
if docld then cldfn='cldanom' else cldfn='precip'
ncid=ncdf_open('/cru/u2/f055/data/obs/grid/surface/'+cldfn+'_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'
print,fntem
ncid=ncdf_open('/cru/u2/f055/data/obs/grid/surface/'+fntem)
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
moir=fltarr(nchron)*z                 ; corr with moisture deficit
moir2=fltarr(nchron)*z                ; corr Apr-Sep T vs moist deficit
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
  ;
  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
  ;
  ; Compute the monthly PET year by year
  ;
  ; Define PET array
  monpet=fltarr(nmon,nyr)*!values.f_nan
  ; Get absolute temperature & precip for 1961-90
  tabs=reform(lat6190(locx,locy,*))
  pabs=reform(pre6190(locx,locy,*))
  ; Repeat for each year
  for j = 0 , nyr-1 do begin
    onet=reform(tem12(*,j))+tabs(*)
    if finite(total(onet)) then begin
      onepet=thornthwaite(onet,g.y(locy))
      monpet(*,j)=onepet(*)
    endif
  endfor
  ;
  ; Compute April to July sum of PET
  ;
  j2apet=monpet(3:6,*)
  j2apet=total(j2apet,1)
  ;
  ; Compute November to July precipitation (TOTAL not MEAN)
  ;
  pre1=mkseason(pre12,10,6,datathresh=8)      ; need 8 out of 9 for mean
  pre1=pre1*9.                                ; convert to 9 month total
  ;
  ; Moisture deficit at end of July
  ;
  mdef=j2apet-pre1
  ;
  save,filename='datastore/'+savtem+'example_chron'+string(i,format='(I3.3)')+'.idlsave',$
    mxd1,pre12,tem12,nmon,nyr,timey,i,location
  ;
  ; Extract period for analysis
  ;
  mxd1=mxd1(klperiod)
  pre12=pre12(*,klperiod)
  tem12=tem12(*,klperiod)
  mdef=mdef(klperiod)
  ;
  ; Correlate moisture with MXD
  ;
  kl=where(finite(mxd1+mdef),nkeep)
  if nkeep gt 5 then begin
    filter_cru,10.,/nan,tsin=mxd1,tshigh=th1
    filter_cru,10.,/nan,tsin=mdef,tshigh=th2
;    moir(i)=correlate(mxd1(kl),mdef(kl))       ; do raw series
    moir(i)=correlate(th1(kl),th2(kl))       ; do high-pass series
  endif
  ;
  ; 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
    ; (1) With previous year's months (June to December)
    ;
    icorr=0
    for imon = 5 , nmon-1 do begin
      ; Full correlation
      pre1=reform(var12(imon,*))
      pre1=shift(pre1,1)
      pre1(0)=!values.f_nan
      kl=where(finite(mxd1+pre1),nkeep)
      if nkeep gt 5 then begin
        allr(i,icorr,ivar)=correlate(mxd1(kl),pre1(kl))
        filter_cru,thalf,/nan,tsin=mxd1,tslow=tl1
        filter_cru,thalf,/nan,tsin=pre1,tslow=tl2
        lowr(i,icorr,ivar)=correlate(tl1(kl),tl2(kl))
      endif
      ; Partial correlation
      con1=reform(con12(imon,*))
      con1=shift(con1,1)
      con1(0)=!values.f_nan
      kl=where(finite(mxd1+pre1+con1),nkeep)
      if nkeep gt 5 then begin
  allp(i,icorr,ivar)=p_correlate(pre1(kl),mxd1(kl),reform(con1(kl),1,nkeep))
        filter_cru,thalf,/nan,tsin=pre1,tslow=tl1
        filter_cru,thalf,/nan,tsin=mxd1,tslow=tl2
        filter_cru,thalf,/nan,tsin=con1,tslow=tl3
        lowp(i,icorr,ivar)=p_correlate(tl1(kl),tl2(kl),reform(tl3(kl),1,nkeep))
      endif
      ;
      icorr=icorr+1
    endfor
    ;
    ; (2) With current year's months (January to September)
    ;
    for imon = 0 , 8 do begin
      ; Full correlation
      pre1=reform(var12(imon,*))
      kl=where(finite(mxd1+pre1),nkeep)
      if nkeep gt 5 then begin
        allr(i,icorr,ivar)=correlate(mxd1(kl),pre1(kl))
        filter_cru,thalf,/nan,tsin=mxd1,tslow=tl1
        filter_cru,thalf,/nan,tsin=pre1,tslow=tl2
        lowr(i,icorr,ivar)=correlate(tl1(kl),tl2(kl))
      endif
      ; Partial correlation
      con1=reform(con12(imon,*))
      kl=where(finite(mxd1+pre1+con1),nkeep)
      if nkeep gt 5 then begin
  allp(i,icorr,ivar)=p_correlate(pre1(kl),mxd1(kl),reform(con1(kl),1,nkeep))
        filter_cru,thalf,/nan,tsin=pre1,tslow=tl1
        filter_cru,thalf,/nan,tsin=mxd1,tslow=tl2
        filter_cru,thalf,/nan,tsin=con1,tslow=tl3
        lowp(i,icorr,ivar)=p_correlate(tl1(kl),tl2(kl),reform(tl3(kl),1,nkeep))
      endif
      ;
      icorr=icorr+1
    endfor
    ;
    ; (3) With Oct to Sep annual mean and Oct to Mar winter mean
    ;
    ; Full correlation
    pre1=mkseason(var12,9,8,datathresh=10)
    kl=where(finite(mxd1+pre1),nkeep)
    if nkeep gt 5 then begin
      allr(i,icorr,ivar)=correlate(mxd1(kl),pre1(kl))
      filter_cru,thalf,/nan,tsin=mxd1,tslow=tl1
      filter_cru,thalf,/nan,tsin=pre1,tslow=tl2
      lowr(i,icorr,ivar)=correlate(tl1(kl),tl2(kl))
    endif
    ; Partial correlation
    con1=mkseason(con12,9,8,datathresh=10)
    kl=where(finite(mxd1+pre1+con1),nkeep)
    if nkeep gt 5 then begin
  allp(i,icorr,ivar)=p_correlate(pre1(kl),mxd1(kl),reform(con1(kl),1,nkeep))
      filter_cru,thalf,/nan,tsin=pre1,tslow=tl1
      filter_cru,thalf,/nan,tsin=mxd1,tslow=tl2
      filter_cru,thalf,/nan,tsin=con1,tslow=tl3
      lowp(i,icorr,ivar)=p_correlate(tl1(kl),tl2(kl),reform(tl3(kl),1,nkeep))
    endif
    ;
    icorr=icorr+1
    ;
    ; Full correlation
    pre1=mkseason(var12,9,2,datathresh=5)
    kl=where(finite(mxd1+pre1),nkeep)
    if nkeep gt 5 then begin
      allr(i,icorr,ivar)=correlate(mxd1(kl),pre1(kl))
      filter_cru,thalf,/nan,tsin=mxd1,tslow=tl1
      filter_cru,thalf,/nan,tsin=pre1,tslow=tl2
      lowr(i,icorr,ivar)=correlate(tl1(kl),tl2(kl))
    endif
    ; Partial correlation
    con1=mkseason(con12,9,2,datathresh=5)
    kl=where(finite(mxd1+pre1+con1),nkeep)
    if nkeep gt 5 then begin
  allp(i,icorr,ivar)=p_correlate(pre1(kl),mxd1(kl),reform(con1(kl),1,nkeep))
      filter_cru,thalf,/nan,tsin=pre1,tslow=tl1
      filter_cru,thalf,/nan,tsin=mxd1,tslow=tl2
      filter_cru,thalf,/nan,tsin=con1,tslow=tl3
      lowp(i,icorr,ivar)=p_correlate(tl1(kl),tl2(kl),reform(tl3(kl),1,nkeep))
    endif
    ;
    icorr=icorr+1
    ;
    ; (4) April-September, May-August, and June-August means
    ;
    ; Full correlation
    pre1=mkseason(var12,3,8,datathresh=5)
    kl=where(finite(mxd1+pre1),nkeep)
    runlen=runlen+nkeep
    runnum=runnum+1.
    if nkeep gt 5 then begin
      allr(i,icorr,ivar)=correlate(mxd1(kl),pre1(kl))
      filter_cru,thalf,/nan,tsin=mxd1,tslow=tl1
      filter_cru,thalf,/nan,tsin=pre1,tslow=tl2
      lowr(i,icorr,ivar)=correlate(tl1(kl),tl2(kl))
    endif
    ; Partial correlation
    con1=mkseason(con12,3,8,datathresh=5)
    kl=where(finite(mxd1+pre1+con1),nkeep)
    if nkeep gt 5 then begin
  allp(i,icorr,ivar)=p_correlate(pre1(kl),mxd1(kl),reform(con1(kl),1,nkeep))
      filter_cru,thalf,/nan,tsin=pre1,tslow=tl1
      filter_cru,thalf,/nan,tsin=mxd1,tslow=tl2
      filter_cru,thalf,/nan,tsin=con1,tslow=tl3
      lowp(i,icorr,ivar)=p_correlate(tl1(kl),tl2(kl),reform(tl3(kl),1,nkeep))
    endif
    ;
    icorr=icorr+1
    ;
    ; Moisture deficit at end of July and correlate with Apr-Sep temp
    ;
    if ivar eq 1 then begin
      kl=where(finite(pre1+mdef),nkeep)
      if nkeep gt 5 then moir2(i)=correlate(pre1(kl),mdef(kl))
    endif
    ;
    ; Full correlation
    pre1=mkseason(var12,4,7,datathresh=4)
    kl=where(finite(mxd1+pre1),nkeep)
    if nkeep gt 5 then begin
      allr(i,icorr,ivar)=correlate(mxd1(kl),pre1(kl))
      filter_cru,thalf,/nan,tsin=mxd1,tslow=tl1
      filter_cru,thalf,/nan,tsin=pre1,tslow=tl2
      lowr(i,icorr,ivar)=correlate(tl1(kl),tl2(kl))
    endif
    ; Partial correlation
    con1=mkseason(con12,4,7,datathresh=4)
    kl=where(finite(mxd1+pre1+con1),nkeep)
    if nkeep gt 5 then begin
  allp(i,icorr,ivar)=p_correlate(pre1(kl),mxd1(kl),reform(con1(kl),1,nkeep))
      filter_cru,thalf,/nan,tsin=pre1,tslow=tl1
      filter_cru,thalf,/nan,tsin=mxd1,tslow=tl2
      filter_cru,thalf,/nan,tsin=con1,tslow=tl3
      lowp(i,icorr,ivar)=p_correlate(tl1(kl),tl2(kl),reform(tl3(kl),1,nkeep))
    endif
    ;
    icorr=icorr+1
    ;
    ; Full correlation
    pre1=mkseason(var12,5,7,datathresh=2)      ; NEWEST!!! USE JUNE-AUGUST
;    pre1=mkseason(var12,6,7,datathresh=2)      ; NEW!!! USE JULY-AUGUST
;    pre1=mkseason(var12,5,6,datathresh=2)     ; OLD!!! USE JUNE-JULY
    kl=where(finite(mxd1+pre1),nkeep)
    if nkeep gt 5 then begin
      allr(i,icorr,ivar)=correlate(mxd1(kl),pre1(kl))
      filter_cru,thalf,/nan,tsin=mxd1,tslow=tl1
      filter_cru,thalf,/nan,tsin=pre1,tslow=tl2
      lowr(i,icorr,ivar)=correlate(tl1(kl),tl2(kl))
    endif
    ; Partial correlation
    con1=mkseason(con12,5,7,datathresh=2)
    kl=where(finite(mxd1+pre1+con1),nkeep)
    if nkeep gt 5 then begin
  allp(i,icorr,ivar)=p_correlate(pre1(kl),mxd1(kl),reform(con1(kl),1,nkeep))
      filter_cru,thalf,/nan,tsin=pre1,tslow=tl1
      filter_cru,thalf,/nan,tsin=mxd1,tslow=tl2
      filter_cru,thalf,/nan,tsin=con1,tslow=tl3
      lowp(i,icorr,ivar)=p_correlate(tl1(kl),tl2(kl),reform(tl3(kl),1,nkeep))
    endif
    ;
    icorr=icorr+1
    ;
    ; Repeat for next variable
    ;
  endfor
  ;
  ; Repeat for next site
  ;
endfor
print
;
runlen=runlen/runnum
print,'Average Apr-Sep vs '+titadd+' correlation sample length=',runlen
;
; Save all the data
;
if docld eq 0 then fnsave='datastore/'+savtem+fnadd+'_moncorr.idlsave' $
              else fnsave='datastore/'+savtem+fnadd+'_moncorr_cldtem.idlsave'
print,fnsave
save,filename=fnsave,$
  allr,ncorr,nvar,nchron,varname,corrname,statlat,statlon,$
  allp,moir,moir2,lowr,lowp,thalf
;
end