! $Id: ropp_pp_preprocess_GNOS.f90 2021 2009-01-16 10:49:04Z frhl $ SUBROUTINE ropp_pp_preprocess_GNOS(ro_data, config, diag) !****s* Preprocessing/ropp_pp_preprocess_GNOS * ! ! NAME ! ropp_pp_preprocess_GNOS - Mission-specific Level1a data preprocessing ! for GNOS data ! ! SYNOPSIS ! CALL ropp_pp_preprocess_GNOS(ro_data, config, diag) ! ! DESCRIPTION ! 1) Remove profiles with sudden jumps in the LEO position. ! 2) Calculate lowest valid L2 SLTA. ! 3) Calculate some stats on phases and SNRs between 60 and 80 km. ! ! INPUTS ! type(ROprof) :: ro_data ! Radio occultation data strucuture ! type(PPConfig) :: config ! Configuration options ! type(PPDiag) :: diag ! Diagnostic input ! ! OUTPUT ! type(ROprof) :: ro_data ! Corrected radio occultation data ! type(PPConfig) :: config ! Configuration options ! type(PPDiag) :: diag ! Diagnostic output ! ! NOTES ! Requires ROprof data structure type, defined in ropp_io module. This ! routine therefore requires that the ropp_io module is pre-installed before ! compilation. ! ! 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_io_types, ONLY: ROprof USE ropp_pp_preproc, not_this => ropp_pp_preprocess_GNOS USE ropp_pp_types, ONLY: PPconfig, PPdiag USE ropp_utils IMPLICIT NONE ! I/O TYPE(ROprof), INTENT(inout) :: ro_data ! Radio occultation data structure TYPE(PPconfig), INTENT(inout) :: config ! Configuration options TYPE(PPdiag), INTENT(inout) :: diag ! Diagnostic I/O ! Local INTEGER :: n, i, k, nl1, nl2 REAL(wp) :: cos_angle, slta_param REAL(wp), PARAMETER :: cos_10deg=0.985_wp REAL(wp) :: ps1, psN ! impact parameters at beginning and end REAL(wp) :: ave_phsl1 ! L1 statistics REAL(wp) :: ave_phsl2 ! L2 statistics REAL(wp) :: Pleo(3), Pgns(3), Pver(3) ! impact heights REAL(wp) :: min_slta_G REAL(wp), ALLOCATABLE :: slta(:) ! SLTA along profile REAL(wp), PARAMETER :: slta_buffer=6.0E3_wp ! Added to computed SLTA to get min SLTA REAL(wp), PARAMETER :: slta_min=60.0E3_wp ! Minimum SLTA used in averaging REAL(wp), PARAMETER :: slta_max=80.0E3_wp ! Maximum SLTA used in averaging REAL(wp), PARAMETER :: slta_threshold=50.0E3_wp ! Poor quality if SLTA > this REAL(wp), PARAMETER :: phs1_threshold=150.0_wp ! Rising profs are poor quality if || < this REAL(wp), PARAMETER :: phs2_threshold=150.0_wp ! Rising profs are poor quality if || < this LOGICAL, ALLOCATABLE :: lmask(:) ! for averaging operations CHARACTER(LEN = 10) :: delta_angle CHARACTER(LEN = 10) :: str_slta CHARACTER(LEN = 256) :: routine !------------------------------------------------------------------------------- ! 2. Remove any profiles with 'jumps' in the LEO positions, ! diagnosed by a big angle subtended by the LEO over its orbit !------------------------------------------------------------------------------- CALL message_get_routine(routine) CALL message_set_routine('ropp_pp_preprocess_GNOS') n = ro_data%Lev1a%Npoints cos_angle = DOT_PRODUCT(ro_data%Lev1a%r_leo(1, :), ro_data%Lev1a%r_leo(n, :)) / & SQRT( DOT_PRODUCT(ro_data%Lev1a%r_leo(1, :), ro_data%Lev1a%r_leo(1, :)) * & DOT_PRODUCT(ro_data%Lev1a%r_leo(n, :), ro_data%Lev1a%r_leo(n, :)) ) IF (cos_angle < cos_10deg) THEN ! i.e. if the angular displacement is more than about 10 degrees WRITE (delta_angle, FMT='(ES10.3)') ACOS(cos_angle) * 45.0_wp / ATAN(1.0_wp) CALL message(msg_warn, 'Line of sight to LEO changes too much (' // delta_angle // ' deg) ' // & 'during its orbit. Exiting processing.') ro_data%Lev1a%Npoints = 0 config%obs_ok = .FALSE. CALL message_set_routine(routine) RETURN ENDIF !------------------------------------------------------------------------------- ! 3. Calculate (minimum) SLTA !------------------------------------------------------------------------------- ALLOCATE (slta(n), lmask(n)) min_slta_G = HUGE(1.0_wp) DO i=1,n DO k=1,3 Pleo(k) = ro_data%lev1a%r_leo(i, k) - ro_data%georef%r_coc(k) Pgns(k) = ro_data%lev1a%r_gns(i, k) - ro_data%georef%r_coc(k) END DO slta_param = DOT_PRODUCT( -Pleo, Pgns-Pleo) / & ! Where normal to PgnsPleo goes through O. DOT_PRODUCT(Pgns-Pleo, Pgns-Pleo) Pver = Pleo + slta_param * (Pgns - Pleo) slta(i) = SQRT(DOT_PRODUCT(Pver, Pver)) - ro_data%georef%roc IF ( (ro_data%lev1a%snr_L2P(i) >= ropp_MDTV) .AND. & ! amplitude present (ro_data%lev1a%phase_L2(i) >= ropp_MDTV) .AND. & ! phase present (ABS(ro_data%lev1a%phase_L2(i) + 999.0_wp) >= 0.001_wp) ) THEN ! CMA phase present min_slta_G = MIN(min_slta_G, slta(i)) END IF END DO IF ( min_slta_G < (HUGE(1.0_wp) - slta_buffer) ) THEN diag%L2_min_slta = MAX( min_slta_G + slta_buffer, & ! Add 6 km buffer config%hmax_wo + 0.5_wp*config%fw_go_full ) WRITE (str_slta, FMT='(F8.3)') diag%L2_min_slta * 1.E-3_wp CALL message(msg_info, 'Minimum valid L2 SLTA = ' // str_slta // ' km.') ELSE diag%L2_min_slta = ropp_MDFV END IF !------------------------------------------------------------------------------- ! 4. Calculate mean L1 and L2 excess phases between slta_min and slta_max !------------------------------------------------------------------------------- ! 4.1 L1 variables lmask = .FALSE. WHERE ( (slta > slta_min) .AND. & (slta < slta_max) .AND. & (ro_data%lev1a%phase_L1 > ropp_MDTV) .AND. & (ro_data%lev1a%snr_L1ca > ropp_ZERO) ) lmask = .TRUE. nl1 = COUNT(lmask) IF (nl1 > 1) THEN ave_phsl1 = SUM(ro_data%lev1a%phase_L1, MASK=lmask) / REAL(nl1, KIND=wp) ELSE ave_phsl1 = ropp_MDFV END IF ! 4.2 L2 variables lmask = .FALSE. WHERE ( (slta > slta_min) .AND. & (slta < slta_max) .AND. & (ro_data%lev1a%phase_L2 > ropp_MDTV) .AND. & (ro_data%lev1a%snr_L2p > ropp_ZERO) ) lmask = .TRUE. nl2 = COUNT(lmask) IF (nl2 > 1) THEN ave_phsl2 = SUM(ro_data%lev1a%phase_L2, MASK=lmask) / REAL(nl2, KIND=wp) ELSE ave_phsl2 = ropp_MDFV END IF !------------------------------------------------------------------------------- ! 5. Set some preliminary QC on the basis of these results !------------------------------------------------------------------------------- ps1 = impact_parameter(ro_data%lev1a%r_leo(1,:) - ro_data%georef%r_coc(:), & ro_data%lev1a%r_gns(1,:) - ro_data%georef%r_coc(:)) psN = impact_parameter(ro_data%lev1a%r_leo(n,:) - ro_data%georef%r_coc(:), & ro_data%lev1a%r_gns(n,:) - ro_data%georef%r_coc(:)) IF ( ( (psN > ps1) .AND. & (ABS(ave_phsl1) < phs1_threshold) .AND. & (ABS(ave_phsl2) < phs2_threshold) ) .OR. & (min_slta_G > slta_threshold) ) THEN ro_data%overall_qual = 20 ELSE ro_data%overall_qual = 100 END IF !------------------------------------------------------------------------------- ! 6. Clean up !------------------------------------------------------------------------------- DEALLOCATE (lmask, slta) CALL message_set_routine(routine) END SUBROUTINE ropp_pp_preprocess_GNOS