!
      module ampha
      contains
!-------------------------------------------------------------------
      subroutine defamp(amp,id1,id2,id3,plt,nf,spv)
!
! Args:
      integer,intent(in) :: id1,id2,id3,nf
      real,intent(in) :: amp(id1,id2,id3),spv
      real,intent(out) :: plt(id1,id2)
!
! Locals:
      integer :: k,j
!
      do k=1,id2
        if (amp(1,k,nf).ne.spv) then
          plt(1,k) = amp(1,k,nf)
        else
          plt(1,k) = spv
        endif
        do j=2,id1
          if (amp(j,k,nf).ne.spv.and.amp(j-1,k,nf).ne.spv) then
            plt(j,k) = .5*(amp(j,k,nf)+amp(j-1,k,nf))
          else
            plt(j,k) = spv
          endif
        enddo
      enddo
      return
      end subroutine defamp
!-------------------------------------------------------------------
      subroutine defphase(phase,id1,id2,id3,plt,scale,nf,spv)
      parameter (pi = 3.14159)
      dimension phase(id1,id2,id3),plt(id1,id2)
 
      do k=1,id2
        if (phase(1,k,nf).ne.spv) then
          plt(1,k) = phase(1,k,nf)*scale+12.
        else
          plt(1,k) = spv
        endif
        do j=2,id1
          if (phase(j,k,nf).ne.spv.and.phase(j-1,k,nf).ne.spv) then
            plt(j,k) = .5*(phase(j,k,nf)+phase(j-1,k,nf))*scale+12.
            if (abs(phase(j,k,nf)-phase(j-1,k,nf)).gt.pi)
     +        plt(j,k) = plt(j,k)+pi*scale
          else
            plt(j,k) = spv
          endif
        enddo
      enddo
      return
      end subroutine defphase
!-------------------------------------------------------------------
      subroutine getampha(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 subroutine getampha
!-------------------------------------------------------------------
      subroutine fixphase(phase,id1,id2,id3,n2,spv)
      logical,intent(in) :: n2
      dimension phase(id1,id2,id3)
      do k=1,id2
        if (phase(2,k,2).ne.spv) then
          phase(1,k,2) = phase(2,k,2)
        else
          phase(1,k,2) = spv
        endif
        if (phase(3,k,3).ne.spv) then
          phase(1,k,3) = phase(3,k,3)
        else
          phase(1,k,3) = spv
        endif
        if (phase(3,k,3).ne.spv) then
          phase(2,k,3) = phase(3,k,3)
        else 
          phase(2,k,3) = spv
        endif
        if (phase(4,k,4).ne.spv) then
          phase(1,k,4) = phase(4,k,4)
        else
          phase(1,k,4) = spv
        endif
        if (phase(4,k,4).ne.spv) then
          phase(2,k,4) = phase(4,k,4)
        else
          phase(2,k,4) = spv
        endif
        if (phase(4,k,4).ne.spv) then
          phase(3,k,4) = phase(4,k,4)
        else
          phase(3,k,4) = spv
        endif
        if (phase(id1-1,k,2).ne.spv) then
          phase(id1,k,2) = phase(id1-1,k,2)
        else
          phase(id1,k,2) = spv
        endif
        if (phase(id1-2,k,3).ne.spv) then
          phase(id1,k,3) = phase(id1-2,k,3)
        else
          phase(id1,k,3) = spv
        endif
        if (phase(id1-2,k,3).ne.spv) then
          phase(id1-1,k,3) = phase(id1-2,k,3)
        else
          phase(id1-1,k,3) = spv
        endif
        if (phase(id1-3,k,4).ne.spv) then
          phase(id1,k,4) = phase(id1-3,k,4)
        else
          phase(id1,k,4) = spv
        endif
        if (phase(id1-3,k,4).ne.spv) then
          phase(id1-1,k,4) = phase(id1-3,k,4)
        else
          phase(id1-1,k,4) = spv
        endif
        if (phase(id1-3,k,4).ne.spv) then
          phase(id1-2,k,4) = phase(id1-3,k,4)
        else
          phase(id1-2,k,4) = spv
        endif
      enddo
      if (.not.n2) then
        do nf=1,id3
          do j=1,id1
            if (phase(j,2,nf).ne.spv) then
              phase(j,1,nf) = phase(j,2,nf)
            else
              phase(j,1,nf) = spv
            endif
          enddo 
        enddo
      endif
      return
      end subroutine fixphase
      end module ampha
