;+
; .. _fig16_timeseries_nino3.pro:
;
; ==========================
; fig16_timeseries_nino3.pro
; ==========================
;
; DESCRIPTION
; ===========
;
; :file:`${PROJECT_OD}/TropFlux_sst_19890101_20091231.nc`
; containing
; sst corrected on OAFLUX grid
; has been produced by
; :ref:`sst_correction_ncdf.pro`.
;
; :file:`${PROJECT_OD}/netflux_nino3_1989_2007.idl`
; containing
; ++
; has been produced by
; :ref:`++`.

; .. graphviz::
;
;    digraph fig16_timeseries_nino3 {
;
;       tropflux_sst [shape=ellipse,fontname=Courier,label="${PROJECT_OD}/TropFlux_sst_19890101_20091231.nc"];
;       figure [shape=ellipse,fontname=Courier,label="${PROJECT_OD}/WHOTS_validation_net_2004_2007_all_v10.ps"];
;
;       fig16_timeseries_nino3 [shape=box,
;       fontname=Courier,
;       color=blue,
;       URL="http://forge.ipsl.jussieu.fr/tropflux/broswrer/trunk/src/paper01/fig16/fig16_timeseries_nino3.pro",
;       label="${TROPFLUX}/src/paper01/fig16/fig16_timeseries_nino3.pro"];
;
;       {tropflux_sst} -> {fig16_timeseries_nino3} -> {figure}
;    }
;
; SEE ALSO
; ========
;
; :ref:`project_profile.sh`
; :ref:`project_init.pro`
; :ref:`cm_project.pro`
;
; :func:`x_site_location`
; :func:`y_site_location`
;
; EXAMPLES
; ========
;
; ::
;
;  IDL> fig16_timeseries_nino3
;
; TODO
; ====
;
; make it work : missing :file:`${PROJECT_OD}/netflux_nino3_1989_2007.idl`
;
; coding rules
;
; get rid of restore
;
; complete description
;
; handle IO error
;
; EVOLUTIONS
; ==========
;
; $Id$
;
; $URL$
;
; - fplod 20110817T102957Z aedon.locean-ipsl.upmc.fr (Darwin)
;
;   * add PROJECT_OD to restored file
;   * replace PROJECT_ID by PROJECT_OD
;
; - fplod 20110420T123110Z aedon.locean-ipsl.upmc.fr (Darwin)
;
;   * remove hard coding path
;   * add graphviz
;   * externalize functions
;
; - fplod 20110411T142955Z aedon.locean-ipsl.upmc.fr (Darwin)
;
;   * minimal header
;
;-
pro fig16_timeseries_nino3
;
@cm_4cal
@cm_4data
@cm_4mesh
@cm_4data
@cm_project

reinitplt, /z,/invert
key_portrait = 1

openps, FILENAME = project_od_env+'fig16_timeseries_nino3.ps'
date1=19940101 
date2=20071231
box=[240,330,-5,5]

file_trop=project_od_env+'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=project_od_env+'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
end