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

!****s* WaveOpticsPropagator/ropp_pp_wopt_mps_2d *
!
! NAME
!    ropp_pp_wopt_mps_2d
!
! SYNOPSIS
!    CALL ropp_pp_wopt_mps_2d(x, y, radius, y_apodize, kval, dtheta, &
!                             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 (2D) refractivity slice.
!
!   For details, see Karayel and Hinson (1997), Radio Science, vol 32, no 2, 411-423
!
! 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    :: dtheta     ! angular separation between profiles in 2D plane (radians)
!    REAL    :: zed        ! geometric height of model levels (2D)
!    REAL    :: log_refrac ! log(refrac) on model levels (2D)
!    REAL    :: d2logndy2  ! d2logN/dy2 used in spline (2D)
!    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_2d(x, y, radius, y_apodize, kval, dtheta, &
                               zed, log_refrac, d2logndy2, U0, U)

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

  IMPLICIT NONE

! Input
  REAL(wp), INTENT(IN)     :: x(:), y(:)
  REAL(wp), INTENT(IN)     :: radius, y_apodize, kval, dtheta
  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, kk
  INTEGER                  :: nx, ny
  INTEGER                  :: msg_temp

  REAL(wp)                 :: dx, dy
  REAL(wp)                 :: refrac_2d(SIZE(y))   ! refractivity using 2D plane information 
  REAL(wp)                 :: theta, theta_mid, theta_min, theta_max  ! position in plane
  REAL(wp)                 :: wt(2) ! linear spatial interpolation coefficients
  REAL(wp)                 :: radius_window ! locally interpolated value of radius + zed
  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)

! angles in 2D plane

  theta_mid = SIZE(log_refrac, 2) * dtheta / 2.0_wp

  theta_min = -theta_mid

  theta_max =  theta_mid

  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_2d(dtheta, x(i+1), y, radius, zed, log_refrac, d2logndy2, refrac_2d)

! now update the phases

    DO j = 1, ny

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

    ENDDO

! Windowing more complicated for 2D than 1D. Now requires horizontal interpolation,
! because the input refractivities are on model levels, whose heights depend on
! surface pressure (~orography) and therefore vary from profile to profile in 
! the defining 2D refractivity profile.
! Compute radius value appropriate for x(i+1) in 2D windowing case
!
    theta = ASIN(x(i+1)/radius)

    kk = INT((theta+theta_mid)/dtheta) + 1   ! The nearest profile in the 2D refractivity slice

    kk = MAX( MIN(kk, SIZE(zed,2)-1), 1 )    ! Make sure kk is within limits

    IF ( theta < theta_max .AND. theta > theta_min) THEN

      wt(1) = (REAL(kk, wp)*dtheta-theta_mid-theta)/dtheta
      wt(2) = 1.0_wp - wt(1)

    ELSE IF (theta <= theta_min) THEN

      wt(1) = 1.0_wp
      wt(2) = 0.0_wp

    ELSE IF (theta >= theta_max) THEN

      wt(1) = 0.0_wp
      wt(2) = 1.0_wp

    ENDIF

    radius_window = radius + wt(1)*zed(1,kk) + wt(2)*zed(1,kk+1)

    CALL ropp_pp_wopt_window(x(i+1), y, radius_window, y_apodize, U(:))

  ENDDO

  DEALLOCATE (U_tmp)

END SUBROUTINE ropp_pp_wopt_mps_2d