pro zonalfreq1,un,ntime,nx,ny,nz,nshift,ntime2,uspc,nwv,nfrq,ucomp
;; Frequency-zonal wavenumber analyzer at a given time, altitude and latitude.
;; Also reconstruct wavenumber nwv and frequency number nfrq (negative for
;; westward propagating).
;; Hanli Liu, 4/02

;; uspc is the zonal wavenumber-frequency spectrum(amplitude) of un,
;; with first dimension as zonal wavenumber and second dimension as
;; frequency number. Note that the original point (wavenumber 0 and
;; frequency number 0) is with index nshift for wavenumber dimension
;; and 0 for frequency dimension. All components with wave number
;; index (0:nshift-1) are EASTWARD propagating, while (nshift+1:nx-1)
;; are WESTWARD propagating.

uspc = fltarr(nx,ntime2,ny,nz)

;Specified component
ucomp = fltarr(nx,ny,nz,ntime)

;Mask matrices that set components other than the desired one to zero

umask = fltarr(nx,ntime)

if (nfrq lt 0) then begin
   umask(nwv,-nfrq) = 1.
   umask(nx-nwv,ntime+nfrq) = 1.
endif
if (nfrq eq 0) then begin
   umask(nwv,0) = 1.
   umask(nx-nwv,0) = 1.
endif
if (nfrq gt 0) then begin
   umask(nwv,ntime-nfrq) = 1.
   umask(nx-nwv,nfrq) = 1.
endif

for k = 0,nz-1 do begin
    for j = 0,ny-1 do begin
        tmpa = fltarr(nx,ntime)
        for it = 0,ntime-1 do begin
            for i=0,nx-1 do begin
                tmpa(i,it)=un(i,j,k,it)
            endfor
        endfor
        fftmp = fft(tmpa)

        for it = 0,ntime2-1 do begin
            uspc(*,it,j,k) = shift(abs(fftmp(*,it)),nshift)*2.
        endfor

        fftmp1 = complexarr(nx,ntime)
; Specified component
        fftmp1 = fftmp*umask
        rftmp = fft(fftmp1,1)
        for it = 0,ntime-1 do begin
            for i=0,nx-1 do begin
                ucomp(i,j,k,it) = float(rftmp(i,it))
            endfor
        endfor

     endfor
endfor
uspc(*,0,*,*)=uspc(*,0,*,*)/2.

end