;
; Writes an ASCII file with data (gridded, not extended, not corrected,
; possibly ABD-adjusted, calibrated, but then possibly normalised) for input
; to the European temperature experiment of SOAP WP3.
;
doinfill=0         ; use PCR-infilled data or not?
doabd=1            ; use ABD-adjusted data or not?
docorr=0           ; use corrected version or not? (uncorrected only available
                   ;                                for doinfill=doabd=0)
;
doreg=1     ; 0=ALL, 1=Greater Europe
donorm=1    ; 0=leave calibrated, 1=remove calibration by normalisation
;
missval=-9.99
;
; Get the calibrated data
;
print,'Reading reconstructions'
if doabd eq 0 then begin
  if doinfill eq 0 then begin
    restore,'calibmxd5.idlsave'
    ; Gets: g,mxdyear,mxdnyr,fdcalibu,fdcalibc,mxdfd2,timey,fdseas
    if docorr eq 0 then fdcalibc=fdcalibu
  endif else begin
    restore,'calibmxd5_pcr.idlsave'
    ; Gets: g,mxdyear,mxdnyr,fdcalibc,timey,fdseas
  endelse
endif else begin
  if doinfill eq 0 then begin
    print,'THIS IS THE ONE I WANT!'
    restore,'../mann/calibmxd5_abdlow.idlsave'
    ; Gets: g,mxdyear,mxdnyr,fdcalibu
    if docorr eq 0 then fdcalibc=fdcalibu
  endif else begin
    restore,'../mann/calibmxd5_abdlow_pcr.idlsave'
    ; Gets: g,mxdyear,mxdnyr,fdcalibu
  endelse
endelse
;
map=def_map(/npolar)
;
case doreg of
  0: begin
    print,'OUTPUTTING ALL DATA'
  end
  1: begin
    print,'OUTPUTTING GREATER EUROPEAN DATA ONLY'
    missx=where(g.x lt -30)
    fdcalibc[missx,*,*]=!values.f_nan
    missx=where(g.x gt 61)
    fdcalibc[missx,*,*]=!values.f_nan
    missx=where(g.x eq 57.5) & missx=missx[0]
    missy=where(g.y eq 52.5) & missy=missy[0]
    fdcalibc[missx,missy,*]=!values.f_nan
    map.limit=[20,-60,90,90]
  end
endcase
map.origx=mean(map.limit[[1,3]])
;
; Lets plot the location of the boxes with data
;
boxlist=where(total(finite(fdcalibc),3) gt 0,nbox)
fdmask=fltarr(g.nx,g.ny)*!values.f_nan
fdmask[boxlist]=20
inter_boxfd,fdmask,g.x,g.y,map=map
;
; Let's plot the average of all plot
;
pause
xxx=reform(fdcalibc,g.nx*g.ny,mxdnyr)
plot,total(xxx,1,/nan)/float(total(finite(xxx),1))
;
; Now lets output the data to an ASCII file
;
ndec=ceil(mxdnyr/10.)
;
if doabd ne 0 then begin
  openw,1,'treedensity_europe_WITHLOWFREQ.dat'
endif else begin
  openw,1,'treedensity_europe.dat'
endelse
if donorm eq 0 then begin
  printf,1,'These are gridded, calibrated estimates of mean warm-season'
  printf,1,'(April-September) temperature anomalies (degC with respect to the'
  printf,1,'1961-1990 mean), based on tree-ring density records from within'
  printf,1,'each grid box (all boxes contained at least 1 tree-ring'
  printf,1,'chronology).'
endif else begin
  printf,1,'These are gridded, uncalibrated anomalies of maximum latewood'
  printf,1,'density, based on tree-ring density records from within'
  printf,1,'each grid box (all boxes contained at least 1 tree-ring'
  printf,1,'chronology).'
endelse
if doabd ne 0 then begin
  printf,1,'THIS VERSION CONTAINS GRID-BOX TIME SERIES WHOSE LOW-FREQUENCIES'
  printf,1,'HAVE BEEN "RESTORED" TO MATCH THOSE OBTAINED AT REGIONAL SCALES'
  printf,1,'BY BRIFFA ET AL. (2001) JGR, USING AGE-BAND DECOMPOSITION (ABD).'
endif
printf,1
printf,1,nbox,format='("Number of boxes:",I5)'
printf,1,missval,format='("Missing value:",F7.2)'
printf,1
for ix = 0 , g.nx-1 do begin
  for iy = 0 , g.ny-1 do begin
    if finite(fdmask[ix,iy]) then begin
      printf,1,g.x[ix],g.y[iy],format='(2F10.4,"  Longitude & Latitude")'
      ts1=reform(fdcalibc[ix,iy,*])
      if donorm eq 1 then mknormal,ts1,mxdyear,refperiod=[1881,1960]
      ml=where(finite(ts1) eq 0,nmiss)
      ts1[ml]=missval
      ;
      for idec = 0 , ndec-2 do begin    ; don't do last (incomplete) decade
        iyr=idec*10
        jyr=(iyr+9) < (mxdnyr-1)
        printf,1,mxdyear[[iyr,jyr]],ts1[iyr:jyr],format='(I4,":",I4,10F7.2)'
      endfor
    endif
  endfor
endfor
close,1
flush,-1
;
end