;
; Reads in tree ring data from Phil/Keith and converts it to netCDF format,
; with all metadata there, and with number of chronologies extendable.
; JUST DOES SU, BUT BACK FURTHER TO 1400!!!
;
;--------------------------------------------------------------------------
;
; Define files
;
fstem='/cru/keith1/vax/hugswiss/'
fname=['su','su2']
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
  ;
  ; each *.txt file should contain the line with Datens on it at the end.
  ; the first value on this line is equal to the number of chronologies
  ; in the dataset.
  ;
  cmd='grep Datens '+fstem+fname(ifl)+'.txt'
  spawn,cmd,outline
;  print,outline
  neach(ifl)=fix(outline)
  nchron=nchron+neach(ifl)
endfor
print,'Number of chronologies:',nchron
if nchron ne 122 then message,'Unexpected number of chronologies'
;
;Specify times
;
mintime=1400
maxtime=1992
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 then message,'Not found all header information'
;
; 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 then message,'Not found all chronologies'
;
; Create new netCDF file
;
ncid=ncdf_create('tree_dens_su1400.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