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

!****s* WaveOpticsPropagator/ropp_pp_wopt_ref_2d *
!
! NAME
!    ropp_pp_wopt_ref_2d
!
! SYNOPSIS
!    CALL ropp_pp_wopt_ref_2d(dtheta, x, y, radius, zed, logn, d2logndy2, refrac)
!
! DESCRIPTION
!   Reconstructs the refractivity from its spline fitting coefficients.
!
! INPUTS
!   REAL    :: dtheta      ! angular separation
!   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_2d(dtheta, x, y, radius, zed, logn, d2logndy2, refrac)

  USE typesizes, ONLY: wp => EightByteReal

  IMPLICIT NONE

! Input/output
  REAL(wp), INTENT(IN)  :: dtheta              ! angular separation of profiles
  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(:,:) ! 2D z, logN
  REAL(wp), INTENT(IN)  :: d2logndy2(:,:)      ! d2logN/dy2
  REAL(wp), INTENT(OUT) :: refrac(SIZE(y))     ! output y

! Local
  INTEGER               :: ny, nlev, nhoriz, kk, k
  INTEGER               :: i, ipos
  REAL(wp)              :: cos_theta, hval, grad_top
  REAL(wp)              :: h, a, b             ! for spline fit
  REAL(wp)              :: theta,theta_mid,theta_min,theta_max,wt(2)
  REAL(wp)              :: log_refrac

  ny = SIZE(y)

  nlev   = SIZE(zed,1)

  nhoriz = SIZE(zed,2)

  refrac(:) = 0.0_wp

  theta_mid = nhoriz * dtheta / 2.0_wp

  theta_min = -theta_mid

  theta_max =  theta_mid

  ipos = 1

  DO i = 1, ny

! the theta value will change with y

    theta = ATAN(x/(radius + y(i)))

    kk = INT((theta+theta_mid)/dtheta) + 1

    kk = MAX(MIN(kk,nhoriz-1),1)

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

      wt(1) = ((REAL(kk)*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

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

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

    refrac(i) = 0.0_wp  ! initialise

    DO k = 1,2 ! loop over the interpolation coefficients

! the theta value will change with y

      IF (hval < zed(1,kk)) 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) = refrac(i)+wt(k)*1.0E-6_wp*EXP(logn(1,kk))

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

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

        log_refrac = logn(nlev,kk) + grad_top*(hval-zed(nlev,kk))

        refrac(i) = refrac(i)+ wt(k)*MAX(1.0E-6_wp*EXP(log_refrac), 1.0E-18_wp)

      ELSE  ! spline interpolation in between

        ipos = 1  ! should be ascending so store previous

        DO

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

          ipos = ipos + 1

        ENDDO

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

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

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

        b = 1.0_wp - a

! compute LOG(refrac) from the spline coefficients

        log_refrac = a*logn(ipos,kk) + b*logn(ipos+1,kk) + &
          ((a**3-a)*d2logndy2(ipos,kk)+(b**3-b)*d2logndy2(ipos+1,kk))*(h**2)/6.0_wp

        refrac(i) = refrac(i)+ wt(k)*1.0E-6_wp*EXP(log_refrac)

      ENDIF

        kk = kk + 1

    ENDDO ! k = 1, 2

  ENDDO ! i=1, ny

END SUBROUTINE ropp_pp_wopt_ref_2d