! $Id: ropp_pp_wopt_fsi_quick.f90 r4887 2016-05-25 15:48:28Z sti $ !****s* WaveOpticsPropagator/ropp_pp_wopt_fsi_quick * ! ! NAME ! ropp_pp_wopt_fsi_quick ! ! SYNOPSIS ! CALL ropp_pp_wopt_fsi_quick(kval, theta_dot, time, s_geom, & ! phase_leo, amp_leo, x_leo, y_leo, & ! x_gps, y_gps, & ! impact, alpha, amp) ! ! DESCRIPTION ! Calculates bending angles from excess phases at LEO using an FSI method. ! ! REFERENCES ! (1) Jensen, A. S., M. S. Lohmann, H.-H. Benzon, and A. S. Nielsen (2003), ! Full Spectrum Inversion of radio occultation signals, ! Radio Sci., 38, 1040, doi:10.1029/2002RS002763, 3. ! (2) Description of wave optics modelling in ROPP-9 and suggested ! improvements for ROPP-9.1, S. B. Healy, ROM SAF Report 28, ! SAF/ROM/ECMWF/REP/RSR/028 ! ! INPUTS ! REAL :: kval ! wavenumber (m^-1) ! REAL :: theta_dot ! (constant) angular velocity of LEO (rad/s) ! REAL :: time ! times over LEO orbit (s) ! REAL :: s_geom ! straight line distance between GNSS and LEO (m) ! REAL :: phase_leo ! accumulated phase at LEO (m) ! REAL :: amp_leo ! amplitude at LEO ! REAL :: x_leo, yleo ! coordinates of LEO ! REAL :: x_gps, ygps ! coordinates of GPS ! ! OUTPUT ! REAL :: impact ! impact heights (m) ! REAL :: alpha ! bending angle (rad) ! REAL :: amp ! amplitude of signal (rad) ! ! 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_fsi_quick(kval, theta_dot, time, s_geom, & phase_leo, amp_leo, x_leo, y_leo, & x_gps, y_gps, & impact, alpha, amp) USE typesizes, ONLY: wp => EightByteReal USE ropp_pp_constants, ONLY: pi1 => pi USE ropp_pp, ONLY: ropp_pp_fft IMPLICIT NONE ! Interface REAL(wp), INTENT(IN) :: kval ! wave number REAL(wp), INTENT(IN) :: theta_dot ! angular velocity REAL(wp), INTENT(IN) :: time(:) ! time at LEO REAL(wp), INTENT(IN) :: s_geom(:) ! straight line GNSS-to-LEO distance REAL(wp), INTENT(IN) :: x_leo(:),y_leo(:) ! position of LEO REAL(wp), INTENT(IN) :: x_gps(:),y_gps(:) ! position of GPS REAL(wp), INTENT(IN) :: phase_leo(:) ! accumulated phase at LEO REAL(wp), INTENT(IN) :: amp_leo(:) ! amplitude at LEO REAL(wp), INTENT(INOUT) :: impact(:) ! impact height (m) REAL(wp), INTENT(INOUT) :: alpha(:) ! bending angle (rad) REAL(wp), INTENT(INOUT) :: amp(:) ! amplitude (rad) ! Local COMPLEX(wp) :: U_final(SIZE(impact)) REAL(wp) :: phase_acc(SIZE(impact)) REAL(wp) :: arrival_time(SIZE(impact)) REAL(wp) :: twopi REAL(wp), PARAMETER :: amplitude_cutoff=1.0E-2_wp ! amplitudes < amplitude_cutoff*maximum_amp are ignored INTEGER :: ny ! number of points in profile INTEGER :: n_leo ! number of points in LEO orbit ! Really local INTEGER :: i,j,jj,n_interp,n_low,n_up REAL(wp) :: arg1,arg2,delta_t REAL(wp) :: phase_int,amp_int,wt,dopp_test,t_interp,mean_dopp REAL(wp) :: t_const,max_time,time_interval,min_time REAL(wp) :: rad_g,rad_l,theta REAL(wp) :: min_dopp,max_dopp REAL(wp) :: amp_max !------------------------------------------------------------------------------- ! 1.0 initialise !------------------------------------------------------------------------------- ny = SIZE(impact) n_leo = SIZE(x_leo) twopi = 2.0_wp * pi1 U_final(:) = (0.0_wp, 0.0_wp) ! upsample by a factor of up to 40 n_interp = MAX(2, MIN(40, (ny/(2*n_leo)))) t_interp = ABS(time(1)-time(2)) / n_interp !------------------------------------------------------------------------------- ! 2.0 estimate the Doppler offset !------------------------------------------------------------------------------- mean_dopp = 0.0_wp min_dopp = 1.0E18_wp max_dopp = -1.0E18_wp delta_t = ABS(time(1)-time(2)) amp_max = amplitude_cutoff * MAXVAL(amp_leo(:)) ! Bail out if amplitude smaller than this DO i = 2,n_leo IF ( amp_leo(i) < amp_max ) EXIT dopp_test = ( (phase_leo(i)-phase_leo(i-1)) + & (s_geom(i)-s_geom(i-1)) ) / delta_t dopp_test = dopp_test*kval/twopi min_dopp = MIN(min_dopp, dopp_test) max_dopp = MAX(max_dopp, dopp_test) mean_dopp = mean_dopp + dopp_test ENDDO ! just use the smallest Doppler mean_dopp = min_dopp !------------------------------------------------------------------------------- ! 3.0 map LEO signal to ny points by assuming linear variation ! of phase & amp between points !------------------------------------------------------------------------------- jj = 0 DO i = 1,n_leo - 1 jj = jj + 1 arg1 = kval * (phase_leo(i) + s_geom(i)) arg2 = twopi * REAL(jj-1, KIND=wp) * mean_dopp * t_interp U_final(jj) = amp_leo(i) * CMPLX(COS(arg1-arg2), SIN(arg1-arg2), KIND=wp) DO j = 1,n_interp-1 jj = jj + 1 ! linear interpolation between i and i+1 wt = REAL(n_interp-j, KIND=wp) / REAL(n_interp, KIND=wp) phase_int = wt*(phase_leo(i) + s_geom(i) ) + & (1.0_wp-wt)*(phase_leo(i+1) + s_geom(i+1)) amp_int = wt*amp_leo(i) + & (1.0_wp-wt)*amp_leo(i+1) arg1 = kval * phase_int arg2 = twopi * REAL(jj-1, KIND=wp) * mean_dopp * t_interp U_final(jj) = amp_int * CMPLX(COS(arg1-arg2), SIN(arg1-arg2), KIND=wp) ENDDO ENDDO !------------------------------------------------------------------------------- ! 4.0 calculate the FFT !------------------------------------------------------------------------------- CALL ropp_pp_fft(U_final, -1) !------------------------------------------------------------------------------- ! 5.0 identify the interval !------------------------------------------------------------------------------- n_low = 1 time_interval = REAL(ny, KIND=wp) * t_interp n_up = MIN(INT((max_dopp-min_dopp)*time_interval), ny) !------------------------------------------------------------------------------- ! 6.0 accumulate the phase values and differentiate !------------------------------------------------------------------------------- phase_acc(:) = 0.0_wp ! instantaneous phase DO i = n_low,n_up IF ( ABS(U_final(i)) > TINY(0.0_wp) ) phase_acc(i) = ATAN2(AIMAG(U_final(i)), REAL(U_final(i))) ENDDO ! accumulated phase DO i = n_up-1,n_low,-1 phase_acc(i) = phase_acc(i+1) + MODULO(phase_acc(i)-phase_acc(i+1)+pi1, twopi) - pi1 ENDDO ! times time_interval = REAL(ny, KIND=wp) * t_interp t_const = -0.5_wp * time_interval / twopi arrival_time = 0.0_wp ! limits on arrival time min_time = MINVAL(time) max_time = MAXVAL(time) ! differentiate phase DO i = n_low+1,n_up-1 arrival_time(i) = t_const * (phase_acc(i+1) - phase_acc(i-1)) ! apply limits arrival_time(i) = MIN(MAX(arrival_time(i), min_time), max_time) ENDDO ! reset limits n_low = 1 n_up = MIN(INT((max_dopp-min_dopp)*time_interval), ny) !------------------------------------------------------------------------------- ! 7.0 compute the impact heights, bending angles and amplitudes !------------------------------------------------------------------------------- impact(:) = 0.0_wp alpha(:) = 0.0_wp ! circular orbits, so radii are constant rad_l = SQRT(x_leo(1)**2 + y_leo(1)**2) rad_g = SQRT(x_gps(1)**2 + y_gps(1)**2) DO i = n_low, n_up impact(i) = (twopi / kval) * & (mean_dopp + (REAL(i-1, KIND=wp)/time_interval)) / theta_dot j = MAX(1, MIN(INT(arrival_time(i)/delta_t) + 1, n_leo-1)) ! compute the angle between the satellites, theta theta = ACOS( (x_leo(j)*x_gps(j) + y_leo(j)*y_gps(j)) / (rad_l*rad_g) ) theta = theta + theta_dot * (arrival_time(i)-time(j)) ! now compute bending angles alpha(i) = theta + ASIN(impact(i)/rad_l) + ASIN(impact(i)/rad_g) - pi1 amp(i) = ABS(U_final(i)) ENDDO END SUBROUTINE ropp_pp_wopt_fsi_quick