! $Id: ropp_pp_cutoff_amplitude.f90 2021 2009-01-16 10:49:04Z frhl $

SUBROUTINE ropp_pp_cutoff_amplitude(ro_data, LCF, config, amponly)

!****s* Preprocessing/ropp_pp_cutoff_amplitude *
!
! NAME
!    ropp_pp_cutoff_amplitude - Cut off occultation data based on amplitude
!                               and missing data flag
!
! SYNOPSIS
!    call ropp_pp_cutoff_amplitude(ro_data, LCF, config, amponly)
!
! DESCRIPTION
!    Cut off from amplitude (based on config%Acut parameter)
!
! INPUTS
!    type(ROprof)                      :: ro_data ! RO data strucuture
!    integer,  dimension(:)            :: LCF     ! Lost carrier flag
!    type(PPConfig)                    :: config  ! Configuration options
!    logical, optional                 :: amponly ! Flag to cut amplitude only

!
! OUTPUT
!    type(ROprof)                      :: ro_data ! Shrunk RO data strucuture
!    integer,  dimension(:)            :: LCF     ! Shrunk LCF
!    type(PPConfig)                    :: config  ! Configuration options
!
! 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
!
! 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
!
!   Updated;
!   Yong Chen, NOAA/STAR,
!   For full spectral inversion, truncate the signal dynamically using threshold method
!   1. Remove data above straight line impact height of 86 km
!   2. Remove data below L1 amplitude threshold defined below
!   3. Signal cut-off height is dependent on satellite alltitude (01/08/2024)
!
!     Satellite Altitude (km)	<400	400-450	450-500	500-600	>600
!    RFSI cut-off height (km)	-60	-80	-85	-90	-100
!
!
! COPYRIGHT
!   Copyright (c) 1998-2010 Michael Gorbunov <michael.gorbunov@zmaw.de>
!   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: impact_parameter
  USE messages
  USE ropp_io_types, ONLY: ROprof
  USE ropp_io, ONLY: ropp_io_shrink
  USE ropp_pp_preproc, not_this => ropp_pp_cutoff_amplitude
  USE ropp_pp
  USE ropp_pp_types, ONLY: PPconfig
  USE star_fsi, ONLY: runav

  IMPLICIT NONE

  TYPE(ROprof),       INTENT(inout) :: ro_data  ! RO data strucutre
  INTEGER,  DIMENSION(:), POINTER   :: LCF      ! Lost carrier flag
  TYPE(PPconfig),     INTENT(inout) :: config   ! Configuration options
  LOGICAL,  OPTIONAL, INTENT(in)    :: amponly  ! Flag to cut amplitude only
 
  LOGICAL                           :: amp_cut   ! amplitude cut only
  REAL(wp), PARAMETER :: tmin        = 25.0_wp   ! time limit 25 seconds
  REAL(wp), PARAMETER :: Pmax_limit  = 86.0_wp   ! impact height limit in km
!  REAL(wp), PARAMETER :: Pmax_limit  = 150.0_wp   ! impact height limit in km
  REAL(wp), PARAMETER :: ih1         = 60.0_wp   ! impact height limit for reference snr in km
  REAL(wp), PARAMETER :: ih2         = 80.0_wp   ! impact height limit for reference snr in km
!  REAL(wp), PARAMETER :: Minih       = -90.0_wp  ! lowest impact height limit in km for signal cut-off
  INTEGER,  PARAMETER :: nrh         = 500       ! number of impact height below max

  LOGICAL, PARAMETER  :: skipL2      =.FALSE.     ! default to skip L2 cut-off

  REAL(wp)                            :: p1, pN
  INTEGER                             :: i, j, n, iflag1, iflag2, ik
  INTEGER                             :: imax, imin, ibgn, ibgn0, ibgn1
  INTEGER                             :: w_smooth
  INTEGER                             :: ocd ! Occultation direction (1=rising)
  CHARACTER(len=8)                    :: imin_str, imax_str
  CHARACTER(len = 256)                :: routine
  REAL(wp), DIMENSION(:), ALLOCATABLE :: PH    ! Impact parameters [ch, point]
  REAL(wp), DIMENSION(:), ALLOCATABLE :: snrL1 ! Amplitude for L1 [ch, point]
  REAL(wp), DIMENSION(:), ALLOCATABLE :: snrL2 ! Amplitude for L2 [ch, point]
  REAL(wp), DIMENSION(:), ALLOCATABLE :: phL1  ! phase for L1 [ch, point]
  REAL(wp), DIMENSION(:), ALLOCATABLE :: phL2  ! phase for L2 [ch, point]
  REAL(wp), DIMENSION(:), ALLOCATABLE :: imPH  ! Impact parameters [ch, point]
  REAL(wp), DIMENSION(:), ALLOCATABLE :: time  ! time

  REAL(wp), DIMENSION(:), ALLOCATABLE :: ampL1sm   ! smoothed snr
  REAL(wp), DIMENSION(:), ALLOCATABLE :: ampL1sm1  ! smoothed snr
  REAL(wp), DIMENSION(:), ALLOCATABLE :: ampL2sm   ! smoothed snr
  REAL(wp), DIMENSION(:), ALLOCATABLE :: dphLL     ! phase difference between L1 and L2
  REAL(wp), DIMENSION(:), ALLOCATABLE :: dphLLsm   ! smoothed phase difference

  REAL(wp)                            :: Pmin, Pmax, dt, SNRref, avg0, avg1, val0, val1
  REAL(wp)                            :: max_ak0,  avg0_scale
  REAL(wp)                            :: std, tref, ak0, val2, maxL2snr
  INTEGER                             :: nmax1, i1, i2, ii, i20
  INTEGER                             :: nx

  REAL(wp)                            :: Minih     ! lowest impact height limit in km for signal cut-off
  REAL(wp)                            :: MeanSatAlt! Mean satellite altitude
  REAL(wp), DIMENSION(:), ALLOCATABLE :: SatAlt    ! Satellite Altitude [point]
  REAL(wp)                            :: SNRscale  ! SNR scaling factor


  CHARACTER(len=80) :: outstr

  CALL message_get_routine(routine)
  CALL message_set_routine('ropp_pp_cutoff_amplitude')

!-------------------------------------------------------------------------------
! 2. Initialisation
!-------------------------------------------------------------------------------

  n = ro_data%Lev1a%Npoints

  amp_cut = .FALSE.
  IF ( PRESENT(amponly) ) amp_cut = amponly

  p1 = impact_parameter( ro_data%Lev1a%r_leo(1,:)-ro_data%georef%r_coc, &
                         ro_data%Lev1a%r_gns(1,:)-ro_data%georef%r_coc )
  pN = impact_parameter( ro_data%Lev1a%r_leo(n,:)-ro_data%georef%r_coc, &
                         ro_data%Lev1a%r_gns(n,:)-ro_data%georef%r_coc )

  ocd = NINT(SIGN(1.0_wp, pN-p1))

  config%Pmax = MAX(p1, pN)
  config%Pmin = MIN(p1, pN)

  w_smooth = CEILING( config%fw_go_smooth*(n-1)/ABS(config%Pmax-config%Pmin) )

!-------------------------------------------------------------------------------
! FSI method
!-------------------------------------------------------------------------------

  IF ( INDEX(config%occ_method, "FS" ) == 1 .AND. .NOT. amp_cut ) THEN

    SNRscale = 1.0_wp
    ! Yong Chen on 05/10/2024 SNR need to divided by 10 only for C2 and PlanetIQ, others already take cared
    ! Yong Chen on 09/25/2024 adding C1
    ! Yong Chen on 12/19/2025 adding YM08 and STAT processing
    IF (ro_data%processing_centre(1:4) == 'UCAR' .OR. &
        ro_data%processing_centre(1:4) == 'NESD' .OR. &
        ro_data%processing_centre(1:4) == 'NOAA') THEN

       SELECT CASE (ro_data%leo_id(1:2))
         CASE ('C0', 'CO', 'C2', 'GN', 'YM')  ! C1, C2 and PlanetIQ 
            SNRscale = 0.1_wp
         CASE DEFAULT
           SNRscale = 1.0_wp
       END SELECT
    ENDIF

    iflag1 = 0
    iflag2 = 0
    ALLOCATE(PH(n))

    ALLOCATE(SatAlt(n))

    DO i = 1, n

      PH(i) = impact_parameter( ro_data%Lev1a%r_leo(i,:)-ro_data%georef%r_coc, &
                                ro_data%Lev1a%r_gns(i,:)-ro_data%georef%r_coc ) &
                                - ro_data%georef%roc

      ! Satellite altitude
      SatAlt(i) = Sqrt(Dot_Product(ro_data%Lev1a%r_leo(i,:)-ro_data%georef%r_coc, &
                                   ro_data%Lev1a%r_leo(i,:)-ro_data%georef%r_coc)) - ro_data%georef%roc


    ENDDO

    ! Mean satellite altitude (km)
    MeanSatAlt = SUM(SatAlt) / REAL(n, KIND=wp) * 0.001_wp

    ! define the signal cut-off height based on satellite altitude
    ! Yong Chen on 08/29/2024
    ! For Spire data, some satellite altitude can be lower than 250 km
    
    !IF ( MeanSatAlt <= 250.0_wp ) Minih = -60.0_wp
    !IF ( MeanSatAlt > 250.0_wp  .and. MeanSatAlt <= 300.0_wp ) Minih = -65.0_wp
    !IF ( MeanSatAlt > 300.0_wp  .and. MeanSatAlt <= 350.0_wp ) Minih = -70.0_wp
    !IF ( MeanSatAlt > 350.0_wp  .and. MeanSatAlt <= 400.0_wp ) Minih = -75.0_wp
    ! Yong Chen on 09/30/2024
    ! more aggresive to lower the cut-off height for low satellite altitude
    ! Yong Chen on 10/14/2024 for alt <400 and <350
    ! Yong Chen on 02/12/2025
    ! Based on Xinjia Zhou's investigation, the altitude between 350 and 450 km need to
    ! further fine tune
    IF ( MeanSatAlt <= 250.0_wp ) Minih = -25.0_wp
    IF ( MeanSatAlt > 250.0_wp  .and. MeanSatAlt <= 300.0_wp ) Minih = -35.0_wp
    IF ( MeanSatAlt > 300.0_wp  .and. MeanSatAlt <= 350.0_wp ) Minih = -45.0_wp
    IF ( MeanSatAlt > 350.0_wp  .and. MeanSatAlt <= 380.0_wp ) Minih = -60.0_wp
    IF ( MeanSatAlt > 380.0_wp  .and. MeanSatAlt <= 400.0_wp ) Minih = -65.0_wp
    IF ( MeanSatAlt > 400.0_wp  .and. MeanSatAlt <= 420.0_wp ) Minih = -70.0_wp
    IF ( MeanSatAlt > 420.0_wp  .and. MeanSatAlt <= 450.0_wp ) Minih = -80.0_wp
    IF ( MeanSatAlt > 450.0_wp  .and. MeanSatAlt <= 500.0_wp ) Minih = -85.0_wp
    IF ( MeanSatAlt > 500.0_wp  .and. MeanSatAlt <= 600.0_wp ) Minih = -90.0_wp
    IF ( MeanSatAlt > 600.0_wp ) Minih = -100.0_wp

    config%Pmax = MAX(PH(1), PH(n))
    config%Pmin = MIN(PH(1), PH(n))

   ! convert to km
    Pmax = config%Pmax * 0.001_wp
    Pmin = config%Pmin * 0.001_wp

    WRITE( outstr, '(F14.1,1X,F14.1,1X,F14.1,1X,F14.1)') Pmin, Pmax, MeanSatAlt, Minih

    CALL message(msg_info, 'Pmin, Pmax, SatH, cut-off height: ' // TRIM(outstr))

    !-------------------------------------------------------------------------------
    ! 3. Cut-off data from amplitude (snr_L1p) and missing data flag (LCF)
    !-------------------------------------------------------------------------------

    ! 3.1 Determine cut-off limits

    ALLOCATE(snrL1(n))
    ALLOCATE(snrL2(n))
    ALLOCATE(phL1(n))
    ALLOCATE(phL2(n))
    ALLOCATE(imPh(n))
    ALLOCATE(time(n))

    ALLOCATE(ampL1sm(n))
    ALLOCATE(ampL1sm1(n))
    ALLOCATE(ampL2sm(n))
    ALLOCATE(dphLL(n))
    ALLOCATE(dphLLsm(n))

    imin = 1
    imax = n

    IF ( ocd < 0 ) THEN     ! setting occultation

       ! change to rising occultation, reverse the array
       j = 1
       DO i=n,1,-1
          time(j)  = ro_data%Lev1a%dtime(i)
          snrL1(j) = ro_data%Lev1a%snr_L1ca(i)*SNRscale
          snrL2(j) = ro_data%Lev1a%snr_L2p(i)*SNRscale
          phL1(j)  = ro_data%Lev1a%phase_L1(i)
          phL2(j)  = ro_data%Lev1a%phase_L2(i)
          imPh(j)  = PH(i) * 0.001_wp
          j = j + 1

       ENDDO

    ELSE                    ! rising occultation

!       imax = n
       DO i=1,n
          time(i)  = ro_data%Lev1a%dtime(i)
          snrL1(i) = ro_data%Lev1a%snr_L1ca(i)*SNRscale
          snrL2(i) = ro_data%Lev1a%snr_L2p(i)*SNRscale
          phL1(i)  = ro_data%Lev1a%phase_L1(i)
          phL2(i)  = ro_data%Lev1a%phase_L2(i)
          imPh(i)  = PH(i) * 0.001_wp

       ENDDO

    ENDIF

    tref = abs(time(n) - time(1))
    dt = tref/(n - 1)

    IF (Pmax < 35.0_wp .OR. Pmin > 30.0_wp .OR. tref < tmin) THEN
       iflag1 = 1
    ENDIF

    IF (iflag1 == 0) THEN
      IF (Pmax > Pmax_limit) THEN
         nmax1 = MAX(SUM(MINLOC(impH(:), impH(:) <= Pmax_limit - 1.0_wp)),  &
                     SUM(MAXLOC(impH(:), impH(:) <= Pmax_limit - 1.0_wp)))
      ELSE
         nmax1 = n
      ENDIF

      IF (Pmax < ih1) THEN
         i1 = nmax1 - nrh
         i2 = nmax1
      ELSE
         i1 = MIN(SUM(MINLOC(impH(:), ih2 > impH(:)  .AND. impH(:) > ih1)),  &
                  SUM(MAXLOC(impH(:), ih2 > impH(:)  .AND. impH(:) > ih1)))
         i2 = MAX(SUM(MINLOC(impH(:), ih2 > impH(:)  .AND. impH(:) > ih1)),  &
                  SUM(MAXLOC(impH(:), ih2 > impH(:)  .AND. impH(:) > ih1)))
      ENDIF
      ibgn = i1
      ik = i2 - i1 + 1
!      WRITE ( * , * ) '60 km to 80 km: n, nmax1, i1, i2', n, nmax1, i1, i2

      ! calculate reference SNR
      SNRref = SUM( snrL1(i1:i2)) / ik
      std = SQRT( SUM( snrL1(i1:i2)**2) / ik - SNRref**2 )
      ! Apply a 3-seconds smoothing to L1 SNR
      nx = 300
      IF (abs(dt - 0.02) < 1.D-5) nx = 150
      CALL runav(nx,  snrL1(1:nmax1), ampL1sm(1:nmax1))
      CALL runav(nx,  snrL2(1:nmax1), ampL2sm(1:nmax1))

      nx = 100
      IF (abs(dt - 0.02) < 1.D-5) nx = 50
      CALL runav(nx,  snrL1(1:nmax1), ampL1sm1(1:nmax1))
      ! Calculate base (noise) SNR value
      i1 = 50
      i2 = 1000
      IF (abs(dt - 0.02) <1.D-5) THEN
         i1 = 25
         i2 = 500
      ENDIF
      i20 = i1
      IF ( ampL1sm1(i20) < 0.2_wp * SNRref) then
        DO WHILE (( ampL1sm1(i20) < 0.2_wp * SNRref) .AND. (i20 <= i2) )
          i20 = i20 + 1
        ENDDO
        i2 = i20
      ENDIF
      ik = i2 - i1 + 1

      avg0 = SUM( ampL1sm(1:ik)) / ik
      avg1 = SUM( ampL2sm(1:ik)) / ik
      if (avg1 >= 15.0_wp) avg1 = 15.0_wp

      WRITE( outstr, '(I5,1X,I5,1X,F14.1,1X,F14.1,1X,F14.1)') i1, i2, SNRref, avg0, avg1

      CALL message(msg_info, 'i1, i2, reference/base SNR: ' // TRIM(outstr))

      ! Use all the profile if average SNR > threshold
      IF (avg0 > 50.0_wp) THEN
        ibgn0 = 1
        avg0 = 12.5
      ELSE
      ! specification of thresholds
        val0 = 3.0_wp * avg0
        val1 = 1.5_wp * avg0
        IF ( SNRref >= 1500.0_wp) THEN
          val0 = 4.0_wp * avg0
          val1 = 2.0_wp * avg0
        ENDIF

        ! Remove low SNR profiles
        IF (SNRref <= val0) THEN
           iflag2 = 1
        ELSE
           ! First truncation of noise
           DO WHILE ( (ampL1sm1(i1) <= val0) .OR. (snrL1(i1) < avg0) )
             i1 = i1 + 1
           ENDDO
           ibgn0 = i1 - 1
           IF ( ( nmax1 - ibgn0 < 1000) .OR. (ibgn < ibgn0 ) ) THEN
              iflag2 = 2
           ELSE
              DO WHILE ( (ampL1sm1(ibgn0) >= val1) .AND. (snrL1(ibgn0) > 0.0_wp) )
                ibgn0 = ibgn0 - 1
              ENDDO
           ENDIF
           WRITE ( * , * ) 'here, ibgn0 =', ibgn0, ampL1sm1(ibgn0), impH(ibgn0)

        ENDIF

      ENDIF

      ! ychen on 05/07 max_ak0 and avg0_scale are mission dependent
      max_ak0 = 7.5_wp
      avg0_scale = 10.0_wp
      
      ! ychen on 11/4/2024 to refine UCAR metop atmPhs low level SNR cutoff threshold
      ! considering metop satellite atlitude is around 830 km
      IF (ro_data%processing_centre(1:4) == 'UCAR' &
          .AND. (ro_data%leo_id(1:2) == 'ME' .OR. ro_data%leo_id(1:2) == 'MT')) THEN
        max_ak0 = 5.0_wp
        avg0_scale = 6.0_wp
      ENDIF
 
      IF (iflag2 == 0) THEN
        ! Now remove anomalous cases
        ak0 = 2.0_wp
!        WRITE ( * , * ) 'ibgn0 =', ibgn0, impH(ibgn0)

        IF (impH(ibgn0) < Minih) then

           DO WHILE ( ak0 <= max_ak0)
              val0 = avg0_scale * avg0
              val1 = ak0 * avg0
              i1 = ibgn0
              DO WHILE (ampL1sm(i1) <= val0)
               i1 = i1 + 1
              ENDDO
              IF (i1 > 1) ibgn0 = i1 - 1
              IF (ibgn0 <1) ibgn0 = 1
              DO WHILE (ampL1sm(ibgn0) >= val1)
               ibgn0 = ibgn0 - 1
               IF (ibgn0 < 1) EXIT
              ENDDO

              IF (ibgn0 < 1) ibgn0 = 1
              WRITE ( * , * ) 'in while loop ak0, ibgn0 =', ak0, ibgn0, ampL1sm1(ibgn0), impH(ibgn0)

              IF (impH(ibgn0) <= Minih) then
                 ak0 = ak0 + 0.25_wp
!                 WRITE ( * , * ) 'in while loop ak0, ibgn0 =', ak0, ibgn0, impH(ibgn0)
              ELSE
                 EXIT
              ENDIF
           ENDDO

        ENDIF
        IF (ak0 > max_ak0) then
           DO WHILE (impH(ibgn0) <= Minih)
             ibgn0 = ibgn0 + 1
           ENDDO
        ENDIF

        WRITE( outstr, '(I5,1X,F14.1,1X,F14.1)') ibgn0, impH(ibgn0) , impH(nmax1)

        CALL message(msg_info, 'L1 starting point and new pmin/pmax: ' // TRIM(outstr))

        ! Yong Chen on 05/09/2024, don't apply L2 cut-off since L2 if already
        ! corrected before calling this function
        ! Now for L2 data
        IF (.NOT. skipL2) THEN
          nx = 100
          CALL runav(nx,  snrL2(1:nmax1), ampL2sm(1:nmax1))
          val2 = 3.0_wp*avg1
          ibgn1 = ibgn0
          maxL2snr = maxval(ampL2sm(1:nmax1))

          IF (maxL2snr <= val2) THEN
           ibgn1 = ibgn0
          ELSE
           DO WHILE(ampL2sm(ibgn1) < val2)
            ibgn1 = ibgn1+1
           ENDDO
          ENDIF

          WRITE( outstr, '(I5,1X,F14.1)') ibgn1,impH(ibgn1)

          CALL message(msg_info, 'L2 starting point and impact height: ' // TRIM(outstr))

          IF (.NOT. config%opt_DL2) THEN
            ! Now correct L2 phase
            dphLL(1:nmax1) = phL2(1:nmax1) - phL1(1:nmax1)
            nx = 50
            call runav(nx,dphLL(1:nmax1),dphLLsm(1:nmax1))
            phL2(1:nmax1) = phL1(1:nmax1) + dphLLsm(1:nmax1)
            ! Replace L2 values from bgn0 to bgn1 by 0
            if (ibgn1 >= ibgn0) snrL2(ibgn0:ibgn1) = 0.0_wp
          ENDIF
        ENDIF
        ! Determine cut-off limits
        IF ( ocd < 0 ) THEN     ! setting occultation
          imin = n + 1 - nmax1
          imax = n + 1 - ibgn0

          IF (.NOT. config%opt_DL2) THEN
            ro_data%Lev1a%snr_L2p(imin:imax) = 1.0_wp/SNRscale * snrL2(nmax1:ibgn0:ocd)
            ro_data%Lev1a%phase_L2(imin:imax) = phL2(nmax1:ibgn0:ocd)
          ENDIF

          DO i= n + 1 - nmax1, n + 1 - ibgn0
            IF (BTEST(LCF(i),3)) THEN
              imax = MAX(n + 1 - nmax1, MIN(i-1, imax))
              EXIT
            ENDIF
          ENDDO

        ELSE                    ! rising occultation
          imin = ibgn0
          imax = nmax1

          IF (.NOT. config%opt_DL2) THEN
            ro_data%Lev1a%snr_L2p(imin:imax) = 1.0_wp/SNRscale * snrL2(ibgn0:nmax1)
            ro_data%Lev1a%phase_L2(imin:imax) = phL2(ibgn0:nmax1)
          ENDIF

          DO i=nmax1,ibgn0,-1
            IF (BTEST(LCF(i),3)) THEN
              imin = MIN(nmax1, MAX(i+1, imin))
              EXIT
            ENDIF
          ENDDO
        ENDIF

      ENDIF
    ENDIF

    DEALLOCATE(SatAlt)

    DEALLOCATE(PH)
    DEALLOCATE(snrL1)
    DEALLOCATE(snrL2)
    DEALLOCATE(phL1)
    DEALLOCATE(phL2)
    DEALLOCATE(imPh)
    DEALLOCATE(time)
    DEALLOCATE(ampL1sm)
    DEALLOCATE(ampL1sm1)
    DEALLOCATE(ampL2sm)
    DEALLOCATE(dphLL)
    DEALLOCATE(dphLLsm)

  ELSE
  !-------------------------------------------------------------------------------
  ! 3. Cut-off data from amplitude (snr_L1p) and missing data flag (LCF)
  !-------------------------------------------------------------------------------
  ! for othre methods (GO and WO)
    imin = 1
    imax = n

    IF ( ocd < 0 ) THEN     ! setting occultation

       imin = 1
       DO i=n,1,-1
          IF (ro_data%Lev1a%snr_L1ca(i) >                               &
                   MAXVAL(ro_data%Lev1a%snr_L1ca(:))*config%Acut) THEN
             imax = MIN(i + 2*w_smooth, n)
             EXIT
          ENDIF
       ENDDO
       DO i=1,n
         IF (BTEST(LCF(i),3)) THEN
           imax = MAX(1, MIN(i-1, imax))
           EXIT
         ENDIF
       ENDDO

    ELSE                    ! rising occultation

       imax = n
       DO i=1,n
          IF (ro_data%Lev1a%snr_L1ca(i) >                               &
                   MAXVAL(ro_data%Lev1a%snr_L1ca(:))*config%Acut) THEN
             imin = MAX(i - 2*w_smooth, 1)
             EXIT
          ENDIF
       ENDDO
       DO i=n,1,-1
         IF (BTEST(LCF(i),3)) THEN
           imin = MIN(n, MAX(i+1, imin))
           EXIT
         ENDIF
       ENDDO

    ENDIF

  ENDIF

  ! 3.2 Select data subset

  IF ((imin > 1 .OR. imax < n) .AND. (imax >= imin)) THEN
     WRITE(imin_str,  '(i8)') imin
     WRITE(imax_str,  '(i8)') imax
     CALL message(msg_info,"Cut-off (amplitude/LCF criterion). Keep data " // &
                           imin_str  // " to " // imax_str)

     CALL ropp_io_shrink(ro_data%Lev1a, imin, imax, 1)
     CALL shrink_varint(lcf, imin, imax, 1)
     n = ro_data%Lev1a%Npoints

  ENDIF

  !IF ( INDEX(config%occ_method, "FS" ) == 1 ) THEN

  !  open(10,file='truncated_phase1.txt')
  !  write(10,'(F20.9)') ro_data%georef%roc/1000.0_wp
  !  write(10,'(2F20.5)') ro_data%georef%lat,  ro_data%georef%lon
  !  write(10,'(3F20.5)') ro_data%georef%r_coc/1000.0_wp
  !  write(10,*)'time  xleo  yleo  zleo  xgns  ygns  zgns  phL1  '// &
  !    'phL2  phLC  snrL1  snrL2  OLphs'
  !  do i=1, n
  !    write(10,'(F12.8,6F20.11,2F20.7,2F10.2)') &
  !      ro_data%Lev1a%dtime(i), &
  !      (ro_data%Lev1a%r_leo(i,1)-ro_data%georef%r_coc(1))/1000.0_wp, &
  !      (ro_data%Lev1a%r_leo(i,2)-ro_data%georef%r_coc(2))/1000.0_wp, &
  !      (ro_data%Lev1a%r_leo(i,3)-ro_data%georef%r_coc(3))/1000.0_wp, &
  !      (ro_data%Lev1a%r_gns(i,1)-ro_data%georef%r_coc(1))/1000.0_wp, &
  !      (ro_data%Lev1a%r_gns(i,2)-ro_data%georef%r_coc(2))/1000.0_wp, &
  !      (ro_data%Lev1a%r_gns(i,3)-ro_data%georef%r_coc(3))/1000.0_wp, &
  !      ro_data%Lev1a%phase_L1(i), ro_data%Lev1a%phase_L2(i), &
  !      ro_data%Lev1a%snr_L1ca(i)*0.1_wp, ro_data%Lev1a%snr_L2p(i)*0.1_wp
  !
  !  enddo
  !  close(10)
  !
  !ENDIF

  CALL message_set_routine(routine)

CONTAINS

!-------------------------------------------------------------------------------
! 6. Select data subset for additional variables
!-------------------------------------------------------------------------------

  SUBROUTINE shrink_var(var, imin, imax, stride)

    USE typesizes, ONLY: wp => EightByteReal
    USE ropp_utils, ONLY: copy_alloc

    IMPLICIT NONE

    REAL(wp), DIMENSION(:), POINTER       :: var
    INTEGER,                INTENT(in)    :: imin
    INTEGER,                INTENT(in)    :: imax
    INTEGER,                INTENT(in)    :: stride
    REAL(wp), DIMENSION(:), POINTER       :: tmp => null()

    CALL copy_alloc(var(imin:imax:stride), tmp)
    DEALLOCATE(var)

    CALL copy_alloc(tmp, var)
    DEALLOCATE(tmp)

  END SUBROUTINE shrink_var

  SUBROUTINE shrink_varint(var, imin, imax, stride)

    USE ropp_utils, ONLY: copy_alloc

    IMPLICIT NONE

    INTEGER, DIMENSION(:), POINTER      :: var
    INTEGER,               INTENT(in)   :: imin
    INTEGER,               INTENT(in)   :: imax
    INTEGER,               INTENT(in)   :: stride
    INTEGER, DIMENSION(:), POINTER      :: tmp => null()

    CALL copy_alloc(var(imin:imax:stride), tmp)
    DEALLOCATE(var)

    CALL copy_alloc(tmp, var)
    DEALLOCATE(tmp)

  END SUBROUTINE shrink_varint


END SUBROUTINE ropp_pp_cutoff_amplitude