! $Id: ropp_pp_bending_angle_gnos.f90 2021 2009-01-16 10:49:04Z frhl $ SUBROUTINE ropp_pp_bending_angle_gnos(ro_data, config, gnos_ba, bangle, diag) !****s* WaveOptics/ropp_pp_bending_angle_gnos * ! ! NAME ! ropp_pp_bending_angle_gnos - Extrapolate GNOS L2 bending angles from ! fit to L2-L1. ! ! SYNOPSIS ! CALL ropp_pp_bending_angle_gnos(ro_data, config, gnos_ba, bangle, diag) ! ! DESCRIPTION ! This routine calculates L2 by extrapolating L1-L2 by using a simple model. ! ! INPUTS ! TYPE(ROprof) :: ro_data ! Input RO data structure ! TYPE(ppConfig) :: config ! Configuration options (not currently used) ! TYPE(L1btype) :: bangle ! WO-processed bending angles ! TYPE(L1btype) :: gnos_ba ! Non-WO-processed GNOS bending angles ! TYPE(PPDiag) :: diag ! Diagnostic input (includes L2_min_SLTA) ! ! OUTPUT ! TYPE(ROprof) :: ro_data ! RO data structure with modified QC value ! TYPE(L1btype) :: bangle ! overwritten L2 bending angles ! ! REFERENCES ! ROM SAF CDOP-2 VS32 Report ! SAF/ROM/DMI/REP/VS/32 ! ! 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 USE ropp_io_types, ONLY: ROprof, L1btype USE ropp_pp, not_this => ropp_pp_bending_angle_gnos USE ropp_pp_types, ONLY: ppConfig, ppDiag USE messages IMPLICIT NONE ! I/O TYPE(ROprof), INTENT(inout) :: ro_data ! Input RO data structure TYPE(ppConfig), INTENT(in) :: config ! Configuration options (not currently used) TYPE(L1btype), INTENT(in) :: gnos_ba ! Non-WO-processed GNOS bending angles TYPE(L1btype), INTENT(inout) :: bangle ! WO-processed bending angles TYPE(PPdiag), INTENT(in) :: diag ! Diagnostic input ! Local REAL(wp), ALLOCATABLE :: new_bangle_l1(:) ! Working bending angle REAL(wp), ALLOCATABLE :: H_factor(:) ! Working array, expensive to calculate INTEGER :: n ! Number of lev1a points INTEGER :: nl2 ! Number of points used in fitting REAL(wp) :: noise_estimate ! Goodness of least squares fit REAL(wp) :: Hy_sum, HH_sum, J_pen, x_opt ! Params from the LSQ fit REAL(wp) :: r_ion ! Distance to ionosphere peak REAL(wp), PARAMETER :: h_ion_default=300.0E3_wp ! Height of ionosphere peak REAL(wp), PARAMETER :: impact_height_min=20.0E3_wp ! Height above SLTA to start averaging REAL(wp), PARAMETER :: impact_height_max=70.0E3_wp ! Height to stop averaging REAL(wp), PARAMETER :: bangle_diff_max=1.0E-3_wp ! Min allowable bangle diff (rad) REAL(wp), PARAMETER :: noise_max=20.0E-6_wp ! Max allowable noise in fit (rad) LOGICAL, ALLOCATABLE :: lmask(:) ! for averaging operations LOGICAL :: dummy CHARACTER(LEN=5) :: str_nl2 CHARACTER(LEN=256) :: routine !------------------------------------------------------------------------------- ! 2. Overwrite L2 with the extrapolated value !------------------------------------------------------------------------------- CALL message_get_routine(routine) CALL message_set_routine('ropp_pp_bending_angle_gnos') CALL message(msg_info, 'Extrapolating GNOS L2 data from L2 - L1') dummy = config%obs_ok ! dummy use of otherwise unused argument n = ro_data%Lev1a%Npoints ALLOCATE (new_bangle_l1(n), lmask(n), H_factor(n)) ! 2.1 Interpolate pre-WO-processed L1 bangle onto L2 impact parameters CALL ropp_pp_interpol(gnos_ba%impact_l1, gnos_ba%impact_l2, & gnos_ba%bangle_l1, new_bangle_l1) ! 2.2 Define range over which fit will be calculated lmask = .FALSE. WHERE ( (gnos_ba%impact_l2-ro_data%georef%roc >= diag%L2_min_slta) .AND. & (gnos_ba%impact_l2-ro_data%georef%roc < & MIN(diag%L2_min_slta + impact_height_min, impact_height_max)) .AND. & (ABS(gnos_ba%bangle_l2 - new_bangle_l1) < bangle_diff_max) ) lmask = .TRUE. nl2 = COUNT(lmask) ! 2.3 Fit pre-WO L1-L2 and calculate the RMS to the observed data IF (nl2 > 0) THEN r_ion = ro_data%georef%roc + h_ion_default H_factor = r_ion / (r_ion**2 - gnos_ba%impact_l2**2)**1.5_wp Hy_sum = SUM((gnos_ba%bangle_l2 - new_bangle_l1) * H_factor, MASK=lmask) HH_sum = SUM(H_factor * H_factor, MASK=lmask) x_opt = Hy_sum / HH_sum ! Minimises | (gnos_ba%bangle_l2 - new_bangle_l1) - x*H_factor |**2 J_pen = SUM(( (gnos_ba%bangle_l2 - new_bangle_l1) - x_opt*H_factor )**2, MASK=lmask) noise_estimate = SQRT(MAX(J_pen, 1.0E-18_wp) / REAL(nl2, KIND=wp)) ELSE noise_estimate = ropp_MDFV END IF ! 2.4 Extrapolate L2 using the fit IF (nl2 > 10) THEN ! We have enough points ! Interpolate post-WO-processed L1 bangle onto L2 impact parameters CALL ropp_pp_interpol(bangle%impact_l1, bangle%impact_l2, & bangle%bangle_l1, new_bangle_l1) H_factor = r_ion / (r_ion**2 - bangle%impact_l2**2)**1.5_wp WHERE (bangle%impact_l2-ro_data%georef%roc < diag%L2_min_slta) & bangle%bangle_l2 = new_bangle_l1 + x_opt * H_factor ELSE WRITE (str_nl2, FMT='(I2)') nl2 CALL message(msg_warn, 'Not enough points (' // str_nl2 // ') to extrapolate L2') END IF !------------------------------------------------------------------------------- ! 3. Set overall quality value according to the noisiness of the fit !------------------------------------------------------------------------------- IF (noise_estimate > noise_max) ro_data%overall_qual = 20 !------------------------------------------------------------------------------- ! 4. Clean up !------------------------------------------------------------------------------- DEALLOCATE (H_factor, lmask, new_bangle_l1) CALL message_set_routine(routine) END SUBROUTINE ropp_pp_bending_angle_gnos