;
; 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 (or filtered) difference data set.
;
;*** MUST ALTER FUNCT_DECLINE.PRO TO MATCH THE COORDINATES OF THE
; START OF THE DECLINE ***  ALTER THIS EVERY TIME YOU CHANGE ANYTHING ***
;
matchvar=1        ; 0=regression, 1=variance matching
;
doadjust=1        ; 0=do not, 1=show some more results,
                  ; 2=do go on and produce a new data set with adjustment for decline
;
pantit='('+['a','b','c']+') '
panoff=1
;
; Use the calibrate MXD after calibration coefficients estimated for 6
; northern Siberian boxes that had insufficient temperature data.
print,'Reading MXD and temperature data'
if matchvar eq 0 then fnadd='_regress' else fnadd=''
restore,filename='calibmxd2'+fnadd+'.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)
; Where is MXD always missing?
list4=where(outfd eq 0,n4)
; Where is there sometimes some MXD?
list2=where(outfd gt 0,n2)
; Outfd is a mask of 4's and 2's for producing an outline of the region
; with at least some MXD
outfd(list4)=4.
outfd(list2)=2.
; fdmask is the region where there is at least some MXD
fdmask=fltarr(g.nx,g.ny)
fdmask(*,*)=!values.f_nan
fdmask(list2)=1.
;
; Prepare for plotting
;
loadct,39
multi_plot,nrow=3,layout='large'
if !d.name eq 'X' then begin
  window,ysize=750
  !p.font=-1
endif else begin
  !p.font=0
  device,/helvetica,/bold,font_size=16
endelse
def_1color,20,color='red'
def_1color,21,color='orange'
def_1color,22,color='yellow'
def_1color,23,color='green'
def_1color,24,color='lblue'
def_1color,25,color='deepblue'
def_1color,26,color='vlpurple'
def_1color,28,color='black'
def_smearcolor,fromto=[26,28]
cc=27-findgen(8)
;
; 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)
;
; bothmask is the time-dependent mask where both data sets have values
bothmask=float(finite(fddiff))    ; must convert byte array back to float
ml=where(bothmask eq 0,nmiss)
if nmiss gt 0 then bothmask[ml]=!values.f_nan
;
; Compute means over the region north of 50N
;
; (1) The mean of the difference series (implicit time-varying MXD&T mask)
print,'Difference time series'
ky=where(g.y ge 50)
fd=fddiff(*,ky,*)
diffts=globalmean(fd,g.y(ky))
;
; (2) Mean of temperature (explicit time-constant MXD mask, implicit time-varying T mask)
print,'Temperature time series (masked)'
fd=fdseas(*,ky,*)
temts=globalmean(fd,g.y(ky),mask=fdmask(*,ky))
;
; (3) Mean of MXD (implicit time-varying MXD mask)
print,'MXD time series'
fd=fdcalib(*,ky,*)
mxdts=globalmean(fd,g.y(ky))
;
; (4) Difference of the mean series
diffts2=mxdts(kmxd)-temts(ktem)
;
; (5) Mean of MXD (explicit time-varying MXD&T mask)
fd=fdcalib(*,ky,*)
fd=fd[*,*,kmxd]
maskmxdts=globalmean(fd,g.y[ky],mask=bothmask[*,ky,*])
;
; (6) Mean of T (explicit time-varying MXD&T mask)
fd=fdseas[*,ky,*]
fd=fd[*,*,ktem]
masktemts=globalmean(fd,g.y[ky],mask=bothmask[*,ky,*])
;
; (7) Difference of masked mean series
diffts3=maskmxdts-masktemts
;
; (8) Mean of normalised MXD (explicit time-varying MXD&T mask) then calibrated
fd=reform(fdcalib,g.nx*g.ny,mxdnyr)
mkanomaly,fd,mxdyear,refperiod=[1881,1960]
fd=reform(fd,g.nx,g.ny,mxdnyr)
fd=fd[*,ky,*]
fd=fd[*,*,kmxd]
; With implicit time-varying MXD mask
mxdarea=globalmean(fd,g.y[ky])
areacoeff=linfit(mxdarea,temts[ktem])
r=correlate(mxdarea,temts)
print,'Area-average first, then skill (unmasked) =',r
; With explicit time-varying MXD&T mask
maskmxdarea=globalmean(fd,g.y[ky],mask=bothmask[*,ky,*])
maskareacoeff=linfit(maskmxdarea,masktemts)
r=correlate(maskmxdarea,masktemts)
print,'Area-average first, then skill   (masked) =',r
; Now try high-pass filtering before calibration
filter_cru,40.,/nan,tsin=masktemts,tshigh=targts
filter_cru,40.,/nan,tsin=maskmxdarea,tshigh=predts
r1=correlate(predts,targts)
;kcalib=where((timey[ktem] ge 1911) and (timey[ktem] le 1990))
kcalib=where((timey[ktem] ge 1856) and (timey[ktem] le 1960))
r2=correlate(predts[kcalib],targts[kcalib])
print,'hi freq r = (all, calib)',r1,r2
hiareacoeff=linfit(predts[kcalib],targts[kcalib])
;hiareacoeff=linfit(predts,targts)
arearecon=hiareacoeff[1]*maskmxdarea
; Make the means match over the calibration period
mkanomaly,arearecon,mxdyear[kmxd],refperiod=[1911,1990]
;mkanomaly,arearecon,mxdyear[kmxd],refperiod=[1881,1960]
dummy=masktemts
mkanomaly,dummy,timey[ktem],refperiod=[1911,1990],refmean=rfm
;mkanomaly,dummy,timey[ktem],refperiod=[1881,1960],refmean=rfm
arearecon=arearecon+rfm[0]
;
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
filter_cru,thalf,/nan,tsin=diffts3,tslow=difflow3
filter_cru,thalf,/nan,tsin=maskmxdts,tslow=maskmxdlow
filter_cru,thalf,/nan,tsin=masktemts,tslow=masktemlow
filter_cru,thalf,/nan,tsin=arearecon,tslow=areareconlow
print,'Correlation between masked T & MXD',correlate(maskmxdts[kcalib],masktemts[kcalib])
print,'Correlation between smoothed masked T & MXD',correlate(maskmxdlow[kcalib],masktemlow[kcalib])
plot,mxdyear,mxdlow,/nodata,$
  title=pantit[0+panoff],$
;  ymargin=[6,10],$
  /xstyle,xrange=[1850,2000],xtitle='Year',$
  /ystyle,yrange=[-0.8,0.7],ytitle='>50N masked temperature (!Uo!NC wrt 1961-90)'
;oplot,timey,temlow,thick=3,color=20
oplot,mxdyear[kmxd],maskmxdlow,thick=2
oplot,timey[ktem],masktemlow,thick=4,color=20
;oplot,mxdyear[kmxd],areareconlow,thick=2,color=21
oplot,!x.crange,[0.,0.],linestyle=1
legend,horiz=0,$
  ['>50N temperature (masked)','>50N MXD reconstruction (masked)','>50N MXD reconstruction (masked then scaled)'],$
  thick=[4,2,2],color=[20,!p.color,23]
;
; Try regressing mean MXD against mean T to get a better fit
;
kcalib=where((timey[ktem] ge 1856) and (timey[ktem] le 1960))
;bestcoeff=linfit(maskmxdts,masktemts)
bestcoeff=linfit(maskmxdts[kcalib],masktemts[kcalib])
r=correlate(maskmxdts[kcalib],masktemts[kcalib])
print,'Area-average later, then skill   (masked) =',r
print,'Best-fit scaling coefficient=',bestcoeff[1]
oplot,mxdyear[kmxd],bestcoeff[0]+bestcoeff[1]*maskmxdlow,thick=2,color=23
;
if doadjust gt 0 then begin
;
plot,timey(ktem),difflow2,/nodata,$
  title=pantit[1+panoff],$
  /xstyle,xrange=[1850,2000],xtitle='Year',$
  /ystyle,yrange=[-1.,0.5],ytitle='Temperature difference (!Uo!NC)'
oplot,timey(ktem),difflow,thick=3
;oplot,timey(ktem),difflow3,thick=3,color=21
oplot,!x.crange,[0.,0.],linestyle=1
;legend,/horiz,['Mean of differences','Difference of means'],thick=[3,1]
;legend,/horiz,['Difference of means'],thick=[3]
;
; 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.
;
if matchvar eq 0 then ycon=1900 else ycon=1930   ; normally 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
print,'Early constant value of fit',cval,ycon
;
kcurve=where((allx ge ycon) and (allx le 1994),ncurve)
allwt=replicate(1.,ncurve)
acoeff=[0.,1./7200.]
if matchvar eq 0 then funcadd='_regress' else funcadd='_matchvar'
vval=curvefit(allx(kcurve),ally(kcurve),allwt,acoeff,$
  function_name='funct_decline'+funcadd)
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)
;
endif
;
if doadjust gt 1 then begin
;
;!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,$
  c_colors=cc,$
  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
;
!p.multi[0]=0
pause
inter_boxfd,zcoeff,g.x,g.y,$
  labels=labels,map=map,$
  c_colors=cc,$
  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,$
  c_colors=cc,$
  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 a 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'+fnadd+'.idlsave',$
  g,mxdyear,mxdnyr,fdcalib,mxdfd2,fdcorrect
;
endif
;
end