1c1
< ! $Id: ropp_io_read_ncdf_get.f90 6490 2020-09-10 17:36:27Z idculv $
---
> ! $Id$
47,50c47
<                            ThisFmtVer,    &
<                            PCD_summary,   &
<                            PCD_bangle,    &
<                            PCD_refrac
---
>                            ThisFmtVer
67,69c64,65
<   INTEGER                      :: status, varid, ndim, TYPE
<   CHARACTER(len = 256)         :: routine
<   CHARACTER(len = 256)         :: readstr
---
>   INTEGER :: status, varid, ndim, TYPE
>   CHARACTER(len = 256) :: routine
106d101
<   data%bad = ' '               ; CALL ncdf_getatt('bad',  data%bad)
123,124d117
<   CALL ncdf_getvar('occfreq1', data%occfreq1, rec = irec)
<   CALL ncdf_getvar('occfreq2', data%occfreq2, rec = irec)
178c171
<   CALL Date_and_Time_UTC(Values=DT8)
---
>   CALL Date_and_Time_UTC(Values=DT8) ! Set DT8 array to the current local time
239,254d231
<   ! Yong Chen on 05/24/2024, to reset pcd and overall qual based on
<   ! the postprocess quality flag 'bad' for RFSI and ROPP
< 
<   IF(ncdf_isatt('bad')) THEN
<     readstr = ' '
<     CALL ncdf_getatt('bad', readstr)
<     IF ( TRIM(readstr(1:1)) /= "0" ) THEN
<       data%pcd = IBSET(data%pcd, PCD_summary)      ! PCD bit 1 set: non-nominal quality
<       data%pcd = IBSET(data%pcd, PCD_bangle)       ! PCD bit 5 set: non-nominal bangle
<       data%pcd = IBSET(data%pcd, PCD_refrac)       ! PCD bit 6 set: non-nominal refrac
<       data%overall_qual = 0.0
<     ELSE
<       data%overall_qual = 100.0
<     ENDIF
<   ENDIF
< 
275,280d251
<   ! Yong Chen on 05/30/2024 set Undulation to zero if not in the range
<   !write(*, *) 'Undulation range is: ',  data%georef%range%undulation
<   IF ( data%georef%Undulation <  data%georef%range%undulation(1) .or. &
<        data%georef%Undulation >  data%georef%range%undulation(2) )    &
<        data%georef%Undulation = 0.0
< 
289,296d259
<   CALL ncdf_getvar('r_alt',       data%georef%r_alt,             &
<                           units = data%georef%units%r_alt,       &
<                           range = data%georef%range%r_alt,       &
<                           rec   = irec)
<   CALL ncdf_getvar('viewang',     data%georef%viewang,           &
<                           units = data%georef%units%viewang,     &
<                           range = data%georef%range%viewang,     &
<                           rec   = irec)
334,335c297,319
< ! 1.9 Background characterisation (if any)
< ! ----------------------------------------
---
> ! 1.9 Signal characteristics (if present)
> ! ---------------------------------------
> 
>   IF ( ncdf_isvar('freq_L1') ) THEN
>     CALL ncdf_getvar('freq_L1', data%signal1%freq,       &
>                         units = data%signal1%units%freq, &
>                         range = data%signal1%range%freq, &
>                         rec   = irec)
>   ELSE
>     data%signal1%freq = 1.57542E9_wp ! GPS L1
>   ENDIF
> 
>   IF ( ncdf_isvar('freq_L2') ) THEN
>     CALL ncdf_getvar('freq_L2', data%signal2%freq,       &
>                         units = data%signal2%units%freq, &
>                         range = data%signal2%range%freq, &
>                         rec   = irec)
>   ELSE
>     data%signal2%freq = 1.2276E9_wp  ! GPS L2
>   ENDIF
> 
> ! 1.10 Background characterisation (if any)
> ! -----------------------------------------
371,424c355,356
< ! 1.9.1 QA characterisation (if any)
< ! ----------------------------------------
< 
<   CALL ncdf_getvar('snr1avg',     data%QA%snr1avg,           &
<                           units = data%QA%units%snr1avg,     &
<                           range = data%QA%range%snr1avg,     &
<                           rec   = irec)
<   CALL ncdf_getvar('snr2avg',     data%QA%snr2avg,           &
<                           units = data%QA%units%snr2avg,     &
<                           range = data%QA%range%snr2avg,     &
<                           rec   = irec)
<   CALL ncdf_getvar('irs',         data%QA%irs,               &
<                           units = data%QA%units%irs,         &
<                           range = data%QA%range%irs,         &
<                           rec   = irec)
<   CALL ncdf_getvar('smean',       data%QA%smean,             &
<                           units = data%QA%units%smean,       &
<                           range = data%QA%range%smean,       &
<                           rec   = irec)
<   CALL ncdf_getvar('stdv',        data%QA%stdv,              &
<                           units = data%QA%units%stdv,        &
<                           range = data%QA%range%stdv,        &
<                           rec   = irec)
<   CALL ncdf_getvar('reldevmax',   data%QA%reldevmax,         &
<                           units = data%QA%units%reldevmax,   &
<                           range = data%QA%range%reldevmax,   &
<                           rec   = irec)
<   CALL ncdf_getvar('reldevstd',   data%QA%reldevstd,         &
<                           units = data%QA%units%reldevstd,   &
<                           range = data%QA%range%reldevstd,   &
<                           rec   = irec)
<   CALL ncdf_getvar('reldevmax1',  data%QA%reldevmax1,        & 
<                           units = data%QA%units%reldevmax1,  & 
<                           range = data%QA%range%reldevmax1,  & 
<                           rec   = irec)
<   CALL ncdf_getvar('reldevstd1',   data%QA%reldevstd1,       & 
<                           units = data%QA%units%reldevstd1,  & 
<                           range = data%QA%range%reldevstd1,  & 
<                           rec   = irec)
<   CALL ncdf_getvar('difmaxref',   data%QA%difmaxref,         &
<                           units = data%QA%units%difmaxref,   &
<                           range = data%QA%range%difmaxref,   &
<                           rec   = irec)
<   CALL ncdf_getvar('maxdifphase', data%QA%maxdifphase,       &
<                           units = data%QA%units%maxdifphase, &
<                           range = data%QA%range%maxdifphase, &
<                           rec   = irec)
<   CALL ncdf_getvar('finalQA',     data%QA%finalQA,           &
<                           units = data%QA%units%finalQA,     &
<                           range = data%QA%range%finalQA,     &
<                           rec   = irec)
< 
< ! 1.10 Level1a variables (if any)
< ! ------------------------------
---
> ! 1.11 Level1a variables (if any)
> ! -------------------------------
433,434d364
<   WRITE(*, *) 'data%Lev1a%Npoints before reading Level1a', data%Lev1a%Npoints
< 
480c410
< ! 1.10.1 Other attributes
---
> ! 1.11.1 Other attributes
490c420
< ! 1.11 Level1b variables (if any)
---
> ! 1.12 Level1b variables (if any)
501d430
< 
518a448
>     print*,'line 449, before impact_L2'
527,529c457,458
<     ! Yong Chen on 05/24/2024, replace impact with impact_opt
<     IF (ncdf_isvar('impact_opt')) & ! added at v1.1
<      CALL ncdf_getvar('impact_opt',      data%Lev1b%impact_opt,         &
---
>   IF (ncdf_isvar('impact_opt')) & ! added at v1.1
>     CALL ncdf_getvar('impact_opt',       data%Lev1b%impact_opt,         &
546,548c475,476
<     ! Yong Chen on 05/24/2024, replace bangle with bangle_opt
<     IF (ncdf_isvar('bangle_opt')) & ! added at v1.1
<      CALL ncdf_getvar('bangle_opt',      data%Lev1b%bangle_opt,         &
---
>   IF (ncdf_isvar('bangle_opt')) & ! added at v1.1
>     CALL ncdf_getvar('bangle_opt',       data%Lev1b%bangle_opt,         &
565,567c493,494
<     ! Yong Chen on 05/24/2024, replace bangle_sigma with bangle_L1_sigma
<     IF (ncdf_isvar('bangle_opt_sigma')) & ! added at v1.1
<      CALL ncdf_getvar('bangle_opt_sigma', data%Lev1b%bangle_opt_sigma,   &
---
>   IF (ncdf_isvar('bangle_opt_sigma')) & ! added at v1.1
>     CALL ncdf_getvar('bangle_opt_sigma', data%Lev1b%bangle_opt_sigma,   &
592c519
< ! 1.12 Level2a variables (if any)
---
> ! 1.13 Level2a variables (if any)
641c568
< ! 1.13 Level2b variables (if any)
---
> ! 1.14 Level2b variables (if any)
693c620
< ! 1.14 Level2c variables (if any)
---
> ! 1.15 Level2c variables (if any)
1102c1029
< ! 1.15 Level2d variables (if any)
---
> ! 1.16 Level2d variables (if any)
1125a1053
> 
1128c1056
< ! 1.16 Level2e variables (if any)
---
> ! 1.17 Level2e variables (if any)
1216c1144
< ! 1.17 Additional variables (if any)
---
> ! 1.18 Additional variables (if any)
1245c1173
< ! 1.17 Clean up
---
> ! 1.19 Clean up
1324d1251
<   data%bad = ' '               ; CALL ncdf_getatt('bad',  data%bad)
1341,1342d1267
<   CALL ncdf_getvar('occfreq1', data%occfreq1, rec = irec)
<   CALL ncdf_getvar('occfreq2', data%occfreq2, rec = irec)
1477,1482d1401
<   ! Yong Chen on 05/30/2024 set Undulation to zero if not in the range
<   !write(*, *) 'Undulation range is: ',  data%georef%range%undulation
<   IF ( data%georef%Undulation <  data%georef%range%undulation(1) .or. &
<        data%georef%Undulation >  data%georef%range%undulation(2) )    &
<        data%georef%Undulation = 0.0
< 
1491,1498d1409
<   CALL ncdf_getvar('r_alt',       data%georef%r_alt,             &
<                           units = data%georef%units%r_alt,       &
<                           range = data%georef%range%r_alt,       &
<                           rec   = irec)
<   CALL ncdf_getvar('viewang',     data%georef%viewang,           &
<                           units = data%georef%units%viewang,     &
<                           range = data%georef%range%viewang,     &
<                           rec   = irec)
1535a1447,1448
> ! 2.9 Signal characteristics (if present)
> ! ---------------------------------------
1537,1538c1450,1469
< ! 2.9 Background characterisation (if any)
< ! ----------------------------------------
---
>   IF ( ncdf_isvar('freq_L1') ) THEN
>     CALL ncdf_getvar('freq_L1', data%signal1%freq,       &
>                         units = data%signal1%units%freq, &
>                         range = data%signal1%range%freq, &
>                         rec   = irec)
>   ELSE
>     data%signal1%freq = 1.57542E9_wp ! GPS L1
>   ENDIF
> 
>   IF ( ncdf_isvar('freq_L2') ) THEN
>     CALL ncdf_getvar('freq_L2', data%signal2%freq,       &
>                         units = data%signal2%units%freq, &
>                         range = data%signal2%range%freq, &
>                         rec   = irec)
>   ELSE
>     data%signal2%freq = 1.2276E9_wp  ! GPS L2
>   ENDIF
> 
> ! 2.10 Background characterisation (if any)
> ! -----------------------------------------
1574,1627c1505,1506
< ! 2.9.1 QA characterisation
< ! ----------------------------------------
< 
<   CALL ncdf_getvar('snr1avg',     data%QA%snr1avg,           &
<                           units = data%QA%units%snr1avg,     &
<                           range = data%QA%range%snr1avg,     &
<                           rec   = irec)
<   CALL ncdf_getvar('snr2avg',     data%QA%snr2avg,           &
<                           units = data%QA%units%snr2avg,     &
<                           range = data%QA%range%snr2avg,     &
<                           rec   = irec)
<   CALL ncdf_getvar('irs',         data%QA%irs,               &
<                           units = data%QA%units%irs,         &
<                           range = data%QA%range%irs,         &
<                           rec   = irec)
<   CALL ncdf_getvar('smean',       data%QA%smean,             &
<                           units = data%QA%units%smean,       &
<                           range = data%QA%range%smean,       &
<                           rec   = irec)
<   CALL ncdf_getvar('stdv',        data%QA%stdv,              &
<                           units = data%QA%units%stdv,        &
<                           range = data%QA%range%stdv,        &
<                           rec   = irec)
<   CALL ncdf_getvar('reldevmax',   data%QA%reldevmax,         &
<                           units = data%QA%units%reldevmax,   &
<                           range = data%QA%range%reldevmax,   &
<                           rec   = irec)
<   CALL ncdf_getvar('reldevstd',   data%QA%reldevstd,         &
<                           units = data%QA%units%reldevstd,   &
<                           range = data%QA%range%reldevstd,   &
<                           rec   = irec)
<   CALL ncdf_getvar('reldevmax1',  data%QA%reldevmax1,        & 
<                           units = data%QA%units%reldevmax1,  & 
<                           range = data%QA%range%reldevmax1,  & 
<                           rec   = irec)
<   CALL ncdf_getvar('reldevstd1',  data%QA%reldevstd1,        & 
<                           units = data%QA%units%reldevstd1,  & 
<                           range = data%QA%range%reldevstd1,  & 
<                           rec   = irec)
<   CALL ncdf_getvar('difmaxref',   data%QA%difmaxref,         &
<                           units = data%QA%units%difmaxref,   &
<                           range = data%QA%range%difmaxref,   &
<                           rec   = irec)
<   CALL ncdf_getvar('maxdifphase', data%QA%maxdifphase,       &
<                           units = data%QA%units%maxdifphase, &
<                           range = data%QA%range%maxdifphase, &
<                           rec   = irec)
<   CALL ncdf_getvar('finalQA',     data%QA%finalQA,           &
<                           units = data%QA%units%finalQA,     &
<                           range = data%QA%range%finalQA,     &
<                           rec   = irec)
< 
< ! 2.10 Level1a variables (if any)
< ! ------------------------------
---
> ! 2.11 Level1a variables (if any)
> ! -------------------------------
1683c1562
< ! 2.10.1 Other attributes
---
> ! 2.11.1 Other attributes
1692c1571
< ! 2.11 Level1b variables (if any)
---
> ! 2.12 Level1b variables (if any)
1790c1669
< ! 2.12 Level2a variables (if any)
---
> ! 2.13 Level2a variables (if any)
1842c1721
< ! 2.13 Level2b variables (if any)
---
> ! 2.14 Level2b variables (if any)
1896c1775
< ! 2.14 Level2c variables (if any)
---
> ! 2.15 Level2c variables (if any)
1943c1822
< ! 2.15 Level2d variables (if any)
---
> ! 2.16 Level2d variables (if any)
1967c1846
< ! 2.16 Additional variables (if any)
---
> ! 2.17 Additional variables (if any)
1997c1876
< ! 2.17 Clean up
---
> ! 2.18 Clean up
2187,2189c2066
<                                             lgetiono,ldummy)
<          write(*, *) 'data%Lev1a%Npoints after ropp_io_read_ncdf_get_eumdata', DATA%Lev1a%Npoints
< 
---
>                                             lgetiono, ldummy)
2366,2370d2242
<   IF ( ncdf_isatt('occfreq1') ) THEN
<     CALL ncdf_getatt('occfreq1', data%occfreq1)
<     CALL ncdf_getatt('occfreq2', data%occfreq2)
<   ENDIF
< 
2657d2528
<   data%bad = ' '               ; data%bad = "0"
2715d2585
<   INTEGER                      :: nn
2718c2588,2590
<   CHARACTER (len = 256)        :: readstr
---
>   REAL(wp), DIMENSION(2)       :: frequencies      ! Temporary frequency array
> 
>   CHARACTER(len = 256)         :: readstr
2724,2740d2595
<   REAL(wp), DIMENSION(:), ALLOCATABLE :: startTime
<   REAL(wp), DIMENSION(:), ALLOCATABLE :: fulltime
<   REAL(wp), DIMENSION(:), ALLOCATABLE :: orbtime
<   REAL(wp), DIMENSION(:), ALLOCATABLE :: txmitLR
<   REAL(wp), DIMENSION(:), ALLOCATABLE :: txmitHR
<   REAL(wp), DIMENSION(:), ALLOCATABLE :: xLR
<   REAL(wp), DIMENSION(:), ALLOCATABLE :: yLR
<   REAL(wp), DIMENSION(:), ALLOCATABLE :: zLR
<   REAL(wp), DIMENSION(:), ALLOCATABLE :: xHR
<   REAL(wp), DIMENSION(:), ALLOCATABLE :: yHR
<   REAL(wp), DIMENSION(:), ALLOCATABLE :: zHR
<   REAL(wp), DIMENSION(:), ALLOCATABLE :: d2txmit  ! 2nd derivative ln(N)
<   REAL(wp), DIMENSION(:), ALLOCATABLE :: d2xLR    ! 2nd derivative ln(N)
<   REAL(wp), DIMENSION(:), ALLOCATABLE :: d2yLR    ! 2nd derivative ln(N)
<   REAL(wp), DIMENSION(:), ALLOCATABLE :: d2zLR    ! 2nd derivative ln(N)
<   INTEGER                             :: ii      ! Index
< 
2775,2781d2629
<   IF ( ncdf_isatt('occfreq1') ) THEN
<     CALL ncdf_getatt('occfreq1', data%occfreq1)
<     CALL ncdf_getatt('occfreq2', data%occfreq2)
<   ENDIF
< 
<   write(*, *) data%occfreq1, data%occfreq2
< 
2787,2788d2634
<   write(*, *) data%DTocc%year
< 
2817d2662
<     write(*, *) n
2823,2824d2667
<   write (*, *) 'before level1a data:', data%Lev1a%Npoints
<   
2832a2676,2701
>     CALL ncdf_getvar('xLeo', data%Lev1a%r_leo(:,1))
>     CALL ncdf_getvar('yLeo', data%Lev1a%r_leo(:,2))
>     CALL ncdf_getvar('zLeo', data%Lev1a%r_leo(:,3))
>     data%Lev1a%r_leo(:,:)  = data%Lev1a%r_leo(:,:) * 1000.0_wp
>     data%Lev1a%units%r_leo = "metres"
>     data%Lev1a%reference_frame%r_leo = "ECI"
> 
>     CALL ncdf_getvar('xdLeo', data%Lev1a%v_leo(:,1))
>     CALL ncdf_getvar('ydLeo', data%Lev1a%v_leo(:,2))
>     CALL ncdf_getvar('zdLeo', data%Lev1a%v_leo(:,3))
>     data%Lev1a%v_leo(:,:)  = data%Lev1a%v_leo(:,:) * 1000.0_wp
>     data%Lev1a%units%v_leo = "metres / seconds"
> 
>     CALL ncdf_getvar('xGps', data%Lev1a%r_gns(:,1))
>     CALL ncdf_getvar('yGps', data%Lev1a%r_gns(:,2))
>     CALL ncdf_getvar('zGps', data%Lev1a%r_gns(:,3))
>     data%Lev1a%r_gns(:,:)  = data%Lev1a%r_gns(:,:) * 1000.0_wp
>     data%Lev1a%units%r_gns = "metres"
>     data%Lev1a%reference_frame%r_gns = "ECI"
> 
>     CALL ncdf_getvar('xdGps', data%Lev1a%v_gns(:,1))
>     CALL ncdf_getvar('ydGps', data%Lev1a%v_gns(:,2))
>     CALL ncdf_getvar('zdGps', data%Lev1a%v_gns(:,3))
>     data%Lev1a%v_gns(:,:)  = data%Lev1a%v_gns(:,:) * 1000.0_wp
>     data%Lev1a%units%v_gns = "metres / seconds"
> 
2837,3002d2705
<     ! Updated by Yong Chen on 03/20/2023
<     ! Working for both atmPhs and conPhs format
<     ! for new format
<     IF (ncdf_isvar('orbtime')) THEN
< 
<       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
<       ! added by Xinjia, for new UCAR conPhs format
<       CALL ncdf_getsize('orbtime', nn, dim = 1)
<       ALLOCATE(orbtime(nn))
<       ALLOCATE(txmitLR(nn))
< 
<       CALL ncdf_getvar('orbtime', orbtime)
<       CALL ncdf_getvar('txmitLR', txmitLR)
< 
<       ALLOCATE(startTime(1))
<       CALL ncdf_getvar('startTime', startTime)
< 
<       ALLOCATE(fulltime(n))
<       fulltime(:) = startTime(1)+data%Lev1a%dtime
< 
<       ALLOCATE(xLR(nn))
<       ALLOCATE(yLR(nn))
<       ALLOCATE(zLR(nn))
< 
<       ALLOCATE(d2xLR(nn))
<       ALLOCATE(d2yLR(nn))
<       ALLOCATE(d2zLR(nn))
< 
<       ALLOCATE(xHR(n))
<       ALLOCATE(yHR(n))
<       ALLOCATE(zHR(n))
< 
<       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
<       CALL ncdf_getvar('xLeoLR', xLR)
<       CALL ncdf_getvar('yLeoLR', yLR)
<       CALL ncdf_getvar('zLeoLR', zLR)
< 
<       CALL ropp_io_init_spline(orbtime(:), xLR(:), d2xLR)
<       CALL ropp_io_init_spline(orbtime(:), yLR(:), d2yLR)
<       CALL ropp_io_init_spline(orbtime(:), zLR(:), d2zLR)
< 
<       DO ii = 1, n
<         CALL ropp_io_interpol_spline(orbtime, xLR, d2xLR, &
<               fulltime(ii), xHR(ii))
<         CALL ropp_io_interpol_spline(orbtime, yLR, d2yLR, &
<               fulltime(ii), yHR(ii))
<         CALL ropp_io_interpol_spline(orbtime, zLR, d2zLR, &
<               fulltime(ii), zHR(ii))
<       ENDDO
<       data%Lev1a%r_leo(:,1) = xHR(:)
<       data%Lev1a%r_leo(:,2) = yHR(:)
<       data%Lev1a%r_leo(:,3) = zHR(:)
< 
<       data%Lev1a%r_leo(:,:)  = data%Lev1a%r_leo(:,:) * 1000.0_wp
<       data%Lev1a%units%r_leo = "metres"
<       data%Lev1a%reference_frame%r_leo = "ECI"
< 
<       DEALLOCATE(xLR)
<       DEALLOCATE(yLR)
<       DEALLOCATE(zLR)
< 
<       DEALLOCATE(d2xLR)
<       DEALLOCATE(d2yLR)
<       DEALLOCATE(d2zLR)
< 
<       DEALLOCATE(xHR)
<       DEALLOCATE(yHR)
<       DEALLOCATE(zHR)
<       ! updating stopped here
<       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
< 
<       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
<       !data%Lev1a%v_leo(:,:)  = data%Lev1a%r_leo(:,:)
<       data%Lev1a%units%v_leo = "metres / seconds"
< 
<       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
<       ALLOCATE(xLR(nn))
<       ALLOCATE(yLR(nn))
<       ALLOCATE(zLR(nn))
< 
<       ALLOCATE(d2xLR(nn))
<       ALLOCATE(d2yLR(nn))
<       ALLOCATE(d2zLR(nn))
< 
<       ALLOCATE(xHR(n))
<       ALLOCATE(yHR(n))
<       ALLOCATE(zHR(n))
< 
<       ALLOCATE(txmitHR(n))
< 
<       CALL ropp_io_interpol(orbtime, fulltime, txmitLR, txmitHR)
< 
<       CALL ncdf_getvar('xGnssLR', xLR)
<       CALL ncdf_getvar('yGnssLR', yLR)
<       CALL ncdf_getvar('zGnssLR', zLR)
< 
<       CALL ropp_io_init_spline(txmitLR(:), xLR(:), d2xLR)
<       CALL ropp_io_init_spline(txmitLR(:), yLR(:), d2yLR)
<       CALL ropp_io_init_spline(txmitLR(:), zLR(:), d2zLR)
< 
<       DO ii = 1, n
<         CALL ropp_io_interpol_spline(txmitLR, xLR, d2xLR, &
<               txmitHR(ii), xHR(ii))
<         CALL ropp_io_interpol_spline(txmitLR, yLR, d2yLR, &
<               txmitHR(ii), yHR(ii))
<         CALL ropp_io_interpol_spline(txmitLR, zLR, d2zLR, &
<               txmitHR(ii), zHR(ii))
< !        print *, ii, xHR(ii), yHR(ii), zHR(ii)
<       ENDDO
<       data%Lev1a%r_gns(:,1) = xHR(:)
<       data%Lev1a%r_gns(:,2) = yHR(:)
<       data%Lev1a%r_gns(:,3) = zHR(:)
< 
<       data%Lev1a%r_gns(:,:)  = data%Lev1a%r_gns(:,:) * 1000.0_wp
<       data%Lev1a%units%r_gns = "metres"
<       data%Lev1a%reference_frame%r_gns = "ECI"
< 
< !      print *, 'data%Lev1a%r_gns*****************'
< !      print *, data%Lev1a%r_gns
< 
<       DEALLOCATE(xLR)
<       DEALLOCATE(yLR)
<       DEALLOCATE(zLR)
< 
<       DEALLOCATE(d2xLR)
<       DEALLOCATE(d2yLR)
<       DEALLOCATE(d2zLR)
< 
<       DEALLOCATE(xHR)
<       DEALLOCATE(yHR)
<       DEALLOCATE(zHR)
< 
<       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
<       !data%Lev1a%v_gns(:,:)  = data%Lev1a%r_gns(:,:)
<       data%Lev1a%units%v_gns = "metres / seconds"
< 
<     ELSE      ! for old foramt
<       CALL ncdf_getvar('xLeo', data%Lev1a%r_leo(:,1))
<       CALL ncdf_getvar('yLeo', data%Lev1a%r_leo(:,2))
<       CALL ncdf_getvar('zLeo', data%Lev1a%r_leo(:,3))
<       data%Lev1a%r_leo(:,:)  = data%Lev1a%r_leo(:,:) * 1000.0_wp
<       data%Lev1a%units%r_leo = "metres"
<       data%Lev1a%reference_frame%r_leo = "ECI"
< 
<       CALL ncdf_getvar('xdLeo', data%Lev1a%v_leo(:,1))
<       CALL ncdf_getvar('ydLeo', data%Lev1a%v_leo(:,2))
<       CALL ncdf_getvar('zdLeo', data%Lev1a%v_leo(:,3))
<       data%Lev1a%v_leo(:,:)  = data%Lev1a%v_leo(:,:) * 1000.0_wp
<       data%Lev1a%units%v_leo = "metres / seconds"
< 
<       CALL ncdf_getvar('xGps', data%Lev1a%r_gns(:,1))
<       CALL ncdf_getvar('yGps', data%Lev1a%r_gns(:,2))
<       CALL ncdf_getvar('zGps', data%Lev1a%r_gns(:,3))
<       data%Lev1a%r_gns(:,:)  = data%Lev1a%r_gns(:,:) * 1000.0_wp
<       data%Lev1a%units%r_gns = "metres"
<       data%Lev1a%reference_frame%r_gns = "ECI"
< 
<       CALL ncdf_getvar('xdGps', data%Lev1a%v_gns(:,1))
<       CALL ncdf_getvar('ydGps', data%Lev1a%v_gns(:,2))
<       CALL ncdf_getvar('zdGps', data%Lev1a%v_gns(:,3))
<       data%Lev1a%v_gns(:,:)  = data%Lev1a%v_gns(:,:) * 1000.0_wp
<       data%Lev1a%units%v_gns = "metres / seconds"
< 
<     ENDIF
< 
<     
3008c2711
<       CASE ('ME','KO','PA','TS','TD','C2','GO','C0','S0','S1','S2')
---
>       CASE ( 'ME', 'KO' )
3040,3041d2742
<   write(*, *) 'end level1a data'
< 
3054d2754
<   data%bad = ' '               ; data%bad = "0"
3074a2775,2785
> ! Newer versions of atmPhs files contain frequencies.
> ! If not available, let the default be GPS L1 and L2.
>   IF (ncdf_isatt('occfreq1') .AND. ncdf_isatt('occfreq2')) THEN
>     CALL ncdf_getatt('occfreq1', frequencies(1))
>     CALL ncdf_getatt('occfreq2', frequencies(2))
>   ELSE
>     frequencies = (/1.57542E9_wp, 1.2276E9_wp/)
>   END IF
>   data%signal1%freq = frequencies(1)
>   data%signal2%freq = frequencies(2)
> 
3084,3193d2794
<   write (*,*) 'end of ropp_io_read_ucardata_atmPhs'
< END SUBROUTINE ropp_io_read_ucardata_atmPhs
< 
< !-------------------------------------------------------------------------------
< ! 15. Splines: essentially a copy of the ropp_io_spline module
< !-------------------------------------------------------------------------------
< 
< ! 15.1 Generate second derivative of spline
< ! -----------------------------------------
< 
<   SUBROUTINE ropp_io_interpol(x, newx, array, interp)
< 
<     IMPLICIT NONE
< 
<     REAL(wp), DIMENSION(:), INTENT(in)  :: x
<     REAL(wp), DIMENSION(:), INTENT(in)  :: newx
<     REAL(wp), DIMENSION(:), INTENT(in)  :: array
<     REAL(wp), DIMENSION(:), INTENT(out) :: interp
< 
<     INTEGER                             :: i, j, k
< 
<     DO k = 1, SIZE(newx)
<       j = 2
<       DO WHILE (j < SIZE(x) .AND. x(j) < newx(k))
<         j = j + 1
<       ENDDO
<       i = j - 1
< 
<       interp(k) = array(j) + &
<           ( (newx(k) - x(j)) / (x(i) - x(j)) * (array(i) - array(j)) )
< 
<     ENDDO
< 
<   END SUBROUTINE ropp_io_interpol
< 
< 
< ! 15.2 Generate second derivative of spline
< ! -----------------------------------------
< 
<   SUBROUTINE ropp_io_init_spline(x, f, d2)
< 
<     IMPLICIT NONE
< 
<     ! 15.2.1 Declarations
< 
<     REAL(wp), DIMENSION(:), INTENT(in)  :: x  ! Argument grid (monotonous)
<     REAL(wp), DIMENSION(:), INTENT(in)  :: f  ! Gridded function
<     REAL(wp), DIMENSION(:), INTENT(out) :: d2 ! 2nd derivative of spline
< 
<     REAL(wp), DIMENSION(:), ALLOCATABLE :: d1
<     REAL(wp)                     :: dfl, dfr, df, a, b, c
<     INTEGER                      :: i, N
< 
<     ! 15.2.2 Initialisation
< 
<     N = SIZE(x)
<     ALLOCATE(d1(N))
<     d2(:) = 0.0_wp
< 
<     ! 2.3 Drive-through calculation of spline coefficients
< 
<     d1(1) = 0.0_wp
<     d2(N) = 0.0_wp
<     d1(1) = 0.0_wp
< 
<     DO i = 2, N-1
< 
<        dfl   = (f(i) - f(i-1))/(x(i) - x(i-1))
<        dfr   = (f(i+1) - f(i))/(x(i+1) - x(i))
<        df    = (dfr - dfl)/(x(i+1) - x(i-1))
<        a     = (x(i) - x(i-1))/(2*(x(i+1) - x(i-1)))
<        b     = (x(i+1) - x(i))/(2*(x(i+1) - x(i-1)))
<        c     = 1 + a*d1(i-1)
<        d1(i) = -b/c
<        d2(i) = (3*df - a*d2(i-1))/c
<     ENDDO
< 
<     DO i = N-1, 2, -1
<        d2(i) = d1(i)*d2(i+1) + d2(i)
<     ENDDO
< 
<     DEALLOCATE(d1)
< 
<   END SUBROUTINE ropp_io_init_spline
< 
< ! 15.3 Return interpolated function and first two derivatives
< ! -----------------------------------------------------------
< 
<   SUBROUTINE ropp_io_interpol_spline(x, f, d2, x_int, f_int, fd_int, fd2_int)
< 
<     IMPLICIT NONE
< 
<     ! 15.3.1 Declarations
< 
<     REAL(wp), DIMENSION(:), INTENT(in)  :: x  ! Argument grid (monotonous)
<     REAL(wp), DIMENSION(:), INTENT(in)  :: f  ! Gridded function
<     REAL(wp), DIMENSION(:), INTENT(in)  :: d2 ! 2nd derivative of spline
<     REAL(wp),               INTENT(in)  :: x_int  ! Interpolation point
<     REAL(wp),               INTENT(out) :: f_int  ! Interpolated function value
<     REAL(wp), OPTIONAL,     INTENT(out) :: fd_int ! Interpolated 1st derivative
<     REAL(wp), OPTIONAL,     INTENT(out) :: fd2_int ! Interpolated 2nd deriv
< 
<     REAL(wp) :: a1, a2, a3   ! Polynomial coefficients
<     REAL(wp) :: dx           ! Grid interval
<     REAL(wp) :: dx_t         ! Grid-point to interpolation-point distance
<     REAL(wp) :: x_t          ! Interpolation point projected to grid extent
<     REAL(wp) :: fd           ! Interpolated derivative
<     INTEGER  :: i            ! Array index
<     INTEGER  :: N            ! Number of data
<     INTEGER  :: i_int        ! Interpolation interval index
3195,3298c2796
<     ! 15.3.2 Location of interpolation point inside grid
< 
<     N     = SIZE(x)
< 
<     x_t   = MIN(MAX(x_int, MIN(x(1),x(N))), MAX(x(1),x(N)))
< 
<     i_int = ropp_io_seek_index(x, x_t)
<     i     = MAX(i_int, 1)
< 
<     ! 15.3.3 Calculation of interpolation coefficients
< 
<     dx    = x(i+1) - x(i)
<     a2    = d2(i)/2.0_wp
<     a3    = (d2(i+1) - d2(i))/(6.0_wp*dx)
<     a1    = (f(i+1) - f(i))/dx - dx*(a2 + dx*a3)
< 
<     ! 15.3.4 Calculated interpolated value
< 
<     dx_t  = x_t - x(i)
<     fd    = a1 + dx_t*(2*a2 + dx_t*3*a3)
<     f_int = f(i) + dx_t*(a1 + dx_t*(a2 + dx_t*a3)) + &
<          fd*(x_int - x_t)
< 
<     IF (PRESENT(fd_int)) THEN
<        fd_int = fd
<     ENDIF
< 
<     IF (PRESENT(fd2_int)) THEN
<        fd2_int  = 2.0_wp*a2 + dx_t*6.0_wp*a3
<     ENDIF
< 
<   END SUBROUTINE ropp_io_interpol_spline
< 
< ! 15.4 Return index of point in grid
< ! ----------------------------------
< 
<   FUNCTION ropp_io_seek_index(x, xp) RESULT(ip)
< 
<     IMPLICIT NONE
< 
<     REAL(wp), DIMENSION(:), INTENT(in) :: x   ! x-grid (homogeneous)
<     REAL(wp),               INTENT(in) :: xp  ! point to locate inside grid
<     INTEGER                            :: ip
< 
<     INTEGER  :: N      ! Number of grid points
<     INTEGER  :: Imin   ! Upper estimate of index
<     INTEGER  :: Imax   ! Lower estimate of index
<     INTEGER  :: Dir    ! Direction of argument change
<     INTEGER  :: It     ! Iteration count
<     INTEGER  :: di     ! Index increment in iterations
<     INTEGER  :: is     ! Step direction count
< 
<     ! 15.4.1 Grid size and direction calculation
< 
<     N    = SIZE(x)
<     Dir  = NINT(SIGN(1.0_wp, x(N)-x(1)))
< 
<     ! 15.4.2 Checking if point is inside grid
< 
<     IF ((Dir*xp < Dir*x(1)) .OR. (Dir*xp > Dir*x(N))) THEN
<        ip = 0
<        RETURN
<     END IF
< 
<     ! 15.4.3 Initial approximation
< 
<     Imin = 1
<     Imax = N
<     ip   = Imin+FLOOR(REAL(Imax - Imin, wp)*(xp - x(Imin))/(x(Imax) - x(Imin)))
<     ip   = MAX(1, MIN(ip, N-1))
< 
<     ! 15.4.4 Iterative index search
< 
<     It = 0
<     is = 0
< 
<     Search: DO
<        IF ((Dir*x(ip) <= Dir*xp) .AND. (Dir*xp <= Dir*x(ip+1))) THEN
<           EXIT Search
<        END IF
<        IF (ABS(is) > 1) THEN
<           ip = (Imax + Imin)/2
<           is = 0
<        END IF
<        IF (Dir*x(ip+1) < Dir*xp) THEN
<           Imin = ip + 1
<           di   = FLOOR(REAL(Imax - Imin, wp)*(xp - x(Imin))/(x(Imax) - x(Imin)))
<           ip   = Imin + di
<           is   = is + 1
<        ELSE IF (Dir*xp < Dir*x(ip)) THEN
<           Imax = ip
<           di   = FLOOR(REAL(Imin - Imax, wp)*(xp - x(Imax))/(x(Imin) - x(Imax)))
<           ip   = Imax + di
<           is   = is - 1
<        END IF
<        ip = MAX(1, MIN(ip, N-1))
<        It = It + 1
<        IF (It > N) THEN
<           ip = -1
<           EXIT Search
<        END IF
<     END DO Search
< 
<   END FUNCTION ropp_io_seek_index
---
> END SUBROUTINE ropp_io_read_ucardata_atmPhs
3342c2840
<                                          getiono,ldummy)
---
>                                          getiono, ldummy)
3375,3376d2872
<   write(*, *) 'data%Lev1a%Npoints after ropp_io_read_eumdata', DATA%Lev1a%Npoints
< 
3433,3434c2929,2930
<   REAL(wp), DIMENSION(:,:), ALLOCATABLE :: r_gns, v_gns        ! Temporary pos/velocity arrays
<   REAL(wp), DIMENSION(:,:), ALLOCATABLE :: r_leo, v_leo        ! Temporary pos/velocity arrays
---
>   REAL(wp), DIMENSION(:,:), ALLOCATABLE :: r_gns, v_gns          ! Temporary pos/velocity arrays
>   REAL(wp), DIMENSION(:,:), ALLOCATABLE :: r_leo, v_leo          ! Temporary pos/velocity arrays
3439,3443c2935,2939
<   REAL(wp), DIMENSION(:),   ALLOCATABLE :: i_ca_uncorr      ! OL or RS I component [V]
<   REAL(wp), DIMENSION(:),   ALLOCATABLE :: q_ca_uncorr      ! OL or RS Q component [V]
<   REAL(wp), DIMENSION(:),   ALLOCATABLE :: exphase_l1_nco   ! OL or RS NCO excess phase [m]
<   REAL(wp), DIMENSION(:),   ALLOCATABLE :: phase_l1_iq      ! I/Q contribution
<   INTEGER,  DIMENSION(:),   ALLOCATABLE :: tracking_state   ! Tracking state
---
>   REAL(wp), DIMENSION(:),   ALLOCATABLE :: i_ca_uncorr           ! OL or RS I component [V]
>   REAL(wp), DIMENSION(:),   ALLOCATABLE :: q_ca_uncorr           ! OL or RS Q component [V]
>   REAL(wp), DIMENSION(:),   ALLOCATABLE :: exphase_l1_nco        ! OL or RS NCO excess phase [m]
>   REAL(wp), DIMENSION(:),   ALLOCATABLE :: phase_l1_iq           ! I/Q contribution
>   INTEGER,  DIMENSION(:),   ALLOCATABLE :: tracking_state        ! Tracking state
3452a2949,2950
> !  REAL(wp), DIMENSION(:),   ALLOCATABLE :: phase_L1_this
> !  REAL(wp), DIMENSION(:),   ALLOCATABLE :: phase_L2_this
3476,3478c2974,2976
<   INTEGER                               :: groupid, numgrps, status     ! group ID (if applicable)
<   LOGICAL                               :: havegroup   ! found group
<   INTEGER,  DIMENSION(2)                :: ncids
---
>   INTEGER                             :: groupid, numgrps, status     ! group ID (if applicable)
>   LOGICAL                             :: havegroup   ! found group
>   INTEGER,  DIMENSION(2)              :: ncids
3480d2977
<   REAL(wp), PARAMETER                   :: f_L1 = 1.57542e9_wp
3506c3003
<   CALL ncdf_getatt('spacecraft', readstr)
---
>   CALL ncdf_getatt('spacecraft', readstr)  !! global attribute, above /data/
3524c3021
<   data%FmtVersion = TRIM(fver)
---
>   data%FmtVersion = TRIM(fver)  ! '14.0' in MSG example file
3526a3024,3026
>   !print *, 'file:',file
>   !print *, 'data%leo_id:',data%leo_id
> 
3535c3035
<   CALL ncdf_getatt('/data/transmitter/satellite/satellite_prn', readstr)
---
>   CALL ncdf_getatt('/data/transmitter/satellite/satellite_prn', readstr) ! GNSS ID. By Xinjia, 05-05-2025
3566c3066,3067
<   CALL ncdf_getatt('environment', readstr)
---
>   CALL ncdf_getatt('environment', readstr) !! global attribute, above /data/
>                                            !! value is 'Development' in MSG example file
3591a3093
>   !print *, 'before 9.5 Date and time'
3594d3095
< 
3596c3097
<   
---
> 
3628a3130
>   !print *, 'before 9.6 Georeferencing'
3641d3142
< 
3652a3154
>   !print *, 'before 9.6 (velocity_rec_fixed_bfr)'
3674a3177
>   data%georef%leo_pod%vel_frame = 'ECI'
3675a3179
>   !print *, 'before 9.6 (position_gns_fixed)'
3686a3191
>   !print *, 'before 9.6 (velocity_gns)'
3762c3267
<       status = ncdf_getgroupid(ncdf_ncid, '/data/level_1a/combined/signals', &
---
>       status = ncdf_getgroupid(ncdf_ncid, '/data/level_1a/combined/gps_start_absdate', &
3808,3811d3312
<   data%occfreq1 = frequencies(1)
<   data%occfreq2 = frequencies(2)
<   !data%signal1%freq = frequencies(1)
<   !data%signal2%freq = frequencies(2)
3813c3314,3317
< ! 9.7 Level1a variables (if any)
---
>   data%signal1%freq = frequencies(1)
>   data%signal2%freq = frequencies(2)
> 
> ! 9.8 Level1a variables (if any)
3816c3320,3321
< ! 9.7.1 Closed Loop Level1a variables (if requested)
---
>   !print *, 'before 9.8 Level1a'
> ! 9.8.1 Closed Loop Level1a variables (if requested)
3819d3323
< 
4125c3629
< ! 9.7.2 Raw Sampling or Open Loop Level1a variables (if requested)
---
> ! 9.8.2 Raw Sampling or Open Loop Level1a variables (if requested)
4127a3632
>   print*, 'before 9.8.2 Raw Sampling or Open Loop Level1a' ! <--- below not run
4129,4130c3634,3635
<          (TRIM(ADJUSTL(getlevel1a)) == 'cl+ol') .OR. &
<          (TRIM(ADJUSTL(getlevel1a)) == 'cl+rs') ) THEN
---
>          ( (TRIM(ADJUSTL(getlevel1a)) == 'cl+ol') .OR. &
>            (TRIM(ADJUSTL(getlevel1a)) == 'cl+rs') ) ) THEN
4131a3637
>   print*, 'inside 9.8.2 Raw Sampling or Open Loop Level1a' ! checked. not run here
4163a3670
>       print *, 'before ncdf_getvar(/quality/)'
4170a3678
>         print *, 'line 3621, n_olrs=',n_olrs  ! 
4234c3742
< !          CALL ncdf_getvar(TRIM(tdir)//'navbits_internal', navbit_int)
---
>           CALL ncdf_getvar(TRIM(tdir)//TRIM(navbits_internal), navbit_int)
4238,4242c3746,3750
< !          IF ( have_nb > 0 ) THEN
< !            CALL ncdf_getvar(TRIM(tdir)//'navbits_external', navbit_ext)
< !          ELSE
< !            navbit_ext(:) = navbit_int(:)
< !          ENDIF
---
>           IF ( have_nb > 0 ) THEN
>             CALL ncdf_getvar(TRIM(tdir)//TRIM(navbits_external), navbit_ext)
>           ELSE
>             navbit_ext(:) = navbit_int(:)
>           ENDIF
4285c3793
<             (c_light/(2.0_wp*pi1*data%occfreq1))*phase_l1_iq(:)
---
>             (c_light/(2.0_wp*pi1*data%signal1%freq))*phase_l1_iq(:)
4476a3985
>       print*, 'line 3929, inside cl+ol OR cl+rs' ! checked. not run here
4477a3987,3990
>     print*, 'line 3931, outside cl+ol OR cl+rs' ! checked. resume from here
> 
> ! 9.8.3 Combine CL and OL or RS Level1a variables (if requested)
> ! --------------------------------------------------------------
4479,4480c3992
< ! 9.7.3 Combine CL and RS  Level1a variables (if requested)
< ! ---------------------------------------------------------
---
>   !print *, 'before 9.8.3 Combine CL and OL or RS Level1a'
4484d3995
<     write(*, *) 'n_cl, n_olrs, n', n_cl, n_olrs, n
4582c4093
< ! 9.7.4 Missing/invalid data checks
---
> ! 9.8.4 Missing/invalid data checks
4639c4150
< ! 9.8 Level1b variables (if any and if requested)
---
> ! 9.9 Level1b variables (if any and if requested)
4641a4153
>   !print *, 'before 9.9 Level1b'
4649a4162
>     print*, ' line 4106, ELSE' ! checked. <--- reach here, msg and s6a
4653a4167
>   print*, ' line 4111, getlevel1b, getiono',getlevel1b, getiono ! checked,msg and s6a <-- getlevel1b, getiono: F/T, F
4655a4170
>     print*, ' line 4114, read level1b', TRIM(ddir) ! checked,msg and s6a
4662a4178
>     print*, ' line 4122, data%Lev1b%Npoints=',data%Lev1b%Npoints ! checked,msg and s6a
4664a4181
>       print*, 'line 4178, inside Npoints > 0' ! confirmed, and impact_height added, msg and s6a
4672a4190
>       print*, 'line 4187, before impact_height' ! confirmed, msg and s6a
4676a4195
>         print*, 'line 4188, inside getiono' ! not here
4701a4221
>         print*, 'line 4211, get lat_tp' ! confirmed, msg and s6a
4704a4225,4227
>         !print*, data%Lev1b%lat_tp
> 
>         !print*, 'line 4217, get lon_tp' ! confirmed
4707a4231,4233
>         !print*, data%Lev1b%lon_tp
> 
>         !print*, 'line 4223, get bangle'  ! confirmed
4711a4238
>           print*, 'line 4228, get bangle_ca' ! not here
4714a4242
>           print*, 'line 4232, get bangle_l1'  ! confirmed, msg and s6a
4719a4248
>           print*, 'line 4238, getextrap' ! not here
4729a4259
>             print*, 'line 4249, 13.0' ! not here
4733a4264
>               print*, 'line 4253, get bangle_l5' ! confirmed, msg
4736a4268
>               print*, 'line 4266, get bangle_l2' ! confirmed, s6a
4740c4272,4275
<            ENDIF
---
> 
>             !print*, data%Lev1b%bangle_L2 ! value is correct, confirmed
> 
>           ENDIF
4751c4286,4288
<  
---
>         !print*, ' line 4213, data%Lev1b%lat_tp' ! checked, all correct <---
>         !print*, data%Lev1b%lat_tp
> 
4753c4290,4292
<  
---
>         !print*, ' line 4217, data%Lev1b%lon_tp' ! checked, all correct <---
>         !print*, data%Lev1b%lon_tp
> 
4768a4308
>         print*, 'line 4292, ol_data_available' ! confirmed
4772c4312,4313
<           CASE ("SE6A", "SE6B")
---
>           CASE ("SE6A", "SE6B", "MESG")
>             print*, 'line 4297, readbyte1' ! not here
4791a4333
>         print*, 'line 4317, ol_data_used' ! confirmed
4808a4351
>   print*, ' line 4335, after level1b'
4810,4811c4353,4354
< ! 9.9 (Global) Attributes
< ! -----------------------
---
> ! 9.10 (Global) attributes
> ! ------------------------
4827a4371
>     print*, 'line 4287, pod_method'
4828a4373
>     print*, 'line 4289, phase_method'
4829a4375
>     print*, 'line 4291, retrieval_method'
4843a4390
>         print*, 'line 4306, thinner_method'
4860a4408
>   print*, 'line 4324, /status/processing/processor_name'
4862a4411
>   print*, 'line 4324, /status/processing/processor_version'
4878a4428
>   print*, 'line 4344, before SCAN'
4881a4432
>   print*, 'line 4344, after SCAN, i/j = ', i,j
4909a4461
>   print*, 'line 4377, processing_software'
4945a4498
>     print*, 'line 4414, /status/processing/creation_time_utc'
4946a4500
>     print*, 'line 4416, after /status/processing/creation_time_utc:', readreal
4961c4515
< ! 9.10 Occultation ID
---
> ! 9.11 Occultation ID
4966c4520
< ! 9.11 Clean up
---
> ! 9.12 Clean up
4982,4983c4536
< !      Adding +-2*Pi where phase jumps from
< !      +-Pi to -+Pi,
---
> !      Adding +-2*Pi where phase jumps from +-Pi to -+Pi,
4985,4986c4538
< !      Adding +2*Pi where phase jumps from
< !      - to +
---
> !      Adding +2*Pi where phase becomes smaller than previous
4988,4989c4540
< !      Adding -2*Pi where phase jumps from
< !      + to -
---
> !      Adding -2*Pi where phase becomes larger than previous
