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

!****s* WaveOpticsPropagator/ropp_pp_wopt_mps_1d *
!
! NAME
!    ropp_pp_wopt_mps_1d
!
! SYNOPSIS
!    CALL ropp_pp_wopt_mps_1d(x, y, radius, y_apodize, kval, &
!                             zed, log_refrac, d2logndy2, U0, U)
!
! DESCRIPTION
!   Calculates complex signal (phase and amplitude) on final screen in a set,
!   given the signal on the first screen and the (1D) refractivity profiles
!   to be applied throughout.
!
! INPUTS
!    REAL    :: x,y        ! positions of x,y grid
!    REAL    :: radius     ! radius of curvature
!    REAL    :: y_apodise  ! apodisation height to damp signal above this height  
!    REAL    :: kval       ! wave number of signal
!    REAL    :: zed        ! geometric height of model levels (1D)
!    REAL    :: log_refrac ! log(refrac) on model levels (1D)
!    REAL    :: d2logndy2  ! d2logN/dy2 used in spline (1D)
!    COMPLEX :: U0         ! Complex amplitude at first screen
!
! OUTPUT
!    COMPLEX :: U          ! Complex amplitude at final screen
!
! NOTES
!   See RSR 28 for details
!
! REFERENCES
!   Karayel and Hinson (1997), Radio Science, vol 32, no 2, 411-423
!   ROM SAF Report 28, SAF/ROM/ECMWF/REP/RSR/028
!
! AUTHOR
!   ECMWF, Reading, UK.
!   Any comments on this software should be given via the ROM SAF
!   Helpdesk at http://www.romsaf.org
!
! COPYRIGHT
!   (c) EUMETSAT. All rights reserved.
!   For further details please refer to the file COPYRIGHT
!   which you should have received as part of this distribution.
!
!****

SUBROUTINE ropp_pp_wopt_mps_1d(x, y, radius, y_apodize, kval, &
                               zed, log_refrac, d2logndy2, U0, U)

  USE typesizes, ONLY: wp => EightByteReal
  USE messages
  USE ropp_pp_wopt, not_this => ropp_pp_wopt_mps_1d
  USE ropp_pp_constants, ONLY: pi1 => pi

  IMPLICIT NONE

! Input
  REAL(wp), INTENT(IN)     :: x(:), y(:)
  REAL(wp), INTENT(IN)     :: radius, y_apodize, kval
  REAL(wp), INTENT(IN)     :: zed(:), log_refrac(:), d2logndy2(:)
  COMPLEX(wp), INTENT(IN)  :: U0(:)      ! Signal at first screen

! Output
  COMPLEX(wp), INTENT(OUT) :: U(SIZE(y)) ! Signal at final screen

! Local
  INTEGER                  :: i, j
  INTEGER                  :: nx, ny
  INTEGER                  :: msg_temp

  REAL(wp)                 :: dx, dy
  REAL(wp)                 :: refrac(SIZE(y))   ! refractivity on profile
  REAL(wp)                 :: twopi

  COMPLEX(wp), PARAMETER   :: Ci = (0.0_wp, 1.0_wp)  ! Sqrt(-1)
  COMPLEX(wp), ALLOCATABLE :: U_tmp(:)

  CHARACTER (LEN=3)        :: si1, si2, snx  ! For diagnostic messaging

! initialise

  nx   = SIZE(x)
  ny   = SIZE(y)

  twopi = 2.0_wp * pi1

! initialise the field

  U(:) = U0(:)

! loop through the screens

  dy = y(2) - y(1)  ! this should be fixed for the entire screen

  ALLOCATE (U_tmp(ny))

  DO i = 1, nx-1

    WRITE (si1, FMT='(I3)') i
    WRITE (si2, FMT='(I3)') i + 1
    WRITE (snx, FMT='(I3)') nx
    IF ( ((i/10)*10 == i) .AND. (msg_mode /= VerboseMode) ) THEN
      WRITE (si2, FMT='(I3)') MIN(i+10, nx)
      msg_temp = msg_info
    ELSE
      WRITE (si2, FMT='(I3)') i + 1
      msg_temp = msg_diag
    END IF
    CALL message( msg_temp, 'Propagating signal between screens ' // &
                            si1 // ' and ' // si2 // ' (out of ' // snx // ')' )

    U_tmp(:) = U(:)

! step in horizontal

    dx = x(i+1) - x(i)

    CALL ropp_pp_wopt_propagate(U_tmp, kval, dy, dx, U)

! compute the refractivity on the i+1 plane

    CALL ropp_pp_wopt_ref_1d(x(i+1), y, radius, zed, log_refrac, d2logndy2, refrac)

! now update the phases

    DO j = 1, ny

       U(j) = U(j) * EXP( ci * MODULO(kval*dx*refrac(j), twopi) )

    ENDDO

! the window function

    CALL ropp_pp_wopt_window(x(i+1), y, radius+zed(1), y_apodize, U)

  ENDDO

  DEALLOCATE (U_tmp)

END SUBROUTINE ropp_pp_wopt_mps_1d