;
; Reads in tree ring data from Phil/Keith and converts it to netCDF format,
; with all metadata there, and with number of chronologies extendable.
;
; THIS DOES 3 LONG TIMESERIES, AND ADJUSTS THEIR VARIANCES ACCORDING TO RBAR
; AND THE NUMBER OF CORES IN EACH CHRONOLOGY.
;
;--------------------------------------------------------------------------
;
; Define files
;
fstem='/dave3/f023/vax/longvolc/'
fname=['long']
nfile=n_elements(fname)
print,'Files to process:',nfile
;
; Count no. of chronologies
;
nchron=3           ; for this case, I know it is 3!
neach=fltarr(nfile)
neach(0)=3
print,'Number of chronologies:',nchron
if nchron ne 3 then message,'Unexpected number of chronologies'
;
;Specify times
;
mintime=500
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.  - there is none here!
;
;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)
  print,fstem+fname(ifl)+string([mintime,maxtime],format='(I3,I4)')+'.mxd'
  openr,1,fstem+fname(ifl)+string([mintime,maxtime],format='(I3,I4)')+'.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'
;
; Adjust the variance according to sample size
;
maxn=[31,19,42]
rbar=[0.57,0.443,0.4]
for i = 0 , nchron-1 do begin
  xin=density(*,i)
  nin=fraction(*,i)
  misslist=where(nin eq 0,nmiss)
  xout=mkeffective(xin,nin*maxn(i),rbar=rbar(i),neff=nout)
  if nmiss gt 0 then begin
    xout(misslist)=-9.99
    nout(misslist)=0.
  endif
  density(*,i)=xout
  fraction(*,i)=nout
endfor
;
; Create new netCDF file
;
ncid=ncdf_create('longvolc.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