print, 'Running ecmwf_waveamp'

;
; Need to specify number of files, times per file, increment to step through times, and total number
; of times when stepping.  Assumes 4 times per day.
;
  nFiles = 9

  ecmwfFilesTZRW = ["/biptmp/joemci/ecmwf_toga/Dec2005/DSS.Y87610_2005335_2005344_ECMWF.nc", $
                    "/biptmp/joemci/ecmwf_toga/Dec2005/DSS.Y87611_2005345_2005354_ECMWF.nc", $
                    "/biptmp/joemci/ecmwf_toga/Dec2005/DSS.Y87612_2005355_2005365_ECMWF.nc", $
                    "/biptmp/joemci/ecmwf_toga/Dec2005/DSS.D5632_2006001_2006010_ECMWF.nc", $
                    "/biptmp/joemci/ecmwf_toga/Dec2005/DSS.D5634_2006011_2006020_ECMWF.nc", $
                    "/biptmp/joemci/ecmwf_toga/Dec2005/DSS.D5636_2006021_2006031_ECMWF.nc", $
                    "/biptmp/joemci/ecmwf_toga/Dec2005/DSS.D5638_2006032_2006041_ECMWF.nc", $
                    "/biptmp/joemci/ecmwf_toga/Dec2005/DSS.D5640_2006042_2006051_ECMWF.nc", $
                    "/biptmp/joemci/ecmwf_toga/Dec2005/DSS.D5642_2006052_2006059_ECMWF.nc"]
  ecmwfFilenameUV = "/ptmp/joemci/ECMWF/DSS.Y87886_2005355_2005365_ECMWF.nc"

  nFT = [40,40,44,40,40,44,40,40,32]
  nStep = 4
  nTAll = FIX(TOTAL(nFT)/nStep)
;
;  Read first ECMWF file to get dimensions of data
;
  ncdf_ecmwf_read,ecmwfFilesTZRW[0],ecmwfFilenameUV,lon,lat,lev,time,timeH,$
                                               nX,nY,nZ,nT,tN,zG,uN,vN,wN,rH,nameCnst

  tNAll = FLTARR(nX,nY,nZ,nTAll)
;  tN = tNTemp(*,*,*,0:nT-1)
;  zGTemp = zG
  zGAll = FLTARR(nX,nY,nZ,nTAll)
;  zG = zGTemp(*,*,*,0:nT-1)
;  uNTemp = uN
  uNAll = FLTARR(nX,nY,nZ,nTAll)
;  uN = uNTemp(*,*,*,0:nT-1)
;  vNTemp = vN
  vNAll = FLTARR(nX,nY,nZ,nTAll)
;  vN = vNTemp(*,*,*,0:nT-1)
;  wNTemp = wN
  wNAll = FLTARR(nX,nY,nZ,nTAll)
;  wN = wNTemp(*,*,*,0:nT-1)
;  rHTemp = rH
  rHAll = FLTARR(nX,nY,nZ,nTAll)
;  rH = rHTemp(*,*,*,0:nT-1)  

  nTOut = 0
  FOR iFE = 0, nFiles-1 DO BEGIN

    ecmwfFilenameTZRW = ecmwfFilesTZRW[iFE]

    print, 'Reading files ', iFE, ecmwfFilenameTZRW, ecmwfFilenameUV

    ncdf_ecmwf_read,ecmwfFilenameTZRW,ecmwfFilenameUV,lon,lat,lev,time,timeH,$
                                               nX,nY,nZ,nT,tN,zG,uN,vN,wN,rH,nameCnst

;    nTIn = 0
;    FOR iT = nTAll,nTAll+nFT[iFE]-1,nStep DO BEGIN
    FOR iT = 0,nFT[iFE]-1,nStep DO BEGIN

      tNAll[*,*,*,nTOut] = tN[*,*,*,iT] 
      uNAll[*,*,*,nTOut] = uN[*,*,*,iT] 
      vNAll[*,*,*,nTOut] = vN[*,*,*,iT] 
      zGAll[*,*,*,nTOut] = zG[*,*,*,iT] 
;      nTIn = nTIn + 1
      nTOut = nTOut + 1

    ENDFOR

;    nTAll = nTAll + nFT[iFE]
    
  ENDFOR 
;
;  Need to reverse latitude array for ECMWF to match WACCM
;
  latTemp = lat
  FOR iLat=0,nY-1 DO lat[iLat] = latTemp[nY-iLat-1]
  
  nP = 6
  Pres = [1.0,5.0,10.0,50.0,100.0,500.0]  

  tNSub = FLTARR(nX,nY,nP,nTOut)
  uNSub = FLTARR(nX,nY,nP,nTOut)
  vNSub = FLTARR(nX,nY,nP,nTOut)
  zGSub = FLTARR(nX,nY,nP,nTOut)

  FOR iP=0,nP-1 DO BEGIN
;
;  For ECMWF data need to find data for selected vertical levels using 1D vertical array 
;
    PInd = WHERE(lev GE Pres[iP])

    IF PInd[0] EQ -1 THEN MESSAGE, 'Cannot find pressure level index in ECMWF', iP, Pres[iP]
;
;  For ECMWF data need to flip data in latitude dimension 
;  
    FOR iLt=0,nY-1 DO BEGIN
;
;  Need to flip ECMWF data in latitude to match WACCM
;
      iRevLt = nY-iLt-1
     
      tNSub(*,iLt,iP,*) = tNAll(*,iRevLt,PInd[0],*)
      uNSub(*,iLt,iP,*) = uNAll(*,iRevLt,PInd[0],*)
      vNSub(*,iLt,iP,*) = vNAll(*,iRevLt,PInd[0],*)
      zGSub(*,iLt,iP,*) = zGAll(*,iRevLt,PInd[0],*)

    ENDFOR

  ENDFOR  

;lon,lat,level,time zAmp
 
  zAmp    = FLTARR(nX,nY,nP,nTOut)
  tAmp    = FLTARR(nX,nY,nP,nTOut)
  
  zLonMean    = FLTARR(nY,nP,nTOut)
  zLatMean    = FLTARR(nP,nTOut)
;
;  Find pressure levels to process in WACCM data
;
  FOR iT=0,nTOut-1 DO BEGIN

    FOR iLt=0,nY-1 DO BEGIN
;
;  For WACCM data need to find data for selected vertical levels using 3D vertical array
;
      FOR iP=0,nP-1 DO BEGIN

        zAmp[*,iLt,iP,iT] = ABS(FFT(zGSub[*,iLt,iP,iT]/1000.0, -1))
        tAmp[*,iLt,iP,iT] = ABS(FFT(tNSub[*,iLt,iP,iT], -1))
	zLonMean[iLt,iP,iT] = MEAN(zGSub[*,iLt,iP,iT])
	
      ENDFOR

    ENDFOR

    FOR iPP=0,nP-1 DO zLatMean[iPP,iT] = MEAN(zLonMean[*,iPP,iT])

  ENDFOR

TopTitleA = 'FFT Amplitude of Geopotential Height Lat=10 Level=10,Time='
PLOT, /ylog, zAmp[0,*,*,0],xrange=[0,10],xtitle=XTitle,ytitle=YTitle,Title=TopTitleA
; DEVICE, DECOMPOSED = 0
; LOADCT, 39
; contour, zamp[1:10,*,20,0],findgen(10),latw,nlev=140,c_colors=BYTSCL(indgen(140)),/fill,yrange=[-90,90],ystyle=1
;contour, tamp[1:10,*,20,0],findgen(10),latw,nlev=140,c_colors=BYTSCL(indgen(140)),/fill,yrange=[-90,90],ystyle=1


set_plot, 'ps'
device, /color, filename='idl_tamp.ps'
loadct, 39
!p.multi=[2,2]
fgen90 = findgen(90)
fgen40 = findgen(40)
xTitle = 'Time(3hr)'
yTitle = 'Frequency'

Title = ' Lat=-87.5 Level=20 km FFT Amp. Temperature'
contour, reform(transpose(tamp[1:40,1,4,*])),fgen90,fgen40,nlev=140,c_colors=BYTSCL(indgen(140)),/fill,yrange=[0,10],ystyle=1,title=Title, xtitle=xTitle, ytitle=yTitle
Title = ' Lat=-80 Level=20 km'
contour, reform(transpose(tamp[1:40,4,4,*])),fgen90,fgen40,nlev=140,c_colors=BYTSCL(indgen(140)),/fill,yrange=[0,10],ystyle=1,title=Title, xtitle=xTitle, ytitle=yTitle
Title = ' Lat=-70 Level=20 km FFT Amp. Temperature'
contour, reform(transpose(tamp[1:40,8,4,*])),fgen90,fgen40,nlev=140,c_colors=BYTSCL(indgen(140)),/fill,yrange=[0,10],ystyle=1,title=Title, xtitle=xTitle, ytitle=yTitle
Title = ' Lat=-60 Level=20 km'
contour, reform(transpose(tamp[1:40,12,4,*])),fgen90,fgen40,nlev=140,c_colors=BYTSCL(indgen(140)),/fill,yrange=[0,10],ystyle=1,title=Title, xtitle=xTitle, ytitle=yTitle
Title = ' Lat=-50 Level=20 km FFT Amp. Temperature'
contour, reform(transpose(tamp[1:40,16,4,*])),fgen90,fgen40,nlev=140,c_colors=BYTSCL(indgen(140)),/fill,yrange=[0,10],ystyle=1,title=Title, xtitle=xTitle, ytitle=yTitle
Title = ' Lat=-40 Level=20 km'
contour, reform(transpose(tamp[1:40,20,4,*])),fgen90,fgen40,nlev=140,c_colors=BYTSCL(indgen(140)),/fill,yrange=[0,10],ystyle=1,title=Title, xtitle=xTitle, ytitle=yTitle
Title = ' Lat=-30 Level=20 km FFT Amp. Temperature'
contour, reform(transpose(tamp[1:40,24,4,*])),fgen90,fgen40,nlev=140,c_colors=BYTSCL(indgen(140)),/fill,yrange=[0,10],ystyle=1,title=Title, xtitle=xTitle, ytitle=yTitle
Title = ' Lat=-20 Level=20 km'
contour, reform(transpose(tamp[1:40,28,4,*])),fgen90,fgen40,nlev=140,c_colors=BYTSCL(indgen(140)),/fill,yrange=[0,10],ystyle=1,title=Title, xtitle=xTitle, ytitle=yTitle
Title = ' Lat=-10 Level=20 km FFT Amp. Temperature'
contour, reform(transpose(tamp[1:40,32,4,*])),fgen90,fgen40,nlev=140,c_colors=BYTSCL(indgen(140)),/fill,yrange=[0,10],ystyle=1,title=Title, xtitle=xTitle, ytitle=yTitle
Title = ' Lat=0 Level=20 km'
contour, reform(transpose(tamp[1:40,36,4,*])),fgen90,fgen40,nlev=140,c_colors=BYTSCL(indgen(140)),/fill,yrange=[0,10],ystyle=1,title=Title, xtitle=xTitle, ytitle=yTitle
Title = ' Lat=10 Level=20 km FFT Amp. Temperature'
contour, reform(transpose(tamp[1:40,40,4,*])),fgen90,fgen40,nlev=140,c_colors=BYTSCL(indgen(140)),/fill,yrange=[0,10],ystyle=1,title=Title, xtitle=xTitle, ytitle=yTitle
Title = ' Lat=20 Level=20 km'
contour, reform(transpose(tamp[1:40,44,4,*])),fgen90,fgen40,nlev=140,c_colors=BYTSCL(indgen(140)),/fill,yrange=[0,10],ystyle=1,title=Title, xtitle=xTitle, ytitle=yTitle
Title = ' Lat=30 Level=20 km FFT Amp. Temperature'
contour, reform(transpose(tamp[1:40,48,4,*])),fgen90,fgen40,nlev=140,c_colors=BYTSCL(indgen(140)),/fill,yrange=[0,10],ystyle=1,title=Title, xtitle=xTitle, ytitle=yTitle
Title = ' Lat=40 Level=20 km'
contour, reform(transpose(tamp[1:40,52,4,*])),fgen90,fgen40,nlev=140,c_colors=BYTSCL(indgen(140)),/fill,yrange=[0,10],ystyle=1,title=Title, xtitle=xTitle, ytitle=yTitle
Title = ' Lat=50 Level=20 km FFT Amp. Temperature'
contour, reform(transpose(tamp[1:40,56,4,*])),fgen90,fgen40,nlev=140,c_colors=BYTSCL(indgen(140)),/fill,yrange=[0,10],ystyle=1,title=Title, xtitle=xTitle, ytitle=yTitle
Title = ' Lat=60 Level=20 km'
contour, reform(transpose(tamp[1:40,60,4,*])),fgen90,fgen40,nlev=140,c_colors=BYTSCL(indgen(140)),/fill,yrange=[0,10],ystyle=1,title=Title, xtitle=xTitle, ytitle=yTitle
Title = ' Lat=70 Level=20 km FFT Amp. Temperature'
contour, reform(transpose(tamp[1:40,64,4,*])),fgen90,fgen40,nlev=140,c_colors=BYTSCL(indgen(140)),/fill,yrange=[0,10],ystyle=1,title=Title, xtitle=xTitle, ytitle=yTitle
Title = ' Lat=80 Level=20 km'
contour, reform(transpose(tamp[1:40,68,4,*])),fgen90,fgen40,nlev=140,c_colors=BYTSCL(indgen(140)),/fill,yrange=[0,10],ystyle=1,title=Title, xtitle=xTitle, ytitle=yTitle
Title = ' Lat=87.5 Level=20 km FFT Amp. Temperature'
contour, reform(transpose(tamp[1:40,71,4,*])),fgen90,fgen40,nlev=140,c_colors=BYTSCL(indgen(140)),/fill,yrange=[0,10],ystyle=1,title=Title, xtitle=xTitle, ytitle=yTitle

device,/close

device, /color, filename='idl_zamp.ps'
loadct, 39
!p.multi=[2,2]
fgen90 = findgen(90)
fgen40 = findgen(40)
xTitle = 'Time(3hr)'
yTitle = 'Frequency'

Title = ' Lat=-87.5 Level=20 km FFT Amp. Geopot. Ht.'
contour, reform(transpose(zamp[1:40,1,4,*])),fgen90,fgen40,nlev=140,c_colors=BYTSCL(indgen(140)),/fill,yrange=[0,10],ystyle=1,title=Title, xtitle=xTitle, ytitle=yTitle
Title = ' Lat=-80 Level=20 km'
contour, reform(transpose(zamp[1:40,4,4,*])),fgen90,fgen40,nlev=140,c_colors=BYTSCL(indgen(140)),/fill,yrange=[0,10],ystyle=1,title=Title, xtitle=xTitle, ytitle=yTitle
Title = ' Lat=-70 Level=20 km FFT Amp. Geopot. Ht.'
contour, reform(transpose(zamp[1:40,8,4,*])),fgen90,fgen40,nlev=140,c_colors=BYTSCL(indgen(140)),/fill,yrange=[0,10],ystyle=1,title=Title, xtitle=xTitle, ytitle=yTitle
Title = ' Lat=-60 Level=20 km'
contour, reform(transpose(zamp[1:40,12,4,*])),fgen90,fgen40,nlev=140,c_colors=BYTSCL(indgen(140)),/fill,yrange=[0,10],ystyle=1,title=Title, xtitle=xTitle, ytitle=yTitle
Title = ' Lat=-50 Level=20 km FFT Amp. Geopot. Ht.'
contour, reform(transpose(zamp[1:40,16,4,*])),fgen90,fgen40,nlev=140,c_colors=BYTSCL(indgen(140)),/fill,yrange=[0,10],ystyle=1,title=Title, xtitle=xTitle, ytitle=yTitle
Title = ' Lat=-40 Level=20 km'
contour, reform(transpose(zamp[1:40,20,4,*])),fgen90,fgen40,nlev=140,c_colors=BYTSCL(indgen(140)),/fill,yrange=[0,10],ystyle=1,title=Title, xtitle=xTitle, ytitle=yTitle
Title = ' Lat=-30 Level=20 km FFT Amp. Geopot. Ht.'
contour, reform(transpose(zamp[1:40,24,4,*])),fgen90,fgen40,nlev=140,c_colors=BYTSCL(indgen(140)),/fill,yrange=[0,10],ystyle=1,title=Title, xtitle=xTitle, ytitle=yTitle
Title = ' Lat=-20 Level=20 km'
contour, reform(transpose(zamp[1:40,28,4,*])),fgen90,fgen40,nlev=140,c_colors=BYTSCL(indgen(140)),/fill,yrange=[0,10],ystyle=1,title=Title, xtitle=xTitle, ytitle=yTitle
Title = ' Lat=-10 Level=20 km FFT Amp. Geopot. Ht.'
contour, reform(transpose(zamp[1:40,32,4,*])),fgen90,fgen40,nlev=140,c_colors=BYTSCL(indgen(140)),/fill,yrange=[0,10],ystyle=1,title=Title, xtitle=xTitle, ytitle=yTitle
Title = ' Lat=0 Level=20 km'
contour, reform(transpose(zamp[1:40,36,4,*])),fgen90,fgen40,nlev=140,c_colors=BYTSCL(indgen(140)),/fill,yrange=[0,10],ystyle=1,title=Title, xtitle=xTitle, ytitle=yTitle
Title = ' Lat=10 Level=20 km FFT Amp. Geopot. Ht.'
contour, reform(transpose(zamp[1:40,40,4,*])),fgen90,fgen40,nlev=140,c_colors=BYTSCL(indgen(140)),/fill,yrange=[0,10],ystyle=1,title=Title, xtitle=xTitle, ytitle=yTitle
Title = ' Lat=20 Level=20 km'
contour, reform(transpose(zamp[1:40,44,4,*])),fgen90,fgen40,nlev=140,c_colors=BYTSCL(indgen(140)),/fill,yrange=[0,10],ystyle=1,title=Title, xtitle=xTitle, ytitle=yTitle
Title = ' Lat=30 Level=20 km FFT Amp. Geopot. Ht.'
contour, reform(transpose(zamp[1:40,48,4,*])),fgen90,fgen40,nlev=140,c_colors=BYTSCL(indgen(140)),/fill,yrange=[0,10],ystyle=1,title=Title, xtitle=xTitle, ytitle=yTitle
Title = ' Lat=40 Level=20 km'
contour, reform(transpose(zamp[1:40,52,4,*])),fgen90,fgen40,nlev=140,c_colors=BYTSCL(indgen(140)),/fill,yrange=[0,10],ystyle=1,title=Title, xtitle=xTitle, ytitle=yTitle
Title = ' Lat=50 Level=20 km FFT Amp. Geopot. Ht.'
contour, reform(transpose(zamp[1:40,56,4,*])),fgen90,fgen40,nlev=140,c_colors=BYTSCL(indgen(140)),/fill,yrange=[0,10],ystyle=1,title=Title, xtitle=xTitle, ytitle=yTitle
Title = ' Lat=60 Level=20 km'
contour, reform(transpose(zamp[1:40,60,4,*])),fgen90,fgen40,nlev=140,c_colors=BYTSCL(indgen(140)),/fill,yrange=[0,10],ystyle=1,title=Title, xtitle=xTitle, ytitle=yTitle
Title = ' Lat=70 Level=20 km FFT Amp. Geopot. Ht.'
contour, reform(transpose(zamp[1:40,64,4,*])),fgen90,fgen40,nlev=140,c_colors=BYTSCL(indgen(140)),/fill,yrange=[0,10],ystyle=1,title=Title, xtitle=xTitle, ytitle=yTitle
Title = ' Lat=80 Level=20 km'
contour, reform(transpose(zamp[1:40,68,4,*])),fgen90,fgen40,nlev=140,c_colors=BYTSCL(indgen(140)),/fill,yrange=[0,10],ystyle=1,title=Title, xtitle=xTitle, ytitle=yTitle
Title = ' Lat=87.5 Level=20 km FFT Amp. Geopot. Ht.'
contour, reform(transpose(zamp[1:40,71,4,*])),fgen90,fgen40,nlev=140,c_colors=BYTSCL(indgen(140)),/fill,yrange=[0,10],ystyle=1,title=Title, xtitle=xTitle, ytitle=yTitle

device,/close
set_plot,'x'

;  nTAllW2 = nTAllW/2+1
;  nShiftW = (nXW+1)/2-1

;  zonalfreq,uuNSubW,nTAllW,nXW,nYW,nP,nShiftW,nTAllW2,uuSpcW,uuDTW,uuSDTW,uuP10W,uuP20W,uuP12W,uuP21W
;  zonalfreq,vuNSubW,nTAllW,nXW,nYW,nP,nShiftW,nTAllW2,vuSpcW,vuDTW,vuSDTW,vuP10W,vuP20W,vuP12W,vuP21W
;  zonalfreq,tuNSubW,nTAllW,nXW,nYW,nP,nShiftW,nTAllW2,tuSpcW,tuDTW,tuSDTW,tuP10W,tuP20W,tuP12W,tuP21W

;  xwn = indgen(nx)-nshift
;  freq = indgen(nT2)
;  xwnw = nshiftw-indgen(nxw)
;  freqw = indgen(nTAllW2)
; DEVICE, DECOMPOSED = 0
; LOADCT, 39
; nlevels = 20
; clevels = fltarr(nlevels)
; fmin = min(alog10(uuspc[*,*,36,2]))
; fmax = max(alog10(uuspc[*,*,36,2]))
; conint = (fmax-fmin) / (nlevels-1)
; for i=0,nlevels-1 do clevels(i) = fmin+i*conint
; clevelsw = fltarr(nlevels)
; fminw = min(alog10(uuspcw[*,*,36,2]))
; fmaxw = max(alog10(uuspcw[*,*,36,2]))
; conintw = (fmaxw-fminw) / (nlevels-1)
; for i=0,nlevels-1 do clevelsw(i) = fminw+i*conintw
; window, ret=2
;contour, alog10(uuspc[*,*,36,0]), xwn,freq, nlev=21, /cell_fill, xrange=[-10,10], yrange=[0,20]
;contour, alog10(uuspc[*,*,36,2]), xwn,freq, nlev=21, /cell_fill, xrange=[-10,10], yrange=[0,10]
;contour, alog10(uuspcw[*,*,36,2]), xwnw,freqw, nlev=21, /cell_fill, xrange=[-10,10], yrange=[0,6]
;plot, lat, uuspc[70,11,*,0]
;plot, lat, uuspc[70,11,*,1]
;plot, lat, uuspc[70,11,*,2]
;contour, alog10(uuspc[*,*,36,2]), xwn,freq,/noerase,/follow,/fill,xstyle=5,ystyle=5,levels=clevels,c_colors=BYTSCL(indgen(20)),xrange=[-10,10], yrange=[0,10]
;contour, alog10(uuspc[*,*,36,2]), xwn,freq,/noerase,/follow,levels=clevels,xstyle=1, xrange=[-10,10], yrange=[0,10]
;contour, alog10(uuspcw[*,*,36,2]), xwnw,freqw,/noerase,/follow,/fill,xstyle=5,ystyle=5,levels=clevelsw,c_colors=BYTSCL(indgen(20)),xrange=[-10,10], yrange=[0,10]
;contour, alog10(uuspcw[*,*,36,2]), xwnw,freqw,/noerase,/follow,levels=clevelsw,xstyle=1, xrange=[-10,10], yrange=[0,10]
;
;set_plot, 'ps'
;device, /color
;loadct, 39
;!p.multi=[2,2]
;contour, alog10(uuspc[*,*,36,0]), xwn,freq, nlev=21, /cell_fill, xrange=[-10,10],Title='ECMWF 1 hPA'
;contour, alog10(uuspcw[*,*,47,0]), xwnw,freqw, nlev=21, /cell_fill, xrange=[-10,10], yrange=[0,20],Title='WACCM 1 hPA'
;contour, alog10(uuspc[*,*,36,1]), xwn,freq, nlev=21, /cell_fill, xrange=[-10,10],Title='ECMWF 10 hPA'
;contour, alog10(uuspcw[*,*,47,1]), xwnw,freqw, nlev=21, /cell_fill, xrange=[-10,10], yrange=[0,20],Title='WACCM 10 hPA'
;contour, alog10(uuspc[*,*,36,2]), xwn,freq, nlev=21, /cell_fill, xrange=[-10,10],Title='ECMWF 100 hPA'
;contour, alog10(uuspcw[*,*,47,2]), xwnw,freqw, nlev=21, /cell_fill, xrange=[-10,10], yrange=[0,20],Title='WACCM 10 hPA'
;device,/close


END