source: trunk/src/paper01/fig16/fig16_timeseries_nino3.pro @ 43

;+
; .. _fig16_timeseries_nino3.pro:
;
; ==========================
; fig16_timeseries_nino3.pro
; ==========================
;
; DESCRIPTION
; ===========
;
; SEE ALSO
; ========
;
; :ref:`tropflux_profile.sh`
;
; EXAMPLES
; ========
;
; ::
;
;  IDL> @tropflux_init
;  IDL> fig16_timeseries_nino3
;
; EVOLUTIONS
; ==========
;
; $Id$
;
; - fplod 20110411T142955Z aedon.locean-ipsl.upmc.fr (Darwin)
;
;   * minimal header
;
;-
pro fig16_timeseries_nino3
@common

reinitplt, /z,/invert
key_portrait = 1

openps, FILENAME = 'idl.ps'
date1=19940101 & date2=20071231
box=[240,330,-5,5]

file_trop='/Users/pkb/data/TropFlux/TropFlux_sst_19890101_20091231.nc'
initncdf, file_trop
sst=read_ncdf("sst", date1, date2, box=box, file=file_trop,/nostr)
sst=grossemoyenne(sst, "xy",/nan)
restore, filename="/Users/pkb/work/MY_SAXO/netflux_nino3_1989_2007.idl"
help, trop, oafl, nce2, era, nce, sst
;; creating climatology
ny=(2007-1994)+1
nt=365*ny
nsmooth=30.

trop_nclim=total(reform(trop(0:nt-1),365,ny),2)/ny
oafl_nclim=total(reform(oafl(0:nt-1),365,ny),2)/ny
nce2_nclim=total(reform(nce2(0:nt-1),365,ny),2)/ny
era_nclim=total(reform(era(0:nt-1),365,ny),2)/ny
nce_nclim=total(reform(nce(0:nt-1),365,ny),2)/ny
sst_nclim=total(reform(sst(0:nt-1),365,ny),2)/ny

trop_nclim=[trop_nclim ,trop_nclim, trop_nclim] & trop_nclim=smooth(trop_nclim,nsmooth,/nan) & trop_nclim=trop_nclim[365:365+364]
oafl_nclim=[oafl_nclim ,oafl_nclim, oafl_nclim] & oafl_nclim=smooth(oafl_nclim,nsmooth,/nan) & oafl_nclim=oafl_nclim[365:365+364]
nce2_nclim=[nce2_nclim ,nce2_nclim ,nce2_nclim] & nce2_nclim=smooth(nce2_nclim,nsmooth,/nan) & nce2_nclim=nce2_nclim[365:365+364]
era_nclim=[era_nclim ,era_nclim, era_nclim] & era_nclim=smooth(era_nclim,nsmooth,/nan) & era_nclim=era_nclim[365:365+364]
nce_nclim=[nce_nclim ,nce_nclim ,nce_nclim] & nce_nclim=smooth(nce_nclim,nsmooth,/nan) & nce_nclim=nce_nclim[365:365+364]
sst_nclim=[sst_nclim ,sst_nclim ,sst_nclim] & sst_nclim=smooth(sst_nclim,nsmooth,/nan) & sst_nclim=sst_nclim[365:365+364]

time=time & jpt=n_elements(time)

caldat, time,mon,day,yea
trop_clim=trop*0. & era_clim=trop*0.
oafl_clim=trop*0. & nce_clim=trop*0.
nce2_clim=trop*0. & sst_clim=trop*0.

for jt=0,jpt-1 do begin
  jtt=(time(jt)-julday(1,1,yea(jt))) < 364
  trop_clim(jt)=trop_nclim(jtt)
  oafl_clim(jt)=oafl_nclim(jtt)
  nce2_clim(jt)=nce2_nclim(jtt)
  era_clim(jt)=era_nclim(jtt)
  nce_clim(jt)=nce_nclim(jtt)
  sst_clim(jt)=sst_nclim(jtt)
endfor

trop_ano=trop-trop_clim & oafl_ano=oafl-oafl_clim
nce2_ano=nce2-nce2_clim & era_ano=era-era_clim
nce_ano=nce-nce_clim    & sst_ano=sst-sst_clim

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

Thf=1 & 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_ano)
     ff_ts=fft(ts)
     ff_ts_f=ff_ts*msk_spec
     trop_hf=float(fft(ff_ts_f,/inverse))
;
    ts=reform(nce2_ano)
    ff_ts=fft(ts)
    ff_ts_f=ff_ts*msk_spec
    nce2_hf=float(fft(ff_ts_f,/inverse))
;
    ts=reform(era_ano)
    ff_ts=fft(ts)
    ff_ts_f=ff_ts*msk_spec
    era_hf=float(fft(ff_ts_f,/inverse))
;
    ts=reform(nce_ano)
    ff_ts=fft(ts)
    ff_ts_f=ff_ts*msk_spec
    nce_hf=float(fft(ff_ts_f,/inverse))
;
    ts=reform(oafl_ano)
    ff_ts=fft(ts)
    ff_ts_f=ff_ts*msk_spec
    oafl_hf=float(fft(ff_ts_f,/inverse))
;
help, trop_hf, nce2_hf, era_hf, nce_hf

trop_ano=trop_ano-trop_hf & oafl_ano=oafl_ano-oafl_hf
nce2_ano=nce2_ano-nce2_hf & era_ano=era_ano-era_hf
nce_ano=nce_ano-nce_hf

nsmooth=30.
marge=[-3,-3,0,4]

pltt, ts_smooth(sst_ano,nsmooth,/nan), "t",/rempl,/nocolorb, $
      small=[1,3,1], title="a) Nino 3 index ", subtitle="", $
      ytitle="(degree Celcius)", marge=marge, linethick=1.3, charsize=1.1, xmin=-1, ymin=-1

nsmooth=90.

pltt, ts_smooth(trop,nsmooth,/nan), "t",/rempl,/nocolorb, 0,150, $
      small=[1,3,2], title="b) Net Heat Flux", subtitle="", color=30, $
      ytitle="(W/m2)", marge=marge, linethick=1.3, charsize=1.1,/noer, xmin=-1, ymin=-1
pltt, ts_smooth(oafl,nsmooth,/nan), "t",/ov1d, color=145
pltt, ts_smooth(era,nsmooth,/nan), "t",/ov1d, color=190
pltt, ts_smooth(nce2,nsmooth,/nan), "t",/ov1d, color=250
pltt, ts_smooth(nce,nsmooth,/nan), "t",/ov1d, color=65

nsmooth=120
pltt, ts_smooth(trop_ano,nsmooth,/nan), "t",/rempl,/nocolorb, -45,30, color=30, $
      small=[1,3,3], ytitle="(W/m2)", marge=marge, linethick=1.3, xmin=-1, ymin=-1, $
      title="c) Interannual Net Flux",/noer,  charsize=1.1, subtitle=""
pltt, ts_smooth(oafl_ano,nsmooth,/nan), "t",/ov1d, color=145
pltt, ts_smooth(era_ano,nsmooth,/nan), "t",/ov1d, color=190
pltt, ts_smooth(nce2_ano,nsmooth,/nan), "t",/ov1d, color=250
pltt, ts_smooth(nce_ano,nsmooth,/nan), "t",/ov1d, color=65

closeps
fig='fig16_timeseries_nino3.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