! $Id: ropp_pp_invert_lin.f90 3551 2013-02-25 09:51:28Z idculv $ SUBROUTINE ropp_pp_invert_LIN(impact, bangle, nr, refrac, scale) !****s* BendingAngle/ropp_pp_invert_LIN * ! ! NAME ! ropp_pp_invert_LIN - Calculate a one dimensional refractivity profile from ! bending angle / impact parameter profile ! using a Fast Abel Transform ! Assume linear variation of bending angle with height ! between successive impact parameter levels ! ! SYNOPSIS ! call ropp_pp_invert_LIN(impact, bangle, nr, refrac, scale) ! ! DESCRIPTION ! This routine calculates refractivity at a given set of x=nr levels ! from a vertical profile of bending angles given at a set of ! impact parameters. ! ! INPUTS ! real(wp), dimension(:) :: impact ! Input impact parameters ! real(wp), dimension(:) :: bangle ! Bending angles ! real(wp), dimension(:) :: nr ! x=nr product ! real(wp), optional :: scale ! Vertical scale height ! ! OUTPUT ! real(wp), dimension(:) :: refrac ! Refractivity values ! ! NOTES ! The interpolation of refractivity calculated at the input data ! impact height levels to the output geopotential levels is ! carried out assuming that dln(n)/dx varies linearly with x. ! ! SEE ALSO ! ropp_pp_invert_LIN ! ! AUTHOR ! M Gorbunov, Russian Academy of Sciences, Russia. ! Any comments on this software should be given via the ROM SAF ! Helpdesk at http://www.romsaf.org ! ! COPYRIGHT ! Copyright (c) 1998-2010 Michael Gorbunov ! For further details please refer to the file COPYRIGHT ! which you should have received as part of this distribution. ! !**** !------------------------------------------------------------------------------- ! 1. Declarations !------------------------------------------------------------------------------- USE typesizes, ONLY: wp => EightByteReal USE ropp_utils, ONLY: ropp_MDTV, ropp_MDFV, ropp_ZERO USE ropp_pp_constants, ONLY: pi IMPLICIT NONE REAL(wp), DIMENSION(:), INTENT(in) :: impact ! impact parameter REAL(wp), DIMENSION(:), INTENT(in) :: bangle ! bending angle REAL(wp), DIMENSION(:), INTENT(in) :: nr ! x=nr product for output REAL(wp), OPTIONAL, INTENT(in) :: scale ! scale height REAL(wp), DIMENSION(:), INTENT(out) :: refrac ! refractivity REAL(wp), DIMENSION(:), ALLOCATABLE :: kval ! linear decay rate REAL(wp) :: t_upper ! upper bound on integral REAL(wp) :: t_lower ! lower bound on integral REAL(wp) :: bangle_low ! bangle at lower level REAL(wp) :: impact_low ! impact at lower level REAL(wp) :: delta_bangle REAL(wp) :: delta_impact REAL(wp) :: zt REAL(wp) :: erf_up REAL(wp), PARAMETER :: a = 0.3480242_wp REAL(wp), PARAMETER :: b = 0.0958798_wp REAL(wp), PARAMETER :: c = 0.7478556_wp INTEGER :: n_lev, n_lower, n_impact INTEGER :: i, i_bot, i_top, l, k !------------------------------------------------------------------------------- ! 2. Useful variables !------------------------------------------------------------------------------- n_lev = SIZE(nr) n_impact = SIZE(impact) ALLOCATE(kval(n_impact)) !------------------------------------------------------------------------------- ! 3. Calculate lowest usable level (because of superrefraction) !------------------------------------------------------------------------------- n_lower = 1 DO i = n_impact, 2, -1 IF (impact(i) <= impact(i-1)) THEN n_lower = i EXIT ENDIF ENDDO !------------------------------------------------------------------------------- ! 4. Calculate linear decay rate between levels !------------------------------------------------------------------------------- DO i = 1, n_impact - 1 IF (bangle(i) < ropp_ZERO .OR. bangle(i+1) < ropp_ZERO) THEN kval(i) = 1.0e-12_wp ELSE ! Should probably use linear decay rate in this routine ! kval(i) = LOG(bangle(i)/bangle(i+1))/MAX(1.0_wp,(impact(i+1)-impact(i))) kval(i) = (bangle(i) - bangle(i+1))/MAX(1.0_wp,(impact(i+1)-impact(i))) kval(i) = MAX(1.0e-12_wp, kval(i)) ENDIF ENDDO !------------------------------------------------------------------------------- ! 5. Calculate bending angles for observational heights !------------------------------------------------------------------------------- refrac(:) = ropp_MDFV DO l = 1, n_lev IF (nr(l) < impact(n_lower) .OR. nr(l) > impact(n_impact)) THEN CYCLE ENDIF ! 5.1 Find bottom state vector level ! ---------------------------------- i_bot = n_lower DO k = i_bot, n_impact - 1 IF (nr(l) < impact(k + 1)) EXIT ENDDO i_bot = k i_top = n_impact DO WHILE (bangle(i_top) <= ropp_MDTV) i_top = i_top - 1 ENDDO ! 5.2 Loop over all levels above ! ------------------------------ refrac(l) = ropp_ZERO DO i = i_bot, i_top - 1 ! 5.2.1 Values of refractivity and impact parameter at lower level ! ---------------------------------------------------------------- IF (bangle(i) <= ropp_MDTV .OR. bangle(i+1) <= ropp_MDTV) CYCLE IF (i == i_bot) THEN ! Should probably use linear decay rate in this routine ! bangle_low = bangle(i_bot)*EXP(-kval(i_bot)*(nr(l) - impact(i_bot))) bangle_low = bangle(i_bot) - kval(i_bot)*(nr(l) - impact(i_bot)) impact_low = nr(l) ELSE bangle_low = bangle(i) impact_low = impact(i) ENDIF ! 5.2.2 Analytical Abel integral solution ! --------------------------------------- delta_bangle = bangle(i+1) - bangle_low delta_impact = impact(i+1) - impact_low t_upper = SQRT((impact(i+1)-nr(l))*(impact(i+1)+nr(l))) t_lower = SQRT((impact_low-nr(l))*(impact_low+nr(l))) refrac(l) = refrac(l) + & ((bangle_low*impact(i+1)-bangle(i+1)*impact_low) & *LOG((t_upper+impact(i+1))/(t_lower+impact_low)) + & delta_bangle*(t_upper-t_lower))/delta_impact ENDDO refrac(l) = refrac(l)/pi ! 5.3 Asymptotic correction ! ------------------------- IF(PRESENT(scale))THEN t_upper = SQRT((impact(i_top)-nr(l))/scale) zt = 1.0_wp / (1.0_wp + 0.47047_wp * t_upper) erf_up = 1.0_wp - (a-(b-c*zt)*zt) * zt * EXP(-(t_upper*t_upper)) refrac(l) = refrac(l) + & bangle(i_top)*EXP((impact(i_top)-nr(l))/scale)* & (1.0_wp-erf_up)/ SQRT(pi*(impact(i_top)+nr(l))/scale) ENDIF IF(refrac(l) == ropp_ZERO) refrac(l) = ropp_MDFV IF(refrac(l) > ropp_ZERO) refrac(l) = (1E6_wp)*(EXP(refrac(l))-1.0_wp) ENDDO DEALLOCATE(kval) END SUBROUTINE ropp_pp_invert_LIN