;
; Reads in the gridded MXD data, and outputs the (non-infilled)
; version to an ascii file, extracting only those grid boxes
; that have at least some data.  (or optionally, outputting them all).
;
if n_elements(doabd) eq 0 then doabd=1    ; 0=hug, 1=abd
if n_elements(dopcr) eq 0 then dopcr=0    ; 0=boxes with trees, 1=pcr extended
if n_elements(dogrid) eq 0 then dogrid=0    ; 0=just list timeseries, 1=do grid
;
trv=0           ; selects tree-ring-variable: 0=MXD 1=TRW
case trv of
  0: begin
    fnadd='mxd'
    end
  1: begin
    fnadd='trw'
    end
endcase
titadd=strupcase(fnadd)
;
if (trv eq 1) and (doabd ne 0) then message,'Age-banded TRW not available'
;
; Decide what file to read
;
if doabd eq 0 then begin
  fnin=fnadd+'_infill.idlsave'
  print,fnin
  restore,fnin
  ; Gets: ntime,maxbox,boxdens,g,keepi,keepj,timey
  ;       baselist,nbase,mxdfd,mxdyear,mxdnyr,mxdfd2,beglist,endlist
endif else begin
  if dopcr eq 0 then begin
    fnin='../mann/calibmxd5_abdlow.idlsave'
    print,fnin
    restore,fnin
    ; Gets: g,mxdyear,mxdnyr,fdcalibu
    mxdfd=fdcalibu
  endif else begin
    fnin='../mann/calibmxd5_abdlow_pcr.idlsave'
    print,fnin
    restore,fnin
    ; Gets: g,mxdyear,mxdnyr,fdcalibu
    mxdfd=fdcalibc
    ; Because we only have the "corrected" version for this case, we ought
    ; not to output anything after 1960.
    dummy=where(mxdyear le 1960,mxdnyr)
    mxdyear=mxdyear[0:mxdnyr-1]
    mxdfd=mxdfd[*,*,0:mxdnyr-1]
  endelse
endelse
;
; Plot coverage through time
;
numts=total(total(finite(mxdfd),1),1)
plot,mxdyear,numts,/xstyle
;
; Open output file
;
if doabd eq 0 then begin
  fnout='schweingruber_mxd_grid.dat'
endif else begin
  if dopcr eq 0 then fnout='schweingruber_mxdabd_grid.dat' $
                else fnout='tett_mxdabd_grid.dat'
endelse
print,fnout
openw,1,fnout
;
; Either output as columns for each series with data, or as a 2D grid
; through time
;
if dogrid eq 0 then begin
  xfd=fltarr(g.nx,g.ny)
  for iy = 0 , g.ny-1 do xfd[*,iy]=g.x
  yfd=fltarr(g.nx,g.ny)
  for ix = 0 , g.nx-1 do yfd[ix,*]=g.y
  ;
  numfd=total(finite(mxdfd),3)
  kl=where(numfd gt 0,nbox)
  print,'Number of boxes with some MXD data = ',nbox
  if nbox ne 115 then message,'Ooops!'
  ;
  xfd=reform(xfd,g.nx*g.ny)
  yfd=reform(yfd,g.nx*g.ny)
  printf,1,'Longitudes'
  printf,1,xfd[kl],format='(115F8.3)'
  ;
  printf,1,'Latitudes'
  printf,1,yfd[kl],format='(115F8.3)'
  ;
  printf,1,'Gridded MXD data (columns are grid boxes, rows are years)'
  mxdfd=reform(mxdfd,g.nx*g.ny,mxdnyr)
  mxdlist=mxdfd(kl,*)
  ml=where(finite(mxdlist) eq 0,nmiss)
  if nmiss gt 0 then mxdlist(ml)=-9.99
  printf,1,mxdlist,format='(115F8.3)'
endif else begin
  ml=where(finite(mxdfd) eq 0,nmiss)
  if nmiss gt 0 then mxdfd(ml)=-9.99
  printf,1,'Osborn et al. (2004) gridded reconstruction of warm-season'
  printf,1,'(April-September) temperature anomalies (from the 1961-1990 mean).'
  printf,1,'Reconstruction is based on tree-ring density records.'
  printf,1
  printf,1,'NOTE: recent decline in tree-ring density has been ARTIFICIALLY'
  printf,1,'REMOVED to facilitate calibration.  THEREFORE, post-1960 values'
  printf,1,'will be much closer to observed temperatures then they should be,'
  printf,1,'which will incorrectly imply the reconstruction is more skilful'
  printf,1,'than it actually is.  See Osborn et al. (2004).'
  printf,1
  printf,1,'Osborn TJ, Briffa KR, Schweingruber FH and Jones PD (2004)'
  printf,1,'Annually resolved patterns of summer temperature over the Northern'
  printf,1,'Hemisphere since AD 1400 from a tree-ring-density network.'
  printf,1,'Submitted to Global and Planetary Change.'
  ;
  printf,1,'Grid resolution is 5 by 5 degrees, for the Northern Hemisphere only,'
  printf,1,'with the first value at ',g.x[0],' W',g.y[0],' N,'
  printf,1,'then moving eastwards to ',g.x[g.nx-1],' and then southwards to ',g.y[g.ny-1],'.'
  printf,1
  printf,1,'Missing value is -9.99'
  printf,1
  for iyr = 0 , mxdnyr-1 do begin
    print,mxdyear[iyr],format='($,I5)'
    printf,1,mxdyear[iyr],format='(I4)'
    printf,1,mxdfd[*,*,iyr],format='(10F8.2)'
  endfor
  print
endelse
;
close,1
;
end