
      subroutine mkampha(f,amp,phase,imx,kmx,jmx,ut,spv)
      dimension f(imx,kmx,jmx),amp(jmx,kmx,4),phase(jmx,kmx,4)
      dimension w(200),x(imx),field(imx,jmx,kmx)
c
      do i=1,imx
        do j=1,jmx
          field(i,j,:) = f(i,:,j)
        enddo
      enddo
      call ampha(4,imx,imx,jmx,kmx,field,amp,phase,x,w,ut,spv) 
      return
      end
c
      subroutine mkampha0(f,amp,phase,imx,kmx,jmx,ut,spv)
      dimension f(imx,jmx,kmx),amp(jmx,kmx,4),phase(jmx,kmx,4)
      dimension w(200),x(imx)
      call ampha(4,imx,imx,jmx,kmx,f,amp,phase,x,w,ut,spv) 
      return
      end
C
      SUBROUTINE AMPHA(N,IMAXP4,IMAX,JMAX,KMAX,P,AMP,PHASE,X,W,TEE,spv)
      SAVE                                                              
C     ****     CALCULATES AMPLITUDE AND PHASE OF N LOWEST ZONAL FOURIER 
C     ****     COMPONENTS FOR ARRAY P(IMAX,JMAX,KMAX). AMPLITUDE        
C     ****     RETURNED IN AMP(JMAX,KMAX,N), PHASE IN PHASE(JMAX,KMAX,N)
      DIMENSION P(IMAXP4,JMAX,1),AMP(JMAX,KMAX,1),PHASE(JMAX,KMAX,1),X(1
     A),W(1)                                                            
      DATA PI/3.1415926535/                                             
C     ****     INITIALIZE FOURIER TRANSFORM                             
      CALL RFFTI(IMAX,W)                                                
      DO 100 K=1,KMAX                                                     
      DO 100 J=1,JMAX                                                     
      DO 2 I=1,IMAX                                                     
        X(I)=P(I,J,K)                                                     
    2 CONTINUE                                                          
      iok = 1
      do i=1,imax
        if (x(i).eq.spv) iok = 0  
      enddo
      if (iok.gt.0) then
        CALL RFFTF(IMAX,X,W)                                              
        DO 1 L=1,N                                                        
          M=L+1                                                             
          AMP(J,K,L)=SQRT(X(2*M)**2+X(2*M-1)**2)/FLOAT(IMAX)                
C     ****                                                              
C     ****     LOCAL PHASE                                              
C     ****                                                              
          TSHIFT = FLOAT(L)*PI*TEE/12.                                      
          COST = COS(TSHIFT)                                                
          SINT = SIN(TSHIFT)                                                
          AT = X(2*M)*COST+X(2*M-1)*SINT                                    
          BT = X(2*M-1)*COST-X(2*M)*SINT                                    
          IF(BT.EQ.0.)GO TO 3                                               
          PHASE(J,K,L) = ATAN2(AT,BT)                                       
          GO TO 1                                                           
    3     CONTINUE                                                          
          IF(AT.EQ.0.)GO TO 4                                               
          PHASE(J,K,L) = 2.*ATAN(AT/ABS(AT))                                
          GO TO 1                                                           
    4     CONTINUE                                                          
          PHASE(J,K,L) = 0.                                                 
    1   CONTINUE                                                          
      else
        do l=1,n
          amp(j,k,l) = spv
          phase(j,k,l) = spv
        enddo
      endif
 100  continue
      RETURN                                                            
      END                                                               
C                                                                       

