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

!****s* WaveOpticsPropagator/ropp_pp_wopt_ref_1d *
!
! NAME
!    ropp_pp_wopt_ref_1d
!
! SYNOPSIS
!    CALL ropp_pp_wopt_ref_1d(x, y, radius, zed, logn, d2logndy2, refrac)
!
! DESCRIPTION
!   Reconstructs the refractivity from its spline fitting coefficients.
!
! INPUTS
!   REAL    :: x,y         ! positions of x, y grid
!   REAL    :: radius      ! radius of curvature
!   REAL    :: zed, logn   ! input z, logN
!   REAL    :: d2logndy2   ! input d2logN/dy2 (from spline)
!
! OUTPUT
!   REAL    :: refrac      ! fitted refractivity
!
! 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) 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_ref_1d(x, y, radius, zed, logn, d2logndy2, refrac)

  USE typesizes, ONLY: wp => EightByteReal

  IMPLICIT NONE

! Input/output
  REAL(wp), INTENT(IN)  :: x                 ! x-position of screen
  REAL(wp), INTENT(IN)  :: y(:)              ! y-position of points on screen
  REAL(wp), INTENT(IN)  :: radius            ! radius of curvature (+ undulation)
  REAL(wp), INTENT(IN)  :: zed(:), logn(:)   ! 1D z, logN
  REAL(wp), INTENT(IN)  :: d2logndy2(:)      ! d2logN/dy2
  REAL(wp), INTENT(OUT) :: refrac(SIZE(y))   ! output y

! Local
  INTEGER               :: ny, nlev
  INTEGER               :: i, ipos
  REAL(wp)              :: cos_theta, hval, grad_top
  REAL(wp)              :: h, a, b           ! for spline fit

  ny = SIZE(y)

  nlev = SIZE(zed)

  refrac(:) = 0.0_wp

  grad_top = (logn(nlev)-logn(nlev-1))/(zed(nlev)-zed(nlev-1))

  ipos = 1

  DO i = 1, ny

! the theta value will change with y

    cos_theta = COS( ATAN(x/(radius + y(i))) )

    hval = (y(i) - radius*(cos_theta - 1.0_wp))/cos_theta

    IF (hval < zed(1)) THEN  ! constant extrapolation below the surface
                             ! to mitigate the effect of a sharp boundary
                             ! at the surface.  Subsequent Gaussian
                             ! windowing/damping should make the results
                             ! insensitive to the way this is done.

      refrac(i) = 1.0E-6_wp*EXP(logn(1))

    ELSE IF (hval >= zed(nlev)) THEN   ! linear extrapolation of log N above the model top

      refrac(i) = logn(nlev) + grad_top*(hval-zed(nlev))

      refrac(i) = MAX(1.0E-6_wp*EXP(refrac(i)), 1.0E-18_wp)

    ELSE  ! spline interpolation in between

      DO

        IF ( zed(ipos+1) > hval .OR. (ipos+1) == nlev ) EXIT

        ipos = ipos + 1

      ENDDO

      ipos = MIN( MAX(1, ipos), nlev-1 )

      h = zed(ipos+1) - zed(ipos)

      a = (zed(ipos+1) - hval) / h

      b = 1.0_wp - a

! compute LOG(refrac) from the spline coefficients

      refrac(i) = a*logn(ipos) + b*logn(ipos+1) + &
        ((a**3-a)*d2logndy2(ipos)+(b**3-b)*d2logndy2(ipos+1))*(h**2)/6.0_wp

      refrac(i) = 1.0E-6_wp*EXP(refrac(i))

    ENDIF

  ENDDO

END SUBROUTINE ropp_pp_wopt_ref_1d