;
; Plots density 'decline' as a time series of the difference between
; temperature and density averaged over the region north of 50N,
; and an associated pattern in the difference field.
; The difference data set is computed using only boxes and years with
; both temperature and density in them - i.e., the grid changes in time.
; The pattern is computed by correlating and regressing the *filtered*
; time series against the unfiltered difference data set.
;
; Use the calibrate MXD after calibration coefficients estimated for 6
; northern Siberian boxes that had insufficient temperature data.
print,'Reading MXD and temperature data'
restore,filename='calibmxd2.idlsave'
; Gets:  g,mxdyear,mxdnyr,fdcorr,fdalph,fdbeta,fdvexp,fdcalib,mxdfd2,$
;        fdrver,fdvver,fdseas,timey
;
print,'Finding region with data'
outfd=total(finite(mxdfd2),3)
list4=where(outfd eq 0,n4)
list2=where(outfd gt 0,n2)
outfd(list4)=4.
outfd(list2)=2.
fdmask=fltarr(g.nx,g.ny)
fdmask(*,*)=!values.f_nan
fdmask(list2)=1.
;
; Prepare for plotting
;
loadct,39
multi_plot,nrow=4
if !d.name eq 'X' then begin
  window,ysize=850
  !p.font=-1
endif else begin
  !p.font=0
  device,/helvetica,/bold,font_size=9
endelse
;
; Identify the period of overlap
;
print,'Computing difference'
yrstart=timey(0)
yrend=mxdyear(mxdnyr-1)
nlap=yrend-yrstart+1
ktem=where(timey le yrend)
kmxd=where(mxdyear ge yrstart)
;
; Compute the difference dataset
;
fddiff=fdcalib(*,*,kmxd)-fdseas(*,*,ktem)
;
; Compute means over the region north of 50N
;
print,'Difference time series'
ky=where(g.y ge 50)
fd=fddiff(*,ky,*)
diffts=globalmean(fd,g.y(ky))
print,'Temperature time series (masked)'
fd=fdseas(*,ky,*)
temts=globalmean(fd,g.y(ky),mask=fdmask(*,ky))
print,'MXD time series'
fd=fdcalib(*,ky,*)
mxdts=globalmean(fd,g.y(ky))
diffts2=mxdts(kmxd)-temts(ktem)
;
thalf=10
filter_cru,thalf,/nan,tsin=diffts,tslow=difflow
filter_cru,thalf,/nan,tsin=temts,tslow=temlow
filter_cru,thalf,/nan,tsin=mxdts,tslow=mxdlow
filter_cru,thalf,/nan,tsin=diffts2,tslow=difflow2
plot,mxdyear,mxdlow,$
  /xstyle,xrange=[1400,2000],xtitle='Year',$
  /ystyle,yrange=[-0.9,0.9],ytitle='>50N masked temperature (!Uo!NC wrt 1961-90)'
oplot,timey,temlow,thick=3
oplot,!x.crange,[0.,0.],linestyle=1
legend,/horiz,['>50N temperature','>50N MXD reconstruction'],thick=[3,1]
;
plot,timey(ktem),difflow2,$
  /xstyle,xrange=[1400,2000],xtitle='Year',$
  /ystyle,yrange=[-1.1,0.6],ytitle='Temperautre difference (!Uo!NC)'
oplot,timey(ktem),difflow,thick=3
oplot,!x.crange,[0.,0.],linestyle=1
legend,/horiz,['Mean of differences','Difference of means'],thick=[3,1]
;
; Now fit a 2nd degree polynomial to the decline series, and then extend
; it at a constant level back to 1400.  In fact we compute its mean over
; 1856-1930 and use this as the constant level from 1400 to 1930.  The
; polynomial is fitted over 1930-1994, forced to have the constant value
; in 1930.
;
ycon=1930
declinets=fltarr(mxdnyr)*!values.f_nan
allx=timey(ktem)
ally=difflow
;
kconst=where(allx le ycon,nconst)
cval=total(ally(kconst))/float(nconst)
kkk=where(mxdyear le ycon)
declinets(kkk)=cval
;
kcurve=where((allx ge ycon) and (allx le 1994),ncurve)
allwt=replicate(1.,ncurve)
acoeff=[0.,1./7200.]
vval=curvefit(allx(kcurve),ally(kcurve),allwt,acoeff,$
  function_name='funct_decline')
kkk=where(mxdyear ge ycon)
declinets(kkk)=vval(*)
;
oplot,mxdyear,declinets,thick=5
;
; Now compute the correlation and regression coefficients between the
; difference time series and the field of differences.
; Time series is smoothed, and difference can be too
;
print,'Difference pattern'
fd=reform(fddiff,g.nx*g.ny,nlap)
print,n2
for i = 0 , n2-1 do begin
  print,i,format='($,I4)'
  xxx=reform(fd(list2(i),*))
  filter_cru,thalf,/nan,tsin=xxx,tslow=yyy
  fd(list2(i),*)=yyy(*)
endfor
print
ylatdummy=fltarr(g.nx*g.ny)
patterndata,fd,difflow,zcoeff=zcoeff,zcorr=zcorr,ylat=ylatdummy,thresh=0.05
zcoeff=reform(zcoeff,g.nx,g.ny)
zcorr=reform(zcorr,g.nx,g.ny)
;
!p.multi=[0,1,2,0,0]
map=def_map(/npolar) & map.limit(0)=25.
labels=def_labels(/off)
;
pause
inter_boxfd,zcorr,g.x,g.y,$
  labels=labels,map=map,$
  levels=[-1.,0,0.2,0.3,0.4,0.5,0.6,0.7,0.9],/scale,$
  title='Correlation with difference series'
fdin={ fd: outfd, x: g.x, y: g.y, nx: g.nx, ny: g.ny }
whizz_fd,fdin=fdin,fdout=f,limit=map.limit
boxplot,f.fd,f.x,f.y,/overplot,/overmap,highlight=f.fd,thick=1.5
;
inter_boxfd,zcoeff,g.x,g.y,$
  labels=labels,map=map,$
  levels=[-1.5,0,0.2,0.4,0.8,1.2,1.6,2.0,4.0],/scale,$
  title='Regression with difference series'
boxplot,f.fd,f.x,f.y,/overplot,/overmap,highlight=f.fd,thick=1.5
;
; Now need to infill the 6 boxes without regression coefficients.  This is
; done using a gaussian-weighted mean of nearby coefficients.
;
print,'Infilling coefficients'
allx=fltarr(g.nx,g.ny)
for iy = 0 , g.ny-1 do allx(*,iy)=g.x(*)
ally=fltarr(g.nx,g.ny)
for ix = 0 , g.nx-1 do ally(ix,*)=g.y(*)
statx=allx(list2)
staty=ally(list2)
;
statval=zcoeff(list2)
misslist=where(finite(statval) eq 0,nmiss)
if nmiss gt 0 then begin
  fd_extend,statval,statx,staty,search=1000.,wavelen=400.
  zcoeff(list2(misslist))=statval(misslist)
endif
;
pause
inter_boxfd,zcoeff,g.x,g.y,$
  labels=labels,map=map,$
  levels=[-1.5,0,0.2,0.4,0.8,1.2,1.6,2.0,4.0],/scale,$
  title='Regression with difference series (infilled)'
boxplot,f.fd,f.x,f.y,/overplot,/overmap,highlight=f.fd,thick=1.5
;
; Now apply I completely artificial adjustment for the decline
; (only where coefficient is positive!)
;
tfac=declinets-cval
fdcorrect=fdcalib
for iyr = 0 , mxdnyr-1 do begin
  fdcorrect(*,*,iyr)=fdcorrect(*,*,iyr)-tfac(iyr)*(zcoeff(*,*) > 0.)
endfor
;
; Now save the data for later analysis
;
save,filename='calibmxd3.idlsave',$
  g,mxdyear,mxdnyr,fdcalib,mxdfd2,fdcorrect
;
end