pro zonalfreq0,un,ntime,nx,ny,nz,nshift,ntime2,uspc,uphs
;; Frequency-zonal wavenumber analyzer at a given time, altitude and latitude.
;; Hanli Liu, 7/99

uspc = fltarr(nx,ntime2,ny,nz)
uphs = fltarr(nx,ntime2,ny,nz)
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
;            tmpa(*,it) = shift(tmpa(*,it),-nshift) ;;this is to put 0
;            deg at the first point  (for TIME-GCM only)
        endfor
        fftmp = fft(tmpa)

        for it = 0,ntime2-1 do begin
            uspc(*,it,j,k) = shift(abs(fftmp(*,it)),nshift)*2.
;            uphs(*,it,j,k) = -shift(atan(fftmp(*,it),/phase),nshift)*12./!pi
            uphs(*,it,j,k) = -shift(atan(fftmp(*,it),/phase),nshift)
        endfor
        for it = 1,ntime2-1 do begin
            uphs(*,it,j,k) = uphs(*,it,j,k)/it
        endfor
     endfor
endfor
;uphs(where(uphs lt 0))=uphs(where(uphs lt 0))+24.

uphs(where(uphs lt 0)) = uphs(where(uphs lt 0))+2.*!pi
for iwvn = 2,nshift do begin
   lonprd = 2.*!pi/iwvn
   wvnn = nshift+iwvn
   wvnp = nshift-iwvn
   uphs(wvnp,*,*,*) = uphs(wvnp,*,*,*) mod (lonprd)
   uphs(wvnn,*,*,*) = uphs(wvnn,*,*,*) mod (lonprd)
endfor
uspc(*,0,*,*)=uspc(*,0,*,*)/2.

return
end