; ; Reads in tree ring data from Phil/Keith and converts it to netCDF format, ; with all metadata there, and with number of chronologies extendable. ; ;-------------------------------------------------------------------------- ; ; 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 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=[1901,1940] ; 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=[1901,1940] ; 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.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