! $Id: ropp_pp_spectra.f90 2021 2009-01-16 10:49:04Z frhl $ !SUBROUTINE ropp_pp_spectra(time, phase_L1, phase_L2, phase_LM, & ! impact_LM, config, OutRO, filnam) ! Updated by Yong Chen, to output the straight line tangent altitude SUBROUTINE ropp_pp_spectra(time, r_leo, r_gns, r_coc, roc, & phase_L1, phase_L2, snr, phase_LM, & impact_LM, config, OutRO, filnam) !****s* Preprocessing/ropp_pp_spectra * ! ! NAME ! ropp_pp_spectra - Calculation of local spatial spectra ! ! SYNOPSIS ! call ropp_pp_spectra(time, phase_L1, phase_L2, phase_LM, impact_LM, ! config, OutRO, filnam) ! ! DESCRIPTION ! ! INPUTS ! real(wp), dimension(:) :: time ! time of samples (s) ! real(wp), dimension(:,:) :: r_leo ! cartesian LEO coordinates (m) ! real(wp), dimension(:,:) :: r_gns ! cartesian GPS coordinates (m) ! real(wp), dimension(:) :: r_coc ! cartesian centre curvature (m) ! real(wp) :: roc ! radius of curvature (m) ! real(wp), dimension(:) :: phase_L1 ! excess phase L1 (m) ! real(wp), dimension(:) :: phase_L2 ! excess phase L2 (m) ! real(wp), dimension(:) :: phase_LM ! model excess phase (m) ! real(wp), dimension(:,:) :: snr ! L1 and L2 amplitude ! real(wp), dimension(:) :: impact_LM ! model impact parameter (m) ! type(PPConfig) :: config ! Configuration options ! logical, optional :: OutRO ! Flag to output spectra results ! character(len=*), optional :: filnam ! Output file name root ! ! OUTPUT ! Spectra as functions of (doppler,time) written to output ASCII files ! and netCDF files. ! ! NOTES ! ! REFERENCES ! Gorbunov M.E., Lauritsen K.B., Rhodin A., Tomassini M. and Kornblueh L. ! 2006 ! Radio holographic filtering, error estimation, and quality control of ! radio occultation data ! Journal of Geophysical Research (111) D10105 ! ! 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_utils, ONLY: get_io_unit, impact_parameter USE ropp_io_types, ONLY: ROprof USE ropp_io USE ropp_pp, not_this => ropp_pp_spectra USE ropp_pp_types, ONLY: PPConfig USE ropp_pp_constants, ONLY: pi, f_L1, f_L2, c_light IMPLICIT NONE TYPE(ROprof) :: ro_data ! Hold spectra as RO extra data CHARACTER(LEN=256) :: text_filnam ! Name of text file CHARACTER(LEN=256) :: nc_filnam ! Name of netCDF file CHARACTER(LEN=256) :: cline ! Dummy line in spectral file INTEGER, PARAMETER :: nsamp=10 ! Stride of output nc data INTEGER :: id_nx, ix, nx ! INTEGER :: id_ny, iy, ny ! REAL(wp), DIMENSION(:, :), ALLOCATABLE :: x, y, z ! 2d fields from text file REAL(wp), DIMENSION(:), INTENT(in) :: time ! Time of samples (s) REAL(wp), DIMENSION(:,:), INTENT(in) :: r_leo ! Cartesian LEO coords (m) REAL(wp), DIMENSION(:,:), INTENT(in) :: r_gns ! Cartesian GPS coords (m) REAL(wp), DIMENSION(:), INTENT(in) :: r_coc ! Centre curvature coord (m) REAL(wp), INTENT(in) :: roc ! Radius of curvature (m) REAL(wp), DIMENSION(:), INTENT(in) :: phase_L1 ! Excess phase L1 (m) REAL(wp), DIMENSION(:), INTENT(in) :: phase_L2 ! Excess phase L2 (m) REAL(wp), DIMENSION(:,:), INTENT(in) :: snr ! L1 and L2 amplitude REAL(wp), DIMENSION(:), INTENT(in) :: phase_LM ! Model excess phase (m) REAL(wp), DIMENSION(:), INTENT(in) :: impact_LM ! Model impact parameter (m) TYPE(PPconfig), INTENT(in) :: config ! Configuration options LOGICAL, OPTIONAL, INTENT(in) :: OutRO ! Flag to output RO spectra LOGICAL :: l_out ! Output fields? INTEGER :: lun ! Unit number of output text files CHARACTER(LEN=*), OPTIONAL, INTENT(in) :: filnam ! Output file name root COMPLEX(wp), DIMENSION(:), ALLOCATABLE :: U ! Downconverted complx field COMPLEX(wp), DIMENSION(:), ALLOCATABLE :: US ! Sliding spectrum REAL(wp), DIMENSION(:), ALLOCATABLE :: dphase ! Phase deviation REAL(wp), DIMENSION(:), ALLOCATABLE :: doppler ! Doppler frequency REAL(wp), DIMENSION(:), ALLOCATABLE :: window ! Window function REAL(wp), DIMENSION(:), ALLOCATABLE :: HSL ! Straight line tangent altitude REAL(wp) :: k ! L1/L2 wave vectors INTEGER :: i, ic ! Index counters INTEGER :: is ! Sliding window index INTEGER :: is1 ! Start window position INTEGER :: isN ! End window position INTEGER :: di ! Index step INTEGER :: n ! Number of data points INTEGER :: ns ! Number of sliding windows INTEGER :: ocd ! Occulatation direction COMPLEX(wp), PARAMETER :: Ci = (0.0_wp, 1.0_wp) ! Sqrt(-1) INTEGER, PARAMETER :: nc = 2 ! Number of channels ! INTEGER, PARAMETER :: nss = 64 ! Sliding window width ! updated by Yong Chen, 08/09/2021 default value of 64 for 50 Hz Occ data, 128 for 100 Hz INTEGER :: nss = 64 ! Sliding window width REAL(wp) :: difftime ! OCC time interval LOGICAL :: dummy !------------------------------------------------------------------------------- ! 2. Initialization !------------------------------------------------------------------------------- n = SIZE(time) ALLOCATE(dphase(n)) ALLOCATE(HSL(n)) dummy = config%obs_ok ! dummy use of otherwise unused argument l_out = .FALSE. IF ( PRESENT(OutRO) ) l_out = OutRO !------------------------------------------------------------------------------- ! 3. Sliding spectral analysis !------------------------------------------------------------------------------- ! 3.1 Occultation direction (-1 = setting, +1 = rising) ocd = NINT(SIGN(1.0_wp, impact_LM(n)-impact_LM(1))) ! updated by Yong Chen, 08/09/2021, for C2 OCC rate 100Hz difftime = time(2)-time(1) if(difftime < 0.015) nss = 128 SELECT CASE (ocd) CASE (-1) is1 = 1 + nss/2 isN = n - (nss/2 - 1) di = 1 CASE (+1) is1 = n - nss/2 isN = 1 + (nss/2 - 1) di = -1 END SELECT ns = di*isN - di*is1 + 1 ALLOCATE(U(n)) ALLOCATE(US(nss)) ALLOCATE(doppler(nss)) ALLOCATE(window(nss)) ! 3.2 Doppler frequency doppler(1) = 0.0_wp DO i=2,nss/2+1 doppler(i) = (i-1.0_wp)/(nss*(time(2)-time(1))) doppler(nss-i+2) = -doppler(i) END DO doppler(:) = CSHIFT(doppler(:), nss/2) ! 3.3 Window function ! Yong Chen, raise cos function as filter function [0, 1, 0] DO i=1,nss window(i) = (1.0 + COS(2.0_wp*Pi*REAL(i-1-nss/2,wp)/REAL(nss,wp)))/2.0_wp END DO ! 3.4 Compute spectrum DO ic=1,nc IF (l_out) THEN lun = get_io_unit() IF (PRESENT(filnam)) THEN IF (TRIM(filnam) == "ropp_pp_spectra") THEN ! Default in ropp_pp_spectra_tool IF (ic == 1) text_filnam = "ROanalysis_dt_L1.dat" IF (ic == 2) text_filnam = "ROanalysis_dt_L2.dat" ELSE IF (ic == 1) text_filnam = TRIM(filnam)//"_dt_L1.dat" IF (ic == 2) text_filnam = TRIM(filnam)//"_dt_L2.dat" ENDIF ELSE ! Use original default output filenames IF (ic == 1) text_filnam = "ROanalysis_dt_L1.dat" IF (ic == 2) text_filnam = "ROanalysis_dt_L2.dat" ENDIF OPEN(UNIT=lun, FILE=TRIM(ADJUSTL(text_filnam))) WRITE(lun, '(A / A,I5,A,I5,A,A,I2)') & 'VARIABLES = "doppler (Hz)", "time(s)", "ln|U|", "HSL(m)", "SNR1", "SNR2"', & 'ZONE I=', nss, ' J=', ns, ' F=POINT', ' OCD=', ocd ENDIF ! 3.4.1 Compute wave number IF (ic == 1) k = 2.0_wp * pi * f_L1 / C_Light IF (ic == 2) k = 2.0_wp * pi * f_L2 / C_Light ! 3.4.2 Compute phase deviation DO i=1,n HSL(i) = impact_parameter(r_leo(i,:) - r_coc(:), r_gns(i,:) - r_coc(:)) - roc IF (ic == 1) dphase(i) = k*(phase_L1(i) - phase_LM(i)) IF (ic == 2) dphase(i) = k*(phase_L2(i) - phase_LM(i)) END DO ! 3.4.3 Complex field ! exp(i*k*(s-sm)) U(:) = EXP(Ci*(dphase(:))) DO is = is1, isN, di US(:) = U(is - di*nss/2:is + di*nss/2 - di:di)*window(:) CALL ropp_pp_FFT(US, 1) US(:) = US(:)/SQRT(REAL(nss,wp)) ! Normalization US(:) = US(:)/MAXVAL(ABS(US(:))) US(:) = CSHIFT(US(:), nss/2) IF (l_out) THEN DO i=1,nss ! WRITE(lun,'(E13.5,1X,F10.4,1X,E13.5)') & WRITE(lun, '(E20.10,1X,E20.10,1X,E20.10,1X,E20.10,1X,E20.10,1X,E20.10)') & doppler(i), time(is), LOG(ABS(US(i))), HSL(is), snr(1, is), snr(2, is) END DO ENDIF END DO ! 3.4.4 Output (thinned spectra) to netCDF files IF (l_out) THEN ! Simplest to read 2D fields from the test files just written REWIND (UNIT=lun) READ (lun, '(a80)') cline READ (lun, '(a80)') cline id_nx = INDEX( cline, '=') READ ( cline(id_nx+1:id_nx+6), '(i5)' ) nx id_ny = INDEX( cline(id_nx+1:), '=') READ ( cline(id_nx+id_ny+1:id_nx+id_ny+6), '(i5)' ) ny ALLOCATE (x(nx, ny), y(nx, ny), z(nx, ny)) DO iy=1,ny DO ix=1,nx READ (lun, *) x(ix, iy), y(ix, iy), z(ix, iy) END DO END DO CLOSE (UNIT=lun) ! Store as (thinned) ROprof 'extra data' CALL ropp_io_init(ro_data%lev2c, 1) CALL ropp_io_addvar(ro_data, & name='freq', & long_name='Frequency offset', & units='Hz', & range=(/ -1000.0_wp, 1000.0_wp /), & data=x(1:nx, 1)) CALL ropp_io_addvar(ro_data, & name='stime', & long_name='Time since start of occultation', & units='s', & range=(/ -1000.0_wp, 1000.0_wp /), & data=y(1, 1:ny:nsamp)) CALL ropp_io_addvar(ro_data, & name='amp', & long_name='log of modulus of signal amplitude', & units='', & range=(/ -1000.0_wp, 1000.0_wp /), & data=z(1:nx, 1:ny:nsamp)) DEALLOCATE (x, y, z) ! Write out nc_filnam = TRIM(ADJUSTL(text_filnam(1:INDEX(text_filnam, '.', BACK=.TRUE.)))) // 'nc' CALL ropp_io_write(ro_data, nc_filnam, ranchk=.TRUE.) CALL ropp_io_free(ro_data) ENDIF END DO ! loop over channels DEALLOCATE(dphase) DEALLOCATE(HSL) DEALLOCATE(U) DEALLOCATE(US) DEALLOCATE(doppler) DEALLOCATE(window) END SUBROUTINE ropp_pp_spectra