;
; Reads in MXD tree ring data from Phil/Keith, including TORNETRASK and
; JASPER from different sources; sorts out missing data, normalises over
; 1881-1960, and saves for later use.
;
;--------------------------------------------------------------------------
;
trv=0           ; selects tree-ring-variable: 0=MXD 1=TRW
case trv of
  0: fnadd='mxd'
  1: fnadd='trw'
endcase
;
; Define files
;
; I NOW HAVE MY COPY OF THESE - SEE ./COPYOFHUGSWISS
fstem='/cru/keith1/vax/hugswiss/'
fname=['gb','eur','accn','usa','su','skan','su2','su3']
nfile=n_elements(fname)
print,'Files to process:',nfile
;
; Count no. of chronologies
;
nchron=0
neach=fltarr(nfile)
for ifl = 0 , nfile-1 do begin
  cmd='grep Datens '+fstem+fname(ifl)+'.txt'
  spawn,cmd,outline
  neach(ifl)=fix(outline)
  nchron=nchron+neach(ifl)
endfor
print,'Number of chronologies:',nchron
if nchron ne 393 then message,'Unexpected number of chronologies'
nchron=nchron+1      ; extra one for Tornetrask to be added later
;
;Specify times
;
mintime=1400
maxfile=1994
maxtime=2000
filetime=maxfile-mintime+1
ntime=maxtime-mintime+1
time=findgen(ntime)+mintime
;
; Define arrays
;
idno=strarr(nchron)
idname=strarr(nchron)
location=strarr(nchron)
country=strarr(nchron)
tree=strarr(nchron)
yrstart=intarr(nchron)
yrend=intarr(nchron)
statlat=fltarr(nchron)
statlon=fltarr(nchron)
density=fltarr(ntime,nchron)*!values.f_nan
fraction=fltarr(ntime,nchron)
;
; Read header info.
;
dummy=strarr(7)
elest=0
for ifl = 0 , nfile-1 do begin
  print,fname(ifl),neach(ifl)
  indat=strarr(12,neach(ifl))
  openr,1,fstem+fname(ifl)+'.txt'
  readf,1,dummy
  readf,1,indat,format='(A7,2A10,A26,A5,A6,A4,1X,A4,A6,A1,A5,A1)'
  close,1
  eleen=elest+neach(ifl)-1
  idno(elest:eleen)=indat(0,*)
  idname(elest:eleen)=indat(1,*)
  location(elest:eleen)=indat(3,*)
  country(elest:eleen)=indat(4,*)
  tree(elest:eleen)=indat(5,*)
  yrstart(elest:eleen)=fix(indat(6,*))
  yrend(elest:eleen)=fix(indat(7,*))
  statlat(elest:eleen)=latlon(indat(8,*),indat(9,*))
  statlon(elest:eleen)=latlon(indat(10,*),indat(11,*))
  elest=eleen+1
endfor
if elest ne nchron-1 then message,'Not found all header information'
;                 ^---- -1 'cos Tornetrask not read yet!!!
; Now make up something for Tornetrask
;
idno(elest)='9999Z'
idname(elest)='TORNXX'
location(elest)='Tornetrask'
country(elest)='SW'
tree(elest)='PISY'
yrstart(elest)=1400          ; actually starts in 441, but we only want 1400-
yrend(elest)=1980
statlat(elest)=68.3
statlon(elest)=19.5
;
; Now read the chronologies
;
elest=0
for ifl = 0 , nfile-1 do begin
  print,fname(ifl),neach(ifl)
  indat=fltarr(neach(ifl)*2+1,filetime)
  infile=fstem+fname(ifl)+string([mintime,maxfile],format='(2I4)')+'.'+fnadd
  print,infile
  openr,1,infile
  readf,1,indat
  close,1
  eleen=elest+neach(ifl)-1
  for ic = elest , eleen do begin
    colno=(ic-elest)*2+1
    density(0:filetime-1,ic)=indat(colno,*)
    fraction(0:filetime-1,ic)=indat(colno+1,*)
    dummy=where(density(*,ic) ne -9.99,ndum)
    if ndum eq 0 then print,fname(ifl),idname(ic),colno,ndum
  endfor
  elest=eleen+1
endfor
if elest ne nchron-1 then message,'Not found all chronologies'
;                 ^---- -1 'cos Tornetrask not read yet!!!
;
; Now normalise these over 1881-1960 (and set missing data to NaN)
;
for i = 0 , nchron-1 do begin
  xxx=reform(density(*,i))
  misslist=where(xxx eq -9.99,nmiss)
  if nmiss gt 0 then xxx(misslist)=!values.f_nan
  mknormal,xxx,time,refperiod=[1881,1960],refmean=rfm,refsd=rfs
;  print,i,nmiss,reform(rfm),reform(rfs)
  density(*,i)=xxx(*)
endfor
;
; Now need to read Tornetrask from a different file
;
infile='/cru/u2/f055/tree5/torn.winfrin.'+fnadd
print,infile
openr,1,infile
; we know the files start at yr 400 or 440, but we want nothing till 1400, so we
; can skip lines (1400-440 [400 for trw])/10 + 1 header line
nmiss=(1400-440)/10 + 1
if trv eq 1 then nmiss=nmiss+4   ; extra 4 lines for AD400,410,420 & 430
dumst=strarr(nmiss)
readf,1,dumst
; we now want all lines (10 yr per line) from 1400 to 1980 ('81 for trw), which is
; (1980-1400)/10 + 1 lines
nwant=(1980-1400)/10 + 1
rawdat=fltarr(20,nwant)
readf,1,rawdat,format='(59(10X,10(I4,I3),/))'
close,1
; separate into timeseries of MXD and of no. of cores
rawdat=reform(rawdat,20*nwant)
rawdat=reform(rawdat,2,10*nwant)
onemxd=reform(rawdat(0,*))
onencore=reform(rawdat(1,*))
; now extend to end at year 2000
onemxd=[onemxd,replicate(9990.,11)]
onencore=[onencore,replicate(0.,11)]
; now set missing equal to NaN in MXD
ml=where(onemxd eq 9990,nmiss)
if nmiss gt 0 then onemxd(ml)=!values.f_nan
; now normalise over 1881-1960
mknormal,onemxd,time,refperiod=[1881,1960]
; now find fraction of cores available
maxcore=max(onencore)
onencore=onencore/maxcore
; now store them at the end
density(*,elest)=onemxd
fraction(*,elest)=onencore
;
; Now need to read Jasper from a different file and use it to replace
; Sunwapta Pass
;
elest=where(location eq 'Sunwapta Pass, N-Passh.   ',nfind)
if nfind eq 0 then message,'Cant find Sunwapta'
location(elest)='Jasper'
yrstart(elest)=1400
;
infile='/cru/u2/f055/tree5/jasper.winfrin.'+fnadd
print,infile
openr,1,infile
; we know the file starts at yr 1070 (or 1072 for trw), but we want nothing
; till 1400, so we can skip lines (1400-1070)/10 + 1 header line
nmiss=(1400-1070)/10 + 1
dumst=strarr(nmiss)
readf,1,dumst
; we now want all lines (10 yr per line) from 1400 to 1991 (or 1987 for trw),
; which is ( [1980 trw] 1990-1400)/10 + 1 lines
; (since 1991 is on line beginning 1990)
nwant=(1990-1400)/10 + 1
if trv eq 1 then nwant=nwant-1       ; for trw, 1 less line
rawdat=fltarr(20,nwant)
readf,1,rawdat,format='(60(10X,10(I4,I3),/))'
close,1
; separate into timeseries of MXD and of no. of cores
rawdat=reform(rawdat,20*nwant)
rawdat=reform(rawdat,2,10*nwant)
onemxd=reform(rawdat(0,*))
onencore=reform(rawdat(1,*))
; now lengthen to end at year 2000
if trv eq 1 then begin       ; currently ends in 1989, so add 11 years
  onemxd=[onemxd,replicate(9990.,11)]
  onencore=[onencore,replicate(0.,11)]
endif
if trv eq 0 then begin       ; currently ends in 1999, so add 1 year
  onemxd=[onemxd,replicate(9990.,1)]
  onencore=[onencore,replicate(0.,1)]
endif
;onemxd=onemxd(0:ntime-1)
;onencore=onencore(0:ntime-1)
; now set missing equal to NaN in MXD
ml=where(onemxd eq 9990,nmiss)
if nmiss gt 0 then onemxd(ml)=!values.f_nan
; now normalise over 1881-1960
mknormal,onemxd,time,refperiod=[1881,1960]
; now find fraction of cores available
maxcore=max(onencore)
onencore=onencore/maxcore
; now store them at Sunwapta pass
density(*,elest)=onemxd
fraction(*,elest)=onencore
;
; Remove entirely missing chronologies (because they didn't have full data
; over the 1881-1960 normalisation period?)
;
help,density
mtot=total(finite(density),1)
checkl=where(mtot eq 0,ncheck)
print,'Excluded because all missing (maybe due to <5 values during 1881-1960 norm)'
for icheck = 0 , ncheck-1 do begin
 print,location[checkl[icheck]],yrstart[checkl[icheck]],yrend[checkl[icheck]]
endfor
kl=where(mtot gt 0,nkeep)
nchron=nkeep
idno=idno(kl)
idname=idname(kl)
location=location(kl)
country=country(kl)
tree=tree(kl)
yrstart=yrstart(kl)
yrend=yrend(kl)
statlat=statlat(kl)
statlon=statlon(kl)
density=density(*,kl)
fraction=fraction(*,kl)
help,density
;
; Now save dataset
;
timey=time
mxd=density
nyr=ntime
message,'NOT OVERWRITING THE DATA SET!'
save,filename='all'+fnadd+'.idlsave',$
  nchron,idno,idname,location,country,tree,yrstart,yrend,statlat,statlon,$
  mxd,fraction,timey,nyr
;
end