! $Id: ropp_fm_iono_bangle.f90 4010 2014-01-10 11:07:40Z idculv $

!****s* Ionosphere/ropp_fm_iono_bangle *
!
! NAME
!   ropp_fm_iono_bangle - Forward model of ionospheric bending angle.
!
! SYNOPSIS
!   CALL ropp_fm_iono_bangle(Ne_max, R_peak, H_width, R_leo, &
!                            n_L1, impact, bangle)
!
! DESCRIPTION
!   Forward model of ionospheric bending for the L1 and L2 bending angles
!   assuming Chapman layer profiles. Should only be called if
!   x%direct_ion=.TRUE.
!
! INPUTS
!   REAL(wp)               :: Ne_max     ! Peak electron density (m-3)
!   REAL(wp)               :: R_peak     ! Radius of peak value (m)
!   REAL(wp)               :: H_width    ! Width of Chapman Layer (m)
!   REAL(wp)               :: R_leo      ! Radius of LEO (m)
!   INTEGER                :: n_L1       ! No. of L1 bending angles in array
!   REAL(wp), DIMENSION(:) :: impact     ! Impact parameters
!
! INOUTPUTS
!   REAL(wp), DIMENSION(:) :: bangle     ! Total bending angle inc. ionospheric component
!
! NOTES
!  Correction applied to bending angles to match derived bending angles.
!  See Eqn 4.4 of RSR 33.
!
! SEE ALSO
!   SAF/ROM/METO/REP/RSR/033 at http://www.romsaf.org/rsr.php.
!
! AUTHOR
!   Met Office, Exeter, UK and 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_fm_iono_bangle(Ne_max, R_peak, H_width, R_leo, n_L1, impact, bangle)

!-------------------------------------------------------------------------------
! 1. Declarations
!-------------------------------------------------------------------------------

  USE typesizes, ONLY: wp => EightByteReal
  USE ropp_fm,   not_this => ropp_fm_iono_bangle
  USE ropp_utils, ONLY: ropp_MDTV
  USE ropp_fm_constants, ONLY: k4, f_L1, f_L2

  IMPLICIT NONE

  REAL(wp), INTENT(in)                  :: Ne_max  ! Peak electron density (m-3)
  REAL(wp), INTENT(in)                  :: R_peak  ! Radius of peak value (m)
  REAL(wp), INTENT(in)                  :: H_width ! Width of Chapman Layer (m)
  REAL(wp), INTENT(in)                  :: R_leo   ! Radius of LEO satellite (m)
  INTEGER                               :: n_L1    ! Number of L1 bending angles
  REAL(wp), DIMENSION(:), INTENT(in)    :: impact  ! Impact parameters
  REAL(wp), DIMENSION(:), INTENT(inout) :: bangle  ! Updated bending angles

  INTEGER                               :: n_chap, n_impact, i
  REAL(wp)                              :: const_L1,const_L2
  REAL(wp)                              :: Ne_at_leo,uval
  REAL(wp), DIMENSION(:), ALLOCATABLE   :: lg, z, q, bangle_ion, bangle_leo

!-------------------------------------------------------------------------------
! 2. Useful variables
!-------------------------------------------------------------------------------

  n_chap = 1  ! For the moment
  n_impact = SIZE(impact)

  ALLOCATE ( lg(n_impact), z(n_impact), q(n_impact), &
             bangle_ion(n_impact), bangle_leo(n_impact) )

  bangle_ion = 0.0_wp
  bangle_leo = 0.0_wp

!-------------------------------------------------------------------------------
! 3. Calculate the zorro function
!-------------------------------------------------------------------------------

  lg = (R_peak - impact) / H_width

  z = ropp_fm_zorro(lg)

!--------------------------------------------------------------------------------
! 4. Ionospheric bending BUT NOT YET including scaling by frequency
!--------------------------------------------------------------------------------

  q = 2.0_wp * EXP(0.5_wp) * impact * R_peak / &
      SQRT(H_width*(R_peak + impact)**3)

  bangle_ion = bangle_ion + Ne_max * q * z

!--------------------------------------------------------------------------------
! 5. Estimate the electron density at the LEO using the Chapman layer
!--------------------------------------------------------------------------------

  uval = (R_leo - R_peak) / H_width    ! normalised vertical coordinate for Chapman

  Ne_at_leo = Ne_max * EXP( 0.5_wp * (1.0_wp - uval - EXP(-uval)) )

! estmate bending above LEO

  CALL ropp_fm_iono_bangle_above_leo (Ne_max, R_peak, H_width, R_leo, &
                                      impact, bangle_leo)

!--------------------------------------------------------------------------------
! 6. Now calculate total bending = (neutral + ionospheric)
!--------------------------------------------------------------------------------

  const_L1 = k4 / f_L1**2
  const_L2 = k4 / f_L2**2

  n_L1 = n_impact / 2

  DO i = 1, n_impact

    IF (bangle(i) < ropp_MDTV) CYCLE ! neutral is missing

    IF (i <= n_L1) THEN ! the L1 bending angles

      bangle(i) = bangle(i) + &
                  const_L1 * (bangle_ion(i) - bangle_leo(i)) - &  ! iono bending angle up to LEO
                  const_L1 * (Ne_at_leo*impact(i)/ SQRT(R_leo**2-impact(i)**2)) ! correction for electron density at LEO

    ELSE ! the L2 bending angles

      bangle(i) = bangle(i) + &
                  const_L2 * (bangle_ion(i) - bangle_leo(i)) - &  ! iono bending angle up to LEO
                  const_L2 * (Ne_at_leo*impact(i)/ SQRT(R_leo**2-impact(i)**2)) ! correction for electron density at LEO

    END IF

  END DO

!--------------------------------------------------------------------------------
! 7. Deallocate/tidy up
!--------------------------------------------------------------------------------

  DEALLOCATE ( lg, z, q, bangle_ion, bangle_leo )

  RETURN

END SUBROUTINE ropp_fm_iono_bangle