;
; As a check on the SOAP European experiment data, I will now read in the
; gridded MXD, the gridded warm and cold season temperatures, form full
; and two half regional means of each and correlate/regress them against
; each other.
;
maxbox=66
fn=['treedensity_europe_WITHLOWFREQ.dat',$
  'instrumentaltemp_'+['April-September','October-March']+'.dat']
nhead=[  8,   5,   5]
yrst=[1400,1851,1851]
yren=[1989,2000,2000]
;
nvar=n_elements(fn)
minyr=min(yrst)
maxyr=max(yren)
maxlen=maxyr-minyr+1
ndec=ceil(float(yren-yrst+1)/10.)
timey=findgen(maxlen)+minyr
nbox=lonarr(nvar)
xlon=fltarr(maxbox,nvar)*!values.f_nan
ylat=fltarr(maxbox,nvar)*!values.f_nan
fd=fltarr(maxbox,maxlen,nvar)*!values.f_nan
;
for ivar = 0 , nvar-1 do begin
  ;
  print,fn[ivar]
  openr,1,fn[ivar]
  ;
  headst=strarr(nhead[ivar])
  readf,1,headst
  ;
  readf,1,nbox1,format='(16X,I5)'
  readf,1,missval,format='(16X,F5.2)'
  headst=''
  readf,1,headst
  if nbox1 gt maxbox then message,'Too many boxes'
  nbox[ivar]=nbox1
  ;
  for ibox = 0 , nbox1-1 do begin
    readf,1,xlon1,ylat1
    xlon[ibox,ivar]=xlon1
    ylat[ibox,ivar]=ylat1
    for idec = 0 , ndec[ivar]-1 do begin
      dat1=fltarr(10)
      readf,1,dat1,format='(9X,10F7.2)'
      yr1=yrst[ivar]-minyr+idec*10
      yr2=yr1+9
      fd[ibox,yr1:yr2,ivar]=dat1[*]
    endfor
    misslist=where(fd[ibox,*,ivar] eq missval,nmiss)
    if nmiss gt 0 then fd[ibox,misslist,ivar]=!values.f_nan
  endfor
  ;
  close,1
  ;
endfor
;
; Form three regional means from each data set, with cos(latitude) weighting
; because the data represent grid-box-averages and thus the area that they
; represent is proportional to cos(latitude)
;
regname=['ALL-EUROPE','NORTH-EUROPE','SOUTH-EUROPE']
nreg=n_elements(regname)
;
regts=fltarr(maxlen,nvar,nreg)*!values.f_nan
for ivar = 0 , nvar-1 do begin
  vn=nbox[ivar]
  vlon=xlon[0:vn-1,ivar]
  vlat=ylat[0:vn-1,ivar]
  vfd=reform(fd[0:vn-1,*,ivar])
  for ireg = 0 , nreg-1 do begin
    case ireg of
      0: keeplist=lindgen(vn)
      1: keeplist=where(vlat gt 53.)
      2: keeplist=where(vlat le 53.)
    endcase
print,regname[ireg]
help,keeplist
    vlat1=vlat[keeplist]
    vwt=cos(vlat1*!dtor)
    vfd1=vfd[keeplist,*]
    regts[*,ivar,ireg]=var_adjust(vfd1,timey,weight=vwt,/no_anomaly,$
      refperiod=[1881,1960],/td_rbar,/td_sdev,rbar=keeprbar,/mkrbar,$
      corrmat=keepcorr,/mkcorr)
  endfor
endfor
;
; Now prepare for plotting
;
loadct,39
multi_plot,nrow=3,layout='large'
if !d.name eq 'X' then begin
  window,ysize=800,xsize=1400
  !p.font=-1
endif else begin
  !p.font=0
  device,/helvetica,/bold,font_size=12
endelse
def_1color,20,color='red'
;
for ireg = 0 , nreg-1 do begin
  print,regname[ireg]
  ts1=reform(regts[*,0,ireg])
  mknormal,ts1,timey,refperiod=[1881,1960]
  filter_cru,/nan,10.,tsin=ts1,tslow=ts1l
  cpl_barts,timey,ts1,zeroline=ts1l,$
    /xstyle,xtitle='Year',$
    title=regname[ireg]
  oplot,timey,ts1l,thick=2
  oplot,!x.crange,[0,0],linestyle=1
  for ivar = 1 , 2 do begin
    ts2=reform(regts[*,ivar,ireg])
    kl=where(finite(ts1+ts2) and (timey ge 1881) and (timey le 1960))
    print,correlate(ts1[kl],ts2[kl])
    mknormal,ts2,timey,refperiod=[1881,1960]
    if ivar eq 1 then oplot,timey,ts2,color=20
  endfor
endfor
;
end