;
; 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.
;
print,'Reading MXD and temperature data'
restore,filename='calibmxd5.idlsave'
; Gets:  g,mxdyear,mxdnyr,fdcalibu,fdcalibc,mxdfd2,timey,fdseas
;
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
def_1color,10,color='red'
def_1color,11,color='blue'
;
; 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)
;
; Compute means over the region north of 50N
;
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))
print,'Temperature time series (masked)'
fd=fdseas(*,ky,*)
temts=globalmean(fd,g.y(ky),mask=fdmask(*,ky))
print,'MXD time series'
fd=fdcalibu(*,ky,*)
mxdtsu=globalmean(fd,g.y(ky))
fd=fdcalibc(*,ky,*)
mxdtsc=globalmean(fd,g.y(ky))
diffts2u=mxdtsu(kmxd)-temts(ktem)
diffts2c=mxdtsc(kmxd)-temts(ktem)
;
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=diffts2u,tslow=difflow2u
filter_cru,thalf,/nan,tsin=diffts2c,tslow=difflow2c
;
plot,mxdyear,mxdlowu,/nodata,$
  /xstyle,xrange=[1400,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,!x.crange,[0.,0.],linestyle=1
legend,thick=[4,2,2],color=[!p.color,10,11],$
  ['>50N temperature','MXD reconstruction','MXD reconstruction (corrected)']
;
plot,timey(ktem),difflow2u,/nodata,$
  /xstyle,xrange=[1400,2000],xtitle='Year',$
  /ystyle,yrange=[-1.1,0.6],ytitle='Temperautre difference (!Uo!NC)'
oplot,timey(ktem),difflow2u,thick=2,color=10
oplot,timey(ktem),difflowu,thick=4,color=10
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=[10,10,11,11],thick=[2,4,2,4],$
  ['Mean of differences','Difference of means',$
  'Mean of differences (corrected)','Difference of means (corrected)']
;
end