;
; Reads in some long series (single tree chronology, regional tree
; timeseries, or palaeoaverage), and compute the distribution of pre-20th
; century means and trends of various lengths from the series, and finds
; the 95% threshold of trend and mean magnitudes.  Also computes HadCM2 control
; run equivalents.
;
; Read in HadCM2 control run series for equivalent points/areas & seasons
;
restore,filename='hcci_ts.idlsave'
;  Gets:  hccits,timey,nyr,nts
modnyr=nyr
modtimey=timey
;
; Lengths of trends to consider
;
tlen=[6,8,10,15,20,25,30,35,40,45,50,55,60,65,70,75,80,85,90,95,100]
ntry=n_elements(tlen)
maxsize=min(tlen)
;
; Number of series to consider
;
nts=9
alltit=strarr(nts)
binsc=[1.,0.4,0.4,1.9,1.9,1.9,1.9,1.9,1.9]
;
; Define (open-ended) bins to use when generating distributions
; Defined by their centres.
;
nbin=41
bmhalf=0.025    ; half width of mean bins
bthalf=0.05     ; half width of trend bins
binmean=(findgen(nbin)*bmhalf*2.) - ((nbin-1)*bmhalf)
bintren=(findgen(nbin)*bthalf*2.) - ((nbin-1)*bthalf)
freqmean=fltarr(nbin,ntry,nts,2)     ; ,,,0=obs ,,,1=hcci
freqtren=fltarr(nbin,ntry,nts,2)
tren95=fltarr(ntry,nts,2)
mean95=fltarr(ntry,nts,2)
;
; Now get the timeseries to use
;
for i = 0 , nts-1 do begin
  case i of
    0: begin        ; Phil's reconstruction
      alltit(i)="Jones' Northern Hemisphere temperature reconstruction"
      ; Period to consider
      perst=1000      ; was 1400
      peren=1991      ; was 1980
      openr,1,'phil_nhrecon.dat'
      nyr=992
      rawdat=fltarr(4,nyr)
      readf,1,rawdat,format='(I5,F7.2,I3,F7.2)'
      close,1
      timey=reform(rawdat(0,*))
      ts=reform(rawdat(3,*))
      kl=where((timey ge perst) and (timey le peren),nyr)
      timey=timey(kl)
      ts=ts(kl)
      yr=[0.,0.9]
      ; Convert from normalised values to deg C relative to 1961-1990
      ts=ts*0.521 - 0.1134
      openw,lun,/get_lun,'jones1998_nh10.dat'
      printf,lun,'Jones et al. (1998) NH10 reconstruction'
      printf,lun,'Northern Hemisphere June-August mean temperature estimates'
      printf,lun,'degrees C w.r.t. 1961-1990 mean temperature'
      for iii = 0 , nyr-1 do printf,lun,timey(iii),ts(iii),format='(I6,F8.2)'
      free_lun,lun
    end
    1: begin        ; Mike Mann's reconstruction (now use his 1000 yr one)
      alltit(i)="Mann's Northern Hemisphere temperature reconstruction"
      ; Read in Jan-Dec NH reconstruction from Mann et al. 1000 yr version!
      perst=1000
      peren=1980
      openr,1,'/cru/u2/f055/data/paleo/mann9899/mann_nh1000.dat'
      nyr1=1980-1000+1
      rawdat=fltarr(2,nyr1)
      readf,1,rawdat
      close,1
      timey=reform(rawdat(0,*))
      ts=reform(rawdat(1,*))       ; -0.12 to convert to oC wrt 6190 vs. NH all
;      perst=1400
;      peren=1980
;      openr,1,'mann_nhrecon.dat'
;      nyr=596
;      headdat=' '
;      rawdat=fltarr(2,nyr)
;      readf,1,headdat
;      readf,1,rawdat,format='(I6,F11.7)'
;      close,1
;      timey=reform(rawdat(0,*))
;      ts=reform(rawdat(1,*))
      kl=where((timey ge perst) and (timey le peren),nyr)
      timey=timey(kl)
      ts=ts(kl)
      yr=[0.,0.5]
    end
    2: begin        ; NHD1     ***** REPLACE WITH HARRY'S CURVE!!!
      alltit(i)="Age-banded density NH growing-season reconstruction"
      ; Period to consider
      perst=1400
      peren=1960
;      restore,filename='../treeharry/densadj_all(330).idlsave'
;      timey=x
;          ; CONVERSION FACTORS FOR AGE-BANDED MXD, BY REGRESSION ON INSTR.
;      ts=densall*0.151434    ; converts it from density to temperature anom
; NEW STUFF!
      restore,filename='../tree6/bandtemp_calibrated.idlsave'
      ts=calregts(*,9)
;
      kl=where((timey ge perst) and (timey le peren),nyr)
      timey=timey(kl)
      ts=ts(kl)
;      ts=ts(kl)-0.139032    ; to convert it oC wrt 1961-90
      yr=[0.,0.5]
    end
    3: begin        ; Jasper
      alltit(i)="Jasper temperature reconstruction"
      openr,1,'jasper.dat'
      readf,1,nyr,format='(1X,I4)'
      rawdat=fltarr(2,nyr)
      readf,1,rawdat,format='(1X,I5,F8.2)'
      close,1
      timey=reform(rawdat(0,*))
      ts=reform(rawdat(1,*))
      yr=[0.,1.5]
    end
    4: begin        ; Polar Urals
      alltit(i)="Polar Urals temperature reconstruction"
      openr,1,'polarural.dat'
      readf,1,nyr,format='(1X,I4)'
      rawdat=fltarr(2,nyr)
      readf,1,rawdat,format='(1X,I5,F8.2)'
      close,1
      timey=reform(rawdat(0,*))
      ts=reform(rawdat(1,*))
      yr=[0.,2.6]
    end
    5: begin        ; Tornetrask
      alltit(i)="Tornetrask temperature reconstruction"
      openr,1,'tornetrask.dat'
      readf,1,nyr,format='(I4)'
      rawdat=fltarr(2,nyr)
      readf,1,rawdat,format='(I5,F8.2)'
      close,1
      timey=reform(rawdat(0,*))
      ts=reform(rawdat(1,*))
      yr=[0.,1.5]
    end
    6: begin        ; Annual CET
      alltit(i)="Annual Central England Temperature"
      openr,1,'cet.dat'
      readf,1,nyr
      rawdat=fltarr(18,nyr)
      readf,1,rawdat,format='(I4,12F5.1,5F6.2)'
      close,1
      timey=reform(rawdat(0,*))
      moncet=rawdat(1:12,*)
      ts=mkseason(moncet,0,11)
      yr=[0.,1.5]
    end
    7: begin        ; Summer CET
      alltit(i)="Summer Central England Temperature"
      timey=reform(rawdat(0,*))
      ts=mkseason(moncet,3,8)
      yr=[0.,1.5]
    end
    8: begin        ; Winter CET
      alltit(i)="Winter Central England Temperature"
      timey=reform(rawdat(0,*))
      ts=mkseason(moncet,9,2)
      yr=[0.,1.5]
    end
  endcase
  print,alltit(i)
  ;
  ; Compute the trend and mean distributions
  ;
  ; Extract pre-1900 period only
  kl=where(timey le 1900,nyr)
  ts=ts(kl)
  timey=timey(kl)
  ;
  ; Make anomalies relative to their whole length
  mkanomaly,ts,refmean=obsmean
  ;
  ; Get the HCCI series and convert that to anomalies too
  modts=reform(hccits(*,i))
  mkanomaly,modts,refmean=modmean
  print,obsmean,modmean
  ;
  ; Define arrays
  nmax=nyr-maxsize+1
  alltren=fltarr(nmax,ntry)*!values.f_nan
  allmean=fltarr(nmax,ntry)*!values.f_nan
  ;
  ; Process each segment length separately
  for j = 0 , ntry-1 do begin
    ;
    distrib_trendmean,ts,timey,seglen=tlen(j),$
      binmean=binmean*binsc(i),bmhalf=bmhalf*binsc(i),$
      bintren=bintren*binsc(i),bthalf=bthalf*binsc(i),$
      tren95=onet95,mean95=onem95,freqmean=onefmean,freqtren=oneftren
    tren95(j,i,0)=onet95
    mean95(j,i,0)=onem95
    freqmean(*,j,i,0)=onefmean
    freqtren(*,j,i,0)=oneftren
    ;
    distrib_trendmean,modts,modtimey,seglen=tlen(j),$
      binmean=binmean*binsc(i),bmhalf=bmhalf*binsc(i),$
      bintren=bintren*binsc(i),bthalf=bthalf*binsc(i),$
      tren95=onet95,mean95=onem95,freqmean=onefmean,freqtren=oneftren
    tren95(j,i,1)=onet95
    mean95(j,i,1)=onem95
    freqmean(*,j,i,1)=onefmean
    freqtren(*,j,i,1)=oneftren
    ;
  endfor
  ;
endfor
;
; Now save results for later plotting
;
save,filename='trendmean_distrib.idlsave',$
  ntry,tlen,nts,alltit,tren95,mean95,binmean,bintren,freqmean,freqtren,$
  nbin,binsc
;
end