;
; Plots time series of the difference between
; temperature and density averaged over the region north of 50N.
; Uses both corrected and uncorrected data.
; 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.
;
matchvar=1        ; 0=regression, 1=variance matching
;
print,'Reading MXD and temperature data'
if matchvar eq 0 then fnadd='_regress' else fnadd=''
restore,filename='calibmxd5'+fnadd+'.idlsave'
; Gets:  g,mxdyear,mxdnyr,fdcalibu,fdcalibc,mxdfd2,timey,fdseas
;
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=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
def_1color,10,color='red'
def_1color,11,color='blue'
def_1color,12,color='green'
def_1color,13,color='orange'
;
; 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
;
fddiffu=fdcalibu(*,*,kmxd)-fdseas(*,*,ktem)
fddiffc=fdcalibc(*,*,kmxd)-fdseas(*,*,ktem)
;
; bothmask is the time-dependent mask where both data sets have values
bothmask=float(finite(fddiffu))    ; 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)          ; usually ge 50!
fd=fddiffu(*,ky,*)
difftsu=globalmean(fd,g.y(ky))
fd=fddiffc(*,ky,*)
difftsc=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=fdcalibu(*,ky,*)
mxdtsu=globalmean(fd,g.y(ky))
fd=fdcalibc(*,ky,*)
mxdtsc=globalmean(fd,g.y(ky))
;
; (4) Difference of the mean series
diffts2u=mxdtsu(kmxd)-temts(ktem)
diffts2c=mxdtsc(kmxd)-temts(ktem)
;
; (5) Mean of MXD (explicit time-varying MXD&T mask)
fd=fdcalibu(*,ky,*)
fd=fd[*,*,kmxd]
maskmxdtsu=globalmean(fd,g.y[ky],mask=bothmask[*,ky,*])
fd=fdcalibc(*,ky,*)
fd=fd[*,*,kmxd]
maskmxdtsc=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
diffts3u=maskmxdtsu-masktemts
diffts3c=maskmxdtsc-masktemts
;
thalf=10
filter_cru,thalf,/nan,tsin=difftsu,tslow=difflowu
filter_cru,thalf,/nan,tsin=difftsc,tslow=difflowc
filter_cru,thalf,/nan,tsin=temts,tslow=temlow
filter_cru,thalf,/nan,tsin=mxdtsu,tslow=mxdlowu
filter_cru,thalf,/nan,tsin=mxdtsc,tslow=mxdlowc
filter_cru,thalf,/nan,tsin=maskmxdtsu,tslow=maskmxdlowu
filter_cru,thalf,/nan,tsin=maskmxdtsc,tslow=maskmxdlowc
filter_cru,thalf,/nan,tsin=masktemts,tslow=masktemlow
filter_cru,thalf,/nan,tsin=diffts2u,tslow=difflow2u
filter_cru,thalf,/nan,tsin=diffts2c,tslow=difflow2c
filter_cru,thalf,/nan,tsin=diffts3u,tslow=difflow3u
filter_cru,thalf,/nan,tsin=diffts3c,tslow=difflow3c
;
plot,mxdyear,mxdlowu,/nodata,$
  /xstyle,xrange=[1850,2000],xtitle='Year',$
  /ystyle,yrange=[-0.9,0.9],$
  ytitle='>50N masked temperature (!Uo!NC wrt 1961-90)'
;oplot,mxdyear,mxdlowu,thick=2,color=10
;oplot,mxdyear,mxdlowc,thick=2,color=11
;oplot,timey,temlow,thick=4
oplot,mxdyear[kmxd],maskmxdlowu,thick=2
oplot,mxdyear[kmxd],maskmxdlowc,thick=2,color=11
oplot,timey[ktem],masktemlow,thick=4,color=10
oplot,!x.crange,[0.,0.],linestyle=1
;
legend,thick=[4,2,2,2,2],color=[10,!p.color,11,13,12],$
  ['>50N temperature (masked)','>50N MXD reconstruction (masked)',$
    '>50N MXD reconstruction (corrected then masked)','>50N MXD reconstruction (masked then scaled)',$
    '>50N MXD reconstruction (corrected then masked then scaled']
;
; 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(maskmxdtsc,masktemts)
bestcoeff=linfit(maskmxdtsu[kcalib],masktemts[kcalib])
;r=correlate(maskmxdtsc,masktemts)
r=correlate(maskmxdtsu[kcalib],masktemts[kcalib])
print,'Area-average later, then skill   (masked) =',r
print,'Best-fit scaling coefficients=',bestcoeff
oplot,mxdyear[kmxd],bestcoeff[0]+bestcoeff[1]*maskmxdlowu,thick=2,color=13
oplot,mxdyear[kmxd],bestcoeff[0]+bestcoeff[1]*maskmxdlowc,thick=2,color=12
;
plot,timey(ktem),difflow2u,/nodata,$
  /xstyle,xrange=[1850,2000],xtitle='Year',$
  /ystyle,yrange=[-1.3,0.4],ytitle='Temperautre difference (!Uo!NC)'
oplot,timey(ktem),difflow2u,thick=2
oplot,timey(ktem),difflowu,thick=4
oplot,timey(ktem),difflow2c,thick=2,color=11
oplot,timey(ktem),difflowc,thick=4,color=11
oplot,!x.crange,[0.,0.],linestyle=1
legend,color=[!p.color,!p.color,11,11],thick=[2,4,2,4],$
  ['Mean of differences','Difference of means',$
  'Mean of differences (corrected)','Difference of means (corrected)']
;
end