;
; Reads in tree ring data from Phil/Keith and converts it to netCDF format,
; with all metadata there, and with number of chronologies extendable.
; INCLUDES TORNETRASK and JASPER
; THIS VERSION ALSO NORMALISES EVERYTHING w.r.t. 1881-1960.
;
;--------------------------------------------------------------------------
;
; Define files
;
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 391 then message,'Unexpected number of chronologies'
nchron=nchron+1      ; extra one for Tornetrask to be added later
;
;Specify times
;
mintime=1400
maxtime=1994
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)
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,ntime)
  openr,1,fstem+fname(ifl)+string([mintime,maxtime],format='(2I4)')+'.mxd'
  readf,1,indat
  close,1
  eleen=elest+neach(ifl)-1
  for ic = elest , eleen do begin
    colno=(ic-elest)*2+1
    density(*,ic)=indat(colno,*)
    fraction(*,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
;
openr,1,'torn.winfrin'
; we know the file starts at yr 440, but we want nothing till 1400, so we
; can skill lines (1400-440)/10 + 1 header line
nmiss=(1400-440)/10 + 1
dumst=strarr(nmiss)
readf,1,dumst
; we now want all lines (10 yr per line) from 1400 to 1980, 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 1994
onemxd=[onemxd,replicate(9990.,5)]
onencore=[onencore,replicate(0.,5)]
; 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 1901-40
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
;
openr,1,'jasper.winfrin'
; we know the file starts at yr 1070, but we want nothing till 1400, so we
; can skill 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, which is
; (1990-1400)/10 + 1 lines  (since 1991 is on line beginning 1990)
nwant=(1990-1400)/10 + 1
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 shorten to end at year 1994
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 1901-40
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
;
; Create new netCDF file
;
ncid=ncdf_create('tree_dens_nh_18811960.nc')
;
; Define dimensions
;
shortid=ncdf_dimdef(ncid,'shortstring',10)
longid=ncdf_dimdef(ncid,'longstring',26)
timid=ncdf_dimdef(ncid,'time',ntime)
staid=ncdf_dimdef(ncid,'station',/unlimited)
;
; Define variables and their attributes
;   NB. IDL converts strings to char (byte) automatically, so no need to
;   do it.  But, the lengths (shortstring and longstring) must be sufficient.
;   When reading it back, it gives byte arrays, but a simple =string()
;   converts it back to strings - including the lengths correctly!
;
statidid=ncdf_vardef(ncid,'statid',[shortid,staid],/char)
;
statabid=ncdf_vardef(ncid,'statabbr',[shortid,staid],/char)
;
statloid=ncdf_vardef(ncid,'statloc',[longid,staid],/char)
;
statcoid=ncdf_vardef(ncid,'country',[shortid,staid],/char)
;
treeid=ncdf_vardef(ncid,'tree',[shortid,staid],/char)
;
yrsid=ncdf_vardef(ncid,'yrstart',[staid],/float)
;
yreid=ncdf_vardef(ncid,'yrend',[staid],/float)
;
latid=ncdf_vardef(ncid,'latitude',[staid],/float)
ncdf_attput,ncid,latid,'Unit','Degrees and hundredths (+ve N, -ve S)'
;
lonid=ncdf_vardef(ncid,'longitude',[staid],/float)
ncdf_attput,ncid,lonid,'Unit','Degrees and hundredths (+ve E, -ve W)'
;
yrid=ncdf_vardef(ncid,'year',[timid],/float)
;
densid=ncdf_vardef(ncid,'density',[timid,staid],/float)
ncdf_attput,ncid,densid,'Source','CRU tree ring chronology database'
ncdf_attput,ncid,densid,'Title','Normalised tree ring density'
ncdf_attput,ncid,densid,'missing_value',-9.99
;
fracid=ncdf_vardef(ncid,'fraction',[timid,staid],/float)
ncdf_attput,ncid,fracid,'Source','CRU tree ring chronology database'
ncdf_attput,ncid,fracid,'Title','Fraction of cores with data'
;
; Switch from define mode to data mode
;
ncdf_control,ncid,/endef
;
; Output the variables
;
ncdf_varput,ncid,statidid,idno
ncdf_varput,ncid,statabid,idname
ncdf_varput,ncid,statloid,location
ncdf_varput,ncid,statcoid,country
ncdf_varput,ncid,treeid,tree
ncdf_varput,ncid,yrsid,yrstart
ncdf_varput,ncid,yreid,yrend
ncdf_varput,ncid,latid,statlat
ncdf_varput,ncid,lonid,statlon
ncdf_varput,ncid,yrid,time
ncdf_varput,ncid,densid,density
ncdf_varput,ncid,fracid,fraction
;
; Close the netCDF file
;
ncdf_close,ncid
;
end