! $Id: ropp_pp_wopt_phase_and_amplitude.f90 r4887 2016-05-25 15:48:28Z sti $

!****s* WaveOpticsPropagator/ropp_pp_wopt_phase_and_amplitude *
!
! NAME
!   ropp_pp_wopt_phase_and_amplitude
!
! SYNOPSIS
!   CALL ropp_pp_wopt_phase_and_amplitude(U, phase, amplitude)
!
! DESCRIPTION
!   Decompose complex signal into amplitude and accumulated phase
!
! INPUTS
!   COMPLEX(wp), INTENT(IN)  :: U  ! Input signal
!
! OUTPUT
!   REAL(wp), INTENT(OUT)    :: phase, amplitude
!
! NOTES
!   See RSR 28 for details
!
! AUTHOR
!   ECMWF, Reading, UK.
!   Any comments on this software should be given via the ROM SAF
!   Helpdesk at http://www.romsaf.org
!
! COPYRIGHT
!   (c) Copyright 1998, M. E. Gorbunov.
!   For further details please refer to the file COPYRIGHT
!   which you should have received as part of this distribution.
!
!****

SUBROUTINE ropp_pp_wopt_phase_and_amplitude(U, phase, amplitude)

  USE typesizes, ONLY: wp => EightByteReal
  USE ropp_pp_constants, ONLY: pi1 => pi

  IMPLICIT NONE

! Input/output
  COMPLEX(wp), INTENT(IN)  :: U(:)
  REAL(wp),    INTENT(OUT) :: phase(SIZE(U)), amplitude(SIZE(U))

! Local
  INTEGER                  :: i, ny
  REAL(wp)                 :: twopi
  REAL(wp), ALLOCATABLE    :: phase_gorb(:)

! initialise

  ny = SIZE(U)

  ALLOCATE (phase_gorb(ny))

  twopi = 2.0_wp * pi1

  phase(:)     = 0.0_wp
  amplitude(:) = 0.0_wp

! amplitude

  amplitude(:) = ABS(U(:))

! instantaneous phase, where calculable (zero elsewhere)

  WHERE ( amplitude > TINY(0.0_wp) ) phase = ATAN2(AIMAG(U), REAL(U))

! accumulated phase

  phase_gorb(:) = 0.0_wp

  DO i = 2,ny

    phase_gorb(i) = phase_gorb(i-1) + MODULO(phase(i)-phase_gorb(i-1)+pi1, twopi) - pi1

  ENDDO

  phase(:) = phase_gorb(:)

  DEALLOCATE (phase_gorb)

END SUBROUTINE ropp_pp_wopt_phase_and_amplitude