source: trunk/serietempECM.pro @ 611

;+
; _serietempECM.pro:
;
; ================
; serietempECM.pro
; ================
;
; DESCRIPTION
; ===========
;
; Analyse fichiers complets 2000 - 2007 ECMWF
; Attention : 2 and de ERA40 et 6 ans d'OPERA
; Hovmoller et series temporelles
;
; EVOLUTIONS
; ==========
;
; $Id$
;
; $URL$
; 
; - fplod 20120413T092359Z cratos (Linux)
; 
;   * add minimal header
;
;-

; Hovmoller 2005 - 2007
;;;;;;;;;;;;;;;;;;;;;;;;;;;
chp='r'
lo1=-10 & lo2=5 & la1=-30 & la2=40
dom=strcompress(string(lo1,lo2,la1,la2),/remove_all)

;geopotentiel : delta entre niveau 700 (7) et 925 hPa
lev1=700
;lecture premier fichier (OPERA)
file1 = '/volumes/temp/CAT.'+chp+'.2002.2007.6TU.nc'
print,file1
initncdf, file1,zaxis='level', glam = [-180,180]
print, gdept
domdef,lo1,lo2,la1,la2,lev1,lev1
data1 = read_ncdf(chp,20050101,20071231, file = file1)
tab1=data1.arr
lev2=850
;lev2=925
domdef,lo1,lo2,la1,la2,lev2,lev2
data1= read_ncdf(chp,20050101,20071231, file = file1)
tab2=data1.arr
;si geopot - sinon 
table1=(tab1+tab2)/2.
;table1=tab1-tab2
dada=smooth(table1,5) ;/100.
min=210 & max=245 & in=2
min=5 & max=90 & in=8
;pltt,window=0, dada, 'yt',title='latitude - time '+chp+'/100.',/rempli,min,max,int=in
pltt,window=0, dada, 'yt',title='latitude - time '+chp+'%',inv=1,/rempli,min,max,int=in

saveimage,'SORTIES/'+chp+'hovy2005-2007'+dom+'.png'

;series temporelles 
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
chp='geopt'
lo1=-5 & lo2=5 & la1=15 & la2=25
dom=strcompress(string(lo1,lo2,la1,la2),/remove_all)

;geopotentiel : delta entre niveau 700 (7) et 925 hPa
;domaine pour geopt: 20-25 N ?

lev1=700
;lecture premier fichier (OPERA)
file1 = '/volumes/temp/CAT.'+chp+'.2002.2007.6TU.nc'
print,file1
initncdf, file1,zaxis='level', glam = [-180,180]
print, gdept
domdef,lo1,lo2,la1,la2,lev1,lev1
data1 = read_ncdf(chp,20050101,20071231, file = file1)
tab1=data1.arr
;lev2=850
lev2=925
domdef,lo1,lo2,la1,la2,lev2,lev2
data1= read_ncdf(chp,20050101,20071231, file = file1)
tab2=data1.arr
;si geopot - sinon 
;table1=(tab1+tab2)/2.
;table1=tab1-tab2
ntime1=max(size(time))
y1=fltarr(ntime1)
serie1=reform(table1)
dimserie1=size(serie1)
nx=dimserie1(1)-1 & ny=dimserie1(2)-1
for i=0,ntime1-1 do y1(i)=mean(serie1(0:nx,0:ny,i))

; si travail sur les trois dernieres annees, on saute ce passage
; lecture deuxieme fichier (ERA)
chp2=chp ;chp2='z';!!!!!attention pour geopotentiel
file2 = '/volumes/temp/CAT.'+chp2+'.2000.2001.6TU.nc'
initncdf, file2,zaxis='level', glam = [-180,180]
print, gdept
domdef,lo1,lo2,la1,la2,lev1,lev1
data2 = read_ncdf(chp2,20000101,20011231, file = file2)
tab1=data2.arr
;domdef,lo1,lo2,la1,la2,lev2,lev2
;data2= read_ncdf(chp,20020101,20071231, file = file2)
;table2=data2.arr
;si geopot - sinon 
table2=tab1
;table2=tab1-tab2
ntime2=max(size(time))
y2=fltarr(ntime2)
serie2=reform(table2)
dimserie2=size(serie2)
nx=dimserie2(1)-1 & ny=dimserie2(2)-1
for i=0,ntime2-1 do y2(i)=mean(serie2(0:nx,0:ny,i))
;construction de la serie temporelle complete et du temps associe
jpt=ntime1+ntime2
print,ntime1,ntime2,jpt
time=timegen(jpt,start=julday(1,1,2000,6,0,0),step_size=1, units='Days')
data=[y2,y1]
y=smooth(data,5)

;si seulement OPERA:
jpt=ntime1
time=timegen(jpt,start=julday(1,1,2005,3,0,0),step_size=1, units='Days')
data=y1
y=data ;smooth(data,3)
pltt,y,'t',/rempli
dimserie=size(y)
; tentative d'application de wavelet
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
ntime=max(size(time))
dt=1
;y=fltarr(ntime)
;nx=dimserie(1)-1 & ny=dimserie(2)-1
;for i=0,ntime-1 do y(i)=mean(serie(0:nx,0:ny,i))
fil=digital_filter(0.04,1,50,80)  ; flow fraction de fnyquist, ici 1/2dt=1/2
;fil=digital_filter(0.02,1,50,120)  ; flow fraction de fnyquist, ici 1/2dt=1/2

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,YRANGE=[MAX(period),MIN(period)],/nodata,ystyle=1

CONTOUR,ABS(wave)^2,time,period, $
       XSTYLE=1,ystyle=1,YTITLE='Period',TITLE='Noise Wavelet '+chp+' 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/waveletECMfil(0.02)'+chp+dom+'.png',/png
saveimage,'sorties/waveletECM2005-2007'+chp+dom+'.png',/png