;
; Computes correlation coefficients between monthly
; temperature and precipitation over all boxes with MXD data in
; them.  Analysis period is 1901-1984.
;
;----------------------------------------------------
;
dotem=1         ; 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
;
; Get chronology locations
;
print,'Reading MXD data'
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 (same period)
;
print,'Reading temperature data'
print,fntem
ncid=ncdf_open('/cru/u2/f055/data/obs/grid/surface/'+fntem)
tmon=crunc_rddata(ncid,tst=[1901,0],tend=[1984,11],grid=gtemp)
ncdf_close,ncid
;
; Define arrays
;
nmon=12
nyr=g.nt/nmon
ncorr=nmon+5      ; 5 seasons, plus individual months
allr=fltarr(nchron,ncorr)
corrname=['J','F','M','A','M','J','J','A','S','O','N','D',$
  'Oct-Sep','Oct-Mar','Apr-Sep','May-Aug','Jun-Jul']
;
; Analyse each chronology grid box (some duplication!)
;
for i = 0 , nchron-1 do begin
  print,i,format='($,I4)'
  ;
  ; 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
  ;
  pre12=reform(pmon(locx,locy,*),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(g.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
  ;
  ; Now correlate T & P month by month
  ;
  icorr=0
  for imon = 0 , nmon-1 do begin
    pre1=reform(pre12(imon,*))
    tem1=reform(tem12(imon,*))
    kl=where(finite(tem1+pre1),nkeep)
    if nkeep gt 5 then allr(i,icorr)=correlate(tem1(kl),pre1(kl))
    icorr=icorr+1
  endfor
  ;
  ; From October to September annual mean and October to March winter mean
  ;
  preshift=reform(pre12,nmon*nyr)
  preshift=shift(preshift,3)         ; shifts October value to start of year
  preshift=reform(preshift,nmon,nyr)
  preshift(*,0)=!values.f_nan
  ;
  temshift=reform(tem12,nmon*nyr)
  temshift=shift(temshift,3)         ; shifts October value to start of year
  temshift=reform(temshift,nmon,nyr)
  temshift(*,0)=!values.f_nan
  ;
  totval=total(preshift,1,/nan)
  totnum=total(finite(preshift),1)
  ml=where(totnum lt 10,nmiss)       ; need 10 out of 12 months
  pre1=totval/float(totnum)
  if nmiss gt 0 then pre1(ml)=!values.f_nan
  ;
  totval=total(temshift,1,/nan)
  totnum=total(finite(temshift),1)
  ml=where(totnum lt 10,nmiss)       ; need 10 out of 12 months
  tem1=totval/float(totnum)
  if nmiss gt 0 then tem1(ml)=!values.f_nan
  ;
  kl=where(finite(tem1+pre1),nkeep)
  if nkeep gt 5 then allr(i,icorr)=correlate(tem1(kl),pre1(kl))
  icorr=icorr+1
  ;
  preshift=preshift(0:5,*)
  totval=total(preshift,1,/nan)
  totnum=total(finite(preshift),1)
  ml=where(totnum lt 5,nmiss)       ; need 5 out of 6 months
  pre1=totval/float(totnum)
  if nmiss gt 0 then pre1(ml)=!values.f_nan
  ;
  temshift=temshift(0:5,*)
  totval=total(temshift,1,/nan)
  totnum=total(finite(temshift),1)
  ml=where(totnum lt 5,nmiss)       ; need 5 out of 6 months
  tem1=totval/float(totnum)
  if nmiss gt 0 then tem1(ml)=!values.f_nan
  ;
  kl=where(finite(tem1+pre1),nkeep)
  if nkeep gt 5 then allr(i,icorr)=correlate(tem1(kl),pre1(kl))
  icorr=icorr+1
  ;
  ; (4) April-September, May-August, and June-July means
  ;
  preshift=pre12(3:8,*)
  totval=total(preshift,1,/nan)
  totnum=total(finite(preshift),1)
  ml=where(totnum lt 5,nmiss)       ; need 5 out of 6 months
  pre1=totval/float(totnum)
  if nmiss gt 0 then pre1(ml)=!values.f_nan
  ;
  temshift=tem12(3:8,*)
  totval=total(temshift,1,/nan)
  totnum=total(finite(temshift),1)
  ml=where(totnum lt 5,nmiss)       ; need 5 out of 6 months
  tem1=totval/float(totnum)
  if nmiss gt 0 then tem1(ml)=!values.f_nan
  ;
  kl=where(finite(tem1+pre1),nkeep)
  if nkeep gt 5 then allr(i,icorr)=correlate(tem1(kl),pre1(kl))
  icorr=icorr+1
  ;
  preshift=pre12(4:7,*)
  totval=total(preshift,1,/nan)
  totnum=total(finite(preshift),1)
  ml=where(totnum lt 4,nmiss)       ; need 4 out of 4 months
  pre1=totval/float(totnum)
  if nmiss gt 0 then pre1(ml)=!values.f_nan
  ;
  temshift=tem12(4:7,*)
  totval=total(temshift,1,/nan)
  totnum=total(finite(temshift),1)
  ml=where(totnum lt 4,nmiss)       ; need 4 out of 4 months
  tem1=totval/float(totnum)
  if nmiss gt 0 then tem1(ml)=!values.f_nan
  ;
  kl=where(finite(tem1+pre1),nkeep)
  if nkeep gt 5 then allr(i,icorr)=correlate(tem1(kl),pre1(kl))
  icorr=icorr+1
  ;
  preshift=pre12(5:6,*)
  totval=total(preshift,1,/nan)
  totnum=total(finite(preshift),1)
  ml=where(totnum lt 2,nmiss)       ; need 2 out of 2 months
  pre1=totval/float(totnum)
  if nmiss gt 0 then pre1(ml)=!values.f_nan
  ;
  temshift=tem12(5:6,*)
  totval=total(temshift,1,/nan)
  totnum=total(finite(temshift),1)
  ml=where(totnum lt 2,nmiss)       ; need 2 out of 2 months
  tem1=totval/float(totnum)
  if nmiss gt 0 then tem1(ml)=!values.f_nan
  ;
  kl=where(finite(tem1+pre1),nkeep)
  if nkeep gt 5 then allr(i,icorr)=correlate(tem1(kl),pre1(kl))
  icorr=icorr+1
  ;
  ; Repeat for next site
  ;
endfor
print
;
; Save all the data
;
save,filename='datastore/'+savtem+'tp_moncorr.idlsave',$
  allr,ncorr,nchron,corrname,statlat,statlon
;
end