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

SUBROUTINE ropp_pp_cutoff_amplitude(ro_data, LCF, config)

!****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)
!
! 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

!
! 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
!
! 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

  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 :: 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
  INTEGER,  PARAMETER :: nrh         = 500       ! number of impact height below max
  
  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)                            :: std, tref, ak0, val2, maxL2snr  
  INTEGER                             :: nmax1, i1, i2, ii, i20 
  INTEGER                             :: nx  

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

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

  n = ro_data%Lev1a%Npoints
 
  iflag1 = 0 
  iflag2 = 0 
  ALLOCATE(PH(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
  
  ENDDO
 
  ocd = NINT(SIGN(1.0_wp, PH(n)-PH(1)))

  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 ( * , * ) 'pmin/pmax =', Pmin, Pmax, ro_data%georef%roc, ro_data%georef%r_coc
  
  w_smooth = CEILING( config%fw_go_smooth*(n-1)/ABS(config%Pmax-config%Pmin) )

  IF ( INDEX(config%occ_method, "FS" ) == 1 ) THEN
!-------------------------------------------------------------------------------
! 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)*0.1_wp
          snrL2(j) = ro_data%Lev1a%snr_L2p(i)*0.1_wp
          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
!          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
          time(i)  = ro_data%Lev1a%dtime(i) 
          snrL1(i) = ro_data%Lev1a%snr_L1ca(i)*0.1_wp
          snrL2(i) = ro_data%Lev1a%snr_L2p(i)*0.1_wp
          phL1(i)  = ro_data%Lev1a%phase_L1(i)
          phL2(i)  = ro_data%Lev1a%phase_L2(i)
          imPh(i)  = PH(i) * 0.001_wp
        
       
!          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

    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 ( * , * ) 'i1, i2', i1, i2 
      WRITE ( * , * ) 'reference/base SNR =', SNRref, avg0, avg1 
      
      ! 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, impH(ibgn0) 
              
        ENDIF 

      ENDIF 
      
      IF (iflag2 == 0) THEN
        ! Now remove anomalous cases                                      
        ak0 = 2.0_wp 
!        WRITE ( * , * ) 'ibgn0 =', ibgn0, impH(ibgn0) 
        
        IF (impH(ibgn0) < -90.0_wp) then 
        
           DO WHILE ( ak0 <= 7.5_wp)
              val0 = 10.0_wp * 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 
              
              IF (impH(ibgn0) <= -90.0_wp) then 
                 ak0 = ak0 + 0.25_wp 
!                 WRITE ( * , * ) 'in while loop ak0, ibgn0 =', ak0, ibgn0, impH(ibgn0) 
              ELSE
                 EXIT
              ENDIF 
           ENDDO 

        ENDIF
        IF (ak0 > 7.5_wp) then                    
           DO WHILE (impH(ibgn0) <= -90.0_wp)    
             ibgn0 = ibgn0 + 1                     
           ENDDO                               
        ENDIF                                  
        WRITE ( * , * ) 'ibgn0 =', ibgn0 
        WRITE ( * , * ) 'new pmin/pmax = ', impH(ibgn0) , impH(nmax1) 

        ! Now for L2 data
        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 ( * , * ) 'ibgn1 =', ibgn1 
        WRITE( *, * )'L2 starting point',ibgn1,impH(ibgn1)
        ! 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
         
        ! Determine cut-off limits
        IF ( ocd < 0 ) THEN     ! setting occultation
          imin = n + 1 - nmax1
          imax = n + 1 - ibgn0  
                 
          ro_data%Lev1a%snr_L2p(imin:imax) = 10.0_wp * snrL2(nmax1:ibgn0:ocd)
          ro_data%Lev1a%phase_L2(imin:imax) = phL2(nmax1:ibgn0:ocd)

          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   
          
          ro_data%Lev1a%snr_L2p(imin:imax) = 10.0_wp * snrL2(ibgn0:nmax1)
          ro_data%Lev1a%phase_L2(imin:imax) = phL2(ibgn0:nmax1)
               
          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(PH)
    DEALLOCATE(snrL1)
    DEALLOCATE(snrL2)
    DEALLOCATE(phL1)
    DEALLOCATE(phL2)
    DEALLOCATE(imPh)
    DEALLOCATE(time)
    DEALLOCATE(ampL1sm)
    DEALLOCATE(ampL1sm1)
    DEALLOCATE(ampL2sm)
    DEALLOCATE(dphLL)
    DEALLOCATE(dphLLsm)
 
  ELSE
  ! 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