;------------------------------------------------------------
pro fig15_timeseries_mjo
@common
;------------------------------------------------------------

reinitplt, /z,/invert
key_portrait = 1

openps, FILENAME = 'idl.ps'
;------------------------------------------------------------
date1=20071001 & date2=20080331
date1=19990101 & date2=20021231
box=[50,80,-12,-5]

file_trop='/Users/pkb/data/TropFlux/TropFlux_19890101_20091231.nc'
initncdf, file_trop
shf_trop=read_ncdf("shf", date1, date2, box=box, file=file_trop,/nostr)
lhf_trop=read_ncdf("lhf", date1, date2, box=box, file=file_trop,/nostr)
swr_trop=read_ncdf("swr", date1, date2, box=box, file=file_trop,/nostr)
lwr_trop=read_ncdf("lwr", date1, date2, box=box, file=file_trop,/nostr) ;; & ref_lw=-1*lwr_trop
trop=swr_trop+(lwr_trop+lhf_trop+shf_trop)
help, trop
trop=grossemoyenne(trop, "xy",/nan)

file_ncep='/Volumes/Iomega_HDD/flux_reconstruction/ncep2/NCEP2_flux_19890101_20090729.nc'
initncdf, file_ncep
shf=read_ncdf("shf", date1, date2, box=box, file=file_ncep,/nostr)
lhf=read_ncdf("lhf", date1, date2, box=box, file=file_ncep,/nostr)
lwr=read_ncdf("lwr", date1, date2, box=box, file=file_ncep,/nostr)
swr=read_ncdf("swr", date1, date2, box=box, file=file_ncep,/nostr)

nce2=swr-(lwr+lhf+shf)
nce2=grossemoyenne(nce2, "xy",/nan)
help, nce2

file_erai='/Volumes/Iomega_HDD/work/flux_reconstruction/gridded_data/erai_fluxes_19930101_20090801_TROP_oafluxgrid.nc'
initncdf, file_erai
shf=read_ncdf("shf", date1, date2, box=box, file=file_erai,/nostr)
lhf=read_ncdf("lhf", date1, date2, box=box, file=file_erai,/nostr)
lwr=read_ncdf("lwr", date1, date2, box=box, file=file_erai,/nostr)
swr=read_ncdf("swr", date1, date2, box=box, file=file_erai,/nostr) & swr=-1*swr
era=swr-(lwr+lhf+shf)
era=grossemoyenne(era, "xy",/nan)
help, era

file='/Volumes/Iomega_HDD/flux_reconstruction/ncep/fluxe_ncep1_19890101_20091231.nc'
initncdf, file
lwr=read_ncdf("lwr", date1, date2, box=box, file=file,/nostr)
swr=read_ncdf("swr", date1, date2, box=box, file=file,/nostr)
lhf=read_ncdf("lhf", date1, date2, box=box, file=file,/nostr)
shf=read_ncdf("shf", date1, date2, box=box, file=file,/nostr)

nce=swr+shf+lhf+lwr
nce=grossemoyenne(nce, "xy",/nan)
help, nce

file_oaf='/Volumes/Iomega_HDD/work/flux_reconstruction/OAFLX_GRID/shf_oafluxgrid_1985_2008.nc'
initncdf, file_oaf
shf=read_ncdf("shf", date1, date2, box=box, file=file_oaf,/nostr)
file_oaf='/Volumes/Iomega_HDD/work/flux_reconstruction/OAFLX_GRID/lhf_oafluxgrid_1985_2008.nc'
initncdf, file_oaf
lhf=read_ncdf("lhf", date1, date2, box=box, file=file_oaf,/nostr)
file_oaf='/Volumes/Iomega_HDD/work/flux_reconstruction/OAFLX_GRID/lwr_oafluxgrid_1985_2007.nc'
initncdf, file_oaf
lwr=read_ncdf("lwr", date1, date2, box=box, file=file_oaf,/nostr)
file_oaf='/Volumes/Iomega_HDD/work/flux_reconstruction/OAFLX_GRID/swr_oafluxgrid_1985_2007.nc'
initncdf, file_oaf
swr=read_ncdf("swr", date1, date2, box=box, file=file_oaf,/nostr) & swr_trop=swr
oafl=swr-(lwr+lhf+shf)
help, oafl
oafl=grossemoyenne(oafl, "xy",/nan)

file="/Users/pkb/data/heat_budget/NOAA_OLR/olr_oafluxgrid_30n30s_19890101_20091231.nc"
initncdf, file
olr=read_ncdf("olr", date1-1, date2, box=box, file=file,/nostr)
olr=grossemoyenne(olr, "xy",/nan)

help, olr

save, olr,trop, nce2, era, nce, filename="netflux_sctr_1999_2002.idl"
;------------------------------------------------------------------------------------------------------------------------

trop_hf=trop*0. & nce2_hf=trop*0.
era_hf=trop*0. & nce_hf=trop*0.

Thf=30 & Tlf=90
tsamp=1.   ; monthly
freq=findgen(jpt/2+1)/(jpt*tsamp)
freq=[freq,-1.*reverse(freq(1:jpt/2))]
peri=1./(abs(freq) > 1.e-20)*(2*(freq ge 0)-1.)
; Variance
     msk_spec=float((abs(peri) ge Thf) and (abs(peri) le Tlf))
     ts=reform(trop)
     ff_ts=fft(ts)
     ff_ts_f=ff_ts*msk_spec
     trop_hf=float(fft(ff_ts_f,/inverse))
;
    ts=reform(nce2)
    ff_ts=fft(ts)
    ff_ts_f=ff_ts*msk_spec
    nce2_hf=float(fft(ff_ts_f,/inverse))
;
    ts=reform(era)
    ff_ts=fft(ts)
    ff_ts_f=ff_ts*msk_spec
    era_hf=float(fft(ff_ts_f,/inverse))
;
    ts=reform(nce)
    ff_ts=fft(ts)
    ff_ts_f=ff_ts*msk_spec
    nce_hf=float(fft(ff_ts_f,/inverse))
;

help, trop_hf, nce2_hf, era_hf, nce_hf
nsmooth=15.

pltt, ts_smooth(olr,nsmooth,/nan), "t",/rempl, small=[1,3,1], $
   charsize=1.1, title="a) Outgoing Longwave Radiation", subtitle=""

pltt, ts_smooth(trop,nsmooth,/nan), "t",/rempl, -70,130, $
      small=[1,3,2],/noer, charsize=1.1, title="b) Net Heat Flux", subtitle=""
pltt, ts_smooth(era,nsmooth,/nan), "t",/ov1d, color=250
pltt, ts_smooth(nce2,nsmooth,/nan), "t",/ov1d, color=150
pltt, ts_smooth(nce,nsmooth,/nan), "t",/ov1d, color=50

pltt, trop_hf, "t",/rempl, -80,80, $
      small=[1,3,3],/noer, title="c) 30-90 day Net Heat Flux", charsize=1.1, subtitle=""
pltt, era_hf, "t",/ov1d, color=250
pltt, nce2_hf, "t",/ov1d, color=150
pltt, nce_hf, "t",/ov1d, color=50

;----------------------------------------------------------
closeps
fig='fig15_timeseries_mjo.ps'
spawn, 'mv '+psdir+'idl.ps '+cpsdir+fig
spawn, 'gv '+cpsdir+fig
return
end
;--------------------------------------------------------------------------
function x_site_location, site
    n1=strpos(site, 's')
if (n1 gt -1) then begin
    ns=-1.
    x=strmid(site, 0, n1)
    x=float(x)*ns
endif else begin
    n1=strpos(site, 'n')
    x=strmid(site, 0, n1)
    ny=1.
    x=float(x)*ny
endelse
return, float(x)
end
;--------------------------------------------------------------------------
function y_site_location, site
    n1=strpos(site, 'e')
if (n1 gt -1) then begin
    n=strpos(site, 's')
    if (n gt -1) then begin
        y=strmid(site, n+1, n1-n-1)
    endif else begin
        n=strpos(site, 'n')
        y=strmid(site, n+1, n1-n-1)
    endelse

endif else begin
    n1=strpos(site, 'w')
    n=strpos(site, 's')
    if (n gt -1) then begin
        y=strmid(site, n+1, n1-n-1)
        y=180+(180-float(y))
    endif else begin
        n=strpos(site, 'n')
        y=strmid(site, n+1, n1-n-1)
        y=180+(180-float(y))
   endelse
endelse
return,float(y)
end

;--------------------------------------------------------------------------