;+
; .. _serieAMSU.pro:
;
; =============
; serieAMSU.pro
; =============
;
; DESCRIPTION
; ============
;
; lecture fichiers complets AMSU pour tracer hovmöllers et ondelettes
;
; EVOLUTIONS
; ==========
;
; $Id$
;
; $URL$
;
; - fplod 20120413T091954Z cratos (Linux)
;
;   * minimal header
;
;-
;
;Hovmöller 2005 - 2007
;;;;;;;;;;;;;;;;;;;;;;

lo1=-10 & lo2=5 & la1=-30 & la2=40
dom=strcompress(string(lo1,lo2,la1,la2),/remove_all)
canal='a5'

j1=366+365+365+365+366
j2=j1+365+365+364
file='stagiaire/'+canal+'.2000_2007.nc'
initncdf, file, xaxisname = 'xlon', yaxisname = 'ylat'

domdef,lo1,lo2,la1,la2
data=read_ncdf('moyenne_tb', j1,j2,/timestep,timevar = 'jours', file = file)
time=julday(1,1,2005)+lindgen(j2-j1+1)
dada=smooth(data.arr,5)
min=260 & max=290 & in=3
min=245 & max=267 & in=1
pltt,window=0, dada, 'yt',title='latitude - time '+canal,/rempli,min,max,int=in
saveimage,'SORTIES/'+canal+'hovy2005-2007.'+dom+'.png'

; tentative d'application de wavelet
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
retall
lo1=-5 & lo2=5 & la1=5 & la2=15
dom=strcompress(string(lo1,lo2,la1,la2),/remove_all)
canal='b5'
;annee=2005
;an=strcompress(string(annee),/remove_all)
;imo=1   ; no du mois
;im1=9
j1=366+365+365+365+366
j2=j1+365+365+364
file='stagiaire/'+canal+'.2000_2007.nc'
initncdf, file, xaxisname = 'xlon', yaxisname = 'ylat'
;j1=0 & j2=364
domdef,lo1,lo2,la1,la2
data=read_ncdf('moyenne_tb', j1,j2,/timestep,timevar = 'jours', file = file)
time=julday(1,1,2005)+lindgen(j2-j1+1)
ntime=max(size(time))
serieb=data.arr
pltt,window=0,serieb, 't',/rempli
dimserie=size(serieb)
y=fltarr(ntime)
nx=dimserie(1)-1 & ny=dimserie(2)-1
for i=0,ntime-1 do y(i)=mean(serieb(0:nx,0:ny,i))

dt=1
; premier param freq de coupure exprimée comme fraction de FNyquist,
; flow fraction de fnyquist, ici 1/2dt=1/2
; freq visee nbjours/ntot, donc fraction=freq visee/ FNyquisrt
;dernier param nb points pour construire le filtre (2xperiode visee)
fil=digital_filter(0.05,1,50,80)  ; 40 jours environ
 fil=digital_filter(0.04,1,50,90)

yout=convol(y,fil)
;yout=y
wave = WAVELET(yout,dt,PERIOD=period,COI=coi,/PAD,SIGNIF=signif)
window,0
nscale = N_ELEMENTS(period)
LOADCT,39
plot,time,period,/ytype,XSTYLE=1,YTITLE='Period',TITLE='Noise Wavelet', $
       YRANGE=[MAX(period),MIN(period)],/nodata,ystyle=1
CONTOUR,ABS(wave)^2,time,period, $
       XSTYLE=1,ystyle=1,YTITLE='Period',TITLE='Noise Wavelet'+canal+' dom'+dom, $
       YRANGE=[MAX(period),MIN(period)], $   ;*** Large-->Small period
       /YTYPE, $                             ;*** make y-axis logarithmic
       NLEVELS=25,/FILL

signif = REBIN(TRANSPOSE(signif),ntime,nscale)
CONTOUR,ABS(wave)^2/signif,time,period, $
           /OVERPLOT,LEVEL=1.0,C_ANNOT='95%'
PLOTS,time,coi,NOCLIP=0   ;*** anything "below" this line is dubious

; le plus incroyable, c'est que ca marche!!!!!
saveimage,'sorties/wavelet-2005-2007'+canal+'dom'+dom+'.png',/png