! $Id: phs2dop.f90 2022-11-02 09:51:28Z ychen $

!****s* fsi/phs2dop *
!
! NAME
!   phs2dop  - doppler from phase 
!
! SYNOPSIS
!    call phs2dop(phase, dop)
!
! DESCRIPTION
!    This subroutine reconstruct the doppler phase 
!
! INPUTS
!    REAL(wp), DIMENSION(:,:),  :: phase     Excess phase (m) [ch, t]
!
! OUTPUT
!    REAL(wp), DIMENSION(:,:),  :: dop       smoothed doppler shift (m) [ch, t]
!
! AUTHOR
!   Yong Chen, NOAA/NESDIS/STAR, yong.chen@noaa.gov
!   Loknath Adhikari, UMD/CISESS, loknath.adhikari@noaa.gov
!
! COPYRIGHT
!   Copyright (c) 2022-2023 Yong Chen <yong.chen@noaa.gov>
!   For further details please refer to the file COPYRIGHT
!   which you should have received as part of this distribution.
!
!****
                            
SUBROUTINE phs2dop(phase, dop) 

      USE typesizes,         ONLY: wp => EightByteReal
      USE ropp_pp_constants, ONLY: c_light, pi, f_L1, f_L2
      
      IMPLICIT none 
      REAL(wp), DIMENSION(:,:), INTENT(in)    :: phase    ! Excess phase (m) [ch, t]
      REAL(wp), DIMENSION(:,:), INTENT(inout) :: dop      ! doppler shift (m) [ch, t]
       
      INTEGER  :: n, nc
      INTEGER  :: ic, i
      REAL(wp), DIMENSION(:),   ALLOCATABLE :: k     ! Wave number [channel]
        
      REAL(wp) :: dp1, dp2    
      
      ! Determination of data sizes

      n  = SIZE(phase, 2)   ! time dimension
      nc = SIZE(phase, 1)   ! 2 frequency channels
                                                                        
      ALLOCATE(k(nc))
      
      ! 5.1 Computation of wave vectors

      k(1)   = 2.0_wp*pi*f_L1/C_Light  !33.018361644831181 
      k(2)   = 2.0_wp*pi*f_L2/C_Light  !25.728593489478847
      
      DO ic=1,nc
         DO i=2,n
!            dop(ic,i) = k(ic) * (phase(ic,i) - phase(ic,i-1))
            dop(ic,i) = (phase(ic,i) - phase(ic,i-1))
          
            dp1 = dop(ic,i) + 2.0_wp*pi
            dp2 = dop(ic,i) - 2.0_wp*pi
            
            IF ( dabs(dp1) < dabs(dop(ic, i)) ) dop(ic, i) = dp1 
            IF ( dabs(dp2) < dabs(dop(ic, i)) ) dop(ic, i) = dp2 
             
         ENDDO
         dop(ic, 1) = dop(ic, 2)
!         dop(ic, :) = dop(ic, :)/k(ic)
         dop(ic, :) = dop(ic, :)
      ENDDO

      DEALLOCATE(k)
                                                                        
                                                                        
END SUBROUTINE phs2dop