! $Id: ropp_pp_fsi.f90 2022 2022-07-05 19:27:28Z ychen $ !****s* bangle/ropp_pp_fsi * ! ! NAME ! ropp_pp_fsi - Calculate L1 and L2 bending angle profiles using ! full spectrum inversion (FSI). ! ! SYNOPSIS ! call ropp_pp_fsi(time, snr, phase, r_leo, r_gns, r_coc, roc, w_ls, ! w_smooth, hmax, filter, opt_DL2, cff, dsh, ! impact, bangle, ba_cov, diag) ! ! DESCRIPTION ! This routine calculates L1 and L2 bending angles using a FSI algorithm. ! ! INPUTS ! real(wp), dimension(:) :: time ! Relative time of samples (s) ! real(wp), dimension(:,:) :: snr ! L1,L2 amplitudes [channel, time] ! real(wp), dimension(:,:) :: phase ! L1,L2 excess phase (m) [ch, time] ! real(wp), dimension(:,:) :: r_leo ! LEO coordinates (m) (ECI or ECF) ! real(wp), dimension(:,:) :: r_gns ! GPS coordinates (m) (ECI or ECF) ! real(wp), dimension(:) :: r_coc ! Centre curvature coordinates (m) ! real(wp) :: roc ! Radius of curvature (m) ! integer :: month ! observation month ! real(wp) :: hmax ! Maximum height for WO processing ! character(len=*) :: filter ! Filter method ('optest','slpoly') ! real(wp), dimension(:,:) :: impact ! L1,L2 impact parameters (m) [ch,t] ! real(wp), dimension(:,:) :: bangle ! L1,L2 bending angles (rad) [ch,t] ! ! OUTPUT ! real(wp), dimension(:,:) :: impact ! L1,L2 impact parameters (m) [ch, t] ! real(wp), dimension(:,:) :: bangle ! L1,L2 bending angles (rad) [ch, ip] ! real(wp), dimension(:,:) :: ba_cov ! Estimate of bangle covariance ! type(PPdiag), optional :: diag ! Additional output diagnostics structure ! ! 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 ! For further details please refer to the file COPYRIGHT ! which you should have received as part of this distribution. ! !**** SUBROUTINE ropp_pp_fsi(time, snr, phase, r_leo, r_gns, r_coc, roc, month, & hmax, opt_DL2, & impact, bangle, ba_cov, diag) !------------------------------------------------------------------------------- ! 1. Declarations !------------------------------------------------------------------------------- USE typesizes, ONLY: wp => EightByteReal USE ropp_pp, not_this => ropp_pp_fsi USE ropp_pp_constants, ONLY: c_light, pi, f_L1, f_L2 USE ropp_pp_spline USE ropp_pp_utils USE ropp_pp_types USE messages USE load_ciraq_model, only: read_ciraq_model USE star_fsi IMPLICIT NONE ! 1.1 Subroutine arguments REAL(wp), DIMENSION(:), INTENT(in) :: time ! Time of samples (s) REAL(wp), DIMENSION(:,:), INTENT(in) :: snr ! Amplitudes [ch, time] REAL(wp), DIMENSION(:,:), INTENT(in) :: phase ! Excess phase (m) [ch, t] REAL(wp), DIMENSION(:,:), INTENT(in) :: r_leo ! LEO coordinates (m) REAL(wp), DIMENSION(:,:), INTENT(in) :: r_gns ! GPS coordinates (m) REAL(wp), DIMENSION(:), INTENT(in) :: r_coc ! Centre of curvature (m) REAL(wp), INTENT(in) :: roc ! Radius of curvature (m) INTEGER, INTENT(in) :: month ! Observation momth REAL(wp), INTENT(in) :: hmax ! Maximum height for FSI LOGICAL, INTENT(in) :: opt_DL2 ! Degraded L2 REAL(wp), DIMENSION(:,:), INTENT(inout) :: impact ! Impact parameters (m) REAL(wp), DIMENSION(:,:), INTENT(inout) :: bangle ! Bending angles (rad) REAL(wp), DIMENSION(:,:), INTENT(inout) :: ba_cov ! Bending angle covariance TYPE(PPdiag), OPTIONAL, INTENT(inout) :: diag ! Additional diagnostics ! 1.2 Local parameters REAL(wp), PARAMETER :: PL = 1700.0_wp ! Lower estimate minimum ray height (meter) LOGICAL, PARAMETER :: opt_SH = .TRUE. ! Automatic shadow determination INTEGER, PARAMETER :: nmaxl1 = 65537 ! number of data point in the extended signal INTEGER, PARAMETER :: np_smooth = 3 ! Polynomial degree for smoothing regression REAL(wp), PARAMETER :: lowHSL = -50000.0_wp ! Low HSL for FSI smoothing (meter) REAL(wp), PARAMETER :: hws = 10000.0_wp ! impact height for different smoothing window size (meter) ! 1.3 Local scalars CHARACTER(len = 256) :: routine CHARACTER(len = 256) :: outstr CHARACTER(len = 256) :: benmod_file REAL(wp) :: Pmin, Pmax, Sec_LowSW ! min and max impact height, seconds using for smoothing below 10 km REAL(wp) :: dt, dx, dp, dE, df, pwmax, pwmin, dx_smooth, dpH INTEGER :: ocd, n, nc, ndf, ndf2, nx, i, ic INTEGER :: icw, imin, imax, ndE, nws, nr, nmaxlo, nmaxlsw REAL(wp), DIMENSION(:), ALLOCATABLE :: k ! Wave number [channel] REAL(wp), DIMENSION(:, :), ALLOCATABLE :: dphdt ! doppler shift [ch, t] REAL(wp), DIMENSION(:, :), ALLOCATABLE :: dphdt_sm ! smoothed doppler shift [ch, t] REAL(wp), DIMENSION(:, :), ALLOCATABLE :: amp_sm ! smoothed amplitude [ch, t] REAL(wp), DIMENSION(:, :), ALLOCATABLE :: dph_sm ! reconstructed excess phase [ch, t] REAL(wp), DIMENSION(:, :), ALLOCATABLE :: full_phase ! full phase after adjustment to circular orbits (m) [ch, t] REAL(wp), DIMENSION(:), ALLOCATABLE :: teta ! central angle between 1st and 2nd radius-vectors REAL(wp) :: r10, r20 ! absolute radius for GNSS (m) for LEO (m) INTEGER :: iflag ! 0~teta increasing with time; 1~teta decreasing LOGICAL :: apply_SH ! flag to apply shadow height determination INTEGER, DIMENSION(:), ALLOCATABLE :: nout ! array size for impact height and bending angle [channel] REAL(wp), DIMENSION(:, :), ALLOCATABLE :: pp ! impact parameter [ch, t] REAL(wp), DIMENSION(:, :), ALLOCATABLE :: ee ! bending angle [ch, t] REAL(wp), DIMENSION(:, :), ALLOCATABLE :: aa ! FSI amplitude [ch, t] INTEGER, DIMENSION(:), ALLOCATABLE :: noutL ! array size for impact height and bending angle [channel] REAL(wp), DIMENSION(:, :), ALLOCATABLE :: ppL ! impact parameter [ch, t] REAL(wp), DIMENSION(:, :), ALLOCATABLE :: eeL ! bending angle [ch, t] REAL(wp), DIMENSION(:, :), ALLOCATABLE :: aaL ! FSI amplitude [ch, t] INTEGER, DIMENSION(:), ALLOCATABLE :: noutf ! array size for impact height and bending angle [channel] REAL(wp), DIMENSION(:, :), ALLOCATABLE :: ppf ! impact parameter [ch, t] REAL(wp), DIMENSION(:, :), ALLOCATABLE :: eef ! bending angle [ch, t] REAL(wp), DIMENSION(:, :), ALLOCATABLE :: aaf ! FSI amplitude [ch, t] REAL(wp), DIMENSION(:), ALLOCATABLE :: spike_th ! spike threshold REAL(wp), DIMENSION(:), ALLOCATABLE :: pp_out ! impact parameter [t] REAL(wp), DIMENSION(:), ALLOCATABLE :: ee_out ! bending angle [t] REAL(wp), DIMENSION(:), ALLOCATABLE :: aa_out ! FSI amplitude [t] INTEGER, DIMENSION(:), ALLOCATABLE :: idx ! array index for the keeped impact parameter and bending angles INTEGER :: n_pe ! aaray size REAL(wp) :: spike_thresh ! spike threshold REAL(wp), DIMENSION(:), ALLOCATABLE :: ppL1 ! impact parameter [t] REAL(wp), DIMENSION(:), ALLOCATABLE :: eeL1 ! bending angle [t] REAL(wp), DIMENSION(:), ALLOCATABLE :: aaL1 ! FSI amplitude [t] INTEGER, DIMENSION(:), ALLOCATABLE :: idxL1 ! array index for the keeped impact parameter and bending angles INTEGER :: npts ! aaray size REAL(wp), DIMENSION(:, :), ALLOCATABLE :: ppLSW ! LSW mean impact height REAL(wp), DIMENSION(:, :), ALLOCATABLE :: eeLSW ! LSW mean bending angle REAL(wp), DIMENSION(:, :), ALLOCATABLE :: LSWMG ! bending angle LSW covariance using squared dalp [ch, t] REAL(wp), DIMENSION(:, :), ALLOCATABLE :: LSWSS ! bending angle LSW covariance using absolute dalp [ch, t] INTEGER, DIMENSION(:), ALLOCATABLE :: noutlws ! array size for LSW impact height and bending REAL(wp), DIMENSION(:), ALLOCATABLE :: PminC ! Estimate of shadow border REAL(wp), DIMENSION(:), ALLOCATABLE :: EW ! Bangle for re-interpolation REAL(wp), DIMENSION(:), ALLOCATABLE :: ECW ! Ba cov for re-interpolation LOGICAL, DIMENSION(:), ALLOCATABLE :: MW ! Mask for computation of dE REAL(wp), DIMENSION(:), ALLOCATABLE :: smooth_ba ! 1.4 Local arrays CALL message_get_routine(routine) CALL message_set_routine('ropp_pp_fsi') !------------------------------------------------------------------------------- ! 2. Initialization !------------------------------------------------------------------------------- ! 2.1 Determination of data sizes n = SIZE(time) nc = SIZE(snr,1) nr = CEILING(REAL(n,wp)/8000.0_wp) ALLOCATE(k(nc)) ! Computation of wave vectors k(1) = 2*pi*f_L1/C_Light k(2) = 2*pi*f_L2/C_Light ! Use CIRA+Q bending angle model to propogate to circular orbits benmod_file='BENMOD' call read_ciraq_model(benmod_file) ! rising/setting flag (rising 1, setting -1) ocd = NINT(SIGN(1.0_wp, impact(1, n) - impact(1, 1))) ! Determine the sampling frequency (normally 50 Hz, for COMSIC-2, 100 Hz) dt = abs(time(n) - time(1))/n ! For impact height above 10 km, using sampling rate as smoothing window ndf = NINT(1.0_wp/dt) ! Yong Chen 11/08/2023 to increasing the smoothing window size below 10 km for RO data rate of 50 Hz ! for RO data rate of 50 Hz IF (ndf < 55 ) THEN Sec_LowSW = 0.2_wp ! for RO data rate of 100 Hz ELSE Sec_LowSW = 0.1_wp ENDIF ndf2 = NINT(1.0_wp*ndf/2) ! lower HSL (-50 km) position for raising/setting profile if (ocd == 1 ) THEN imax = MIN(SUM(MINLOC(impact(1,:), impact(1,:) >= roc + lowHSL)), & SUM(MAXLOC(impact(1,:), impact(1,:) >= roc + lowHSL))) ELSE imax = MAX(SUM(MINLOC(impact(1,:), impact(1,:) >= roc + lowHSL)), & SUM(MAXLOC(impact(1,:), impact(1,:) >= roc + lowHSL))) ENDIF ! array size for raising/setting profile above HSL -50 km IF (ocd == 1) THEN nmaxlo = n - imax + 1 ELSE nmaxlo = imax ENDIF write(*, *) 'n, imax, nmaxlo=', ocd, n, imax, nmaxlo, impact(1,1), impact(1,n), roc + lowHSL ! Determine the Doppler ALLOCATE(dphdt(nc, nmaxlo)) ALLOCATE(dphdt_sm(nc,nmaxlo)) ALLOCATE(amp_sm(nc,nmaxlo)) ALLOCATE(dph_sm(nc,nmaxlo)) ALLOCATE(full_phase(nc,nmaxlo)) ALLOCATE(teta(nmaxlo)) ! Apply Fourier smoothing of the excess phase ! Reconstruct phase from smoothed Doppler ! Rising profile IF (ocd == 1) THEN CALL phs2dop(phase(:,n:n-nmaxlo+1:-1), dphdt(:,:)) DO ic=1,nc CALL fsi_fourier_filter(nmaxlo, ndf2, dphdt(ic, :), dphdt_sm(ic, :)) CALL fsi_fourier_filter(nmaxlo, ndf, snr(ic, n:n-nmaxlo+1:-1), amp_sm(ic, :)) ENDDO CALL dop2phs(phase(:,n:n-nmaxlo+1:-1), dphdt_sm, dph_sm) ELSE ! Setting profile CALL phs2dop(phase(:,1:nmaxlo), dphdt) DO ic=1,nc CALL fsi_fourier_filter(nmaxlo, ndf2, dphdt(ic, :), dphdt_sm(ic, :)) CALL fsi_fourier_filter(nmaxlo, ndf, snr(ic, 1:nmaxlo), amp_sm(ic, :)) ENDDO CALL dop2phs(phase(:,1:nmaxlo), dphdt_sm, dph_sm) ENDIF ! adjustment phase for circular orbits IF (ocd == 1) THEN CALL gpsleofix(nmaxlo, r_gns(n:n-nmaxlo+1:-1, :), r_leo(n:n-nmaxlo+1:-1, :), r_coc, & dph_sm, month, r10, r20, teta, full_phase, iflag) IF (iflag /= 0) THEN CALL message(msg_error,'Rising profile is not recognised. Check gpsleofix. ') ENDIF ELSE CALL gpsleofix(nmaxlo, r_gns(1:nmaxlo, :), r_leo(1:nmaxlo, :), r_coc, & dph_sm, month, r10, r20, teta, full_phase, iflag) IF (iflag /= 0) THEN CALL message(msg_error,'Setting profile is not recognised. Check gpsleofix. ') ENDIF ENDIF ! minimun and maximum impact height from L1 (10 km to hmax (25) km) Pmax = MIN(hmax, MAXVAL(impact(1,:)) - roc) Pmin = PL ! hws ! Full-Spectrum-Inversion (FSI) (using open angle as a basic variable instead of time) ! teta should increasing with time (already treated as setting profile) ALLOCATE(pp(nc,nmaxl1)) ALLOCATE(ee(nc,nmaxl1)) ALLOCATE(aa(nc,nmaxl1)) ALLOCATE(nout(nc)) !apply_SH = .FALSE. apply_SH = .TRUE. CALL fsi_leo(teta, full_phase, amp_sm, r10, r20, pp, ee, aa, Pmin+roc, Pmax+roc, apply_SH, nout) write(*,*)'min max impH after fsi_leo', pp(1, 1)-roc, pp(1, nout(1))-roc, & pp(2, 1)-roc, pp(2, nout(2))-roc DEALLOCATE(dphdt) DEALLOCATE(dphdt_sm) DEALLOCATE(amp_sm) DEALLOCATE(dph_sm) DEALLOCATE(full_phase) DEALLOCATE(teta) ! only apply to L1 signal to remove spikes ALLOCATE(ppL1(nmaxl1)) ALLOCATE(eeL1(nmaxl1)) ALLOCATE(aaL1(nmaxl1)) ALLOCATE(idxL1(nmaxl1)) dx = 250.0_wp dp = abs(pp(1, nout(1)) - pp(1, 1))/(nout(1) -1) nx = NINT(dx/dp) spike_thresh = 0.5_wp IF ( MOD(nx,2) /= 1) THEN nx = nx-1 ENDIF CALL remove_spikes(pp(1, 1:nout(1)), ee(1, 1:nout(1)), & nx, spike_thresh, ppL1, eeL1, idxL1, npts) write(*,*) 'npts', npts, nout(1) DO i=1, npts aaL1(i) = aa(1, idxL1(i)) ENDDO ! For lower level of data below 10 km impact height with smaller smoothing window size !ndf = NINT(0.1_wp/dt) ! Yong Chen 11/08/2023 to change smoothing window size below 10 km based on RO data rate of 50 Hz or 100 Hz ndf = NINT(Sec_LowSW/dt) ndf2 = NINT(1.0_wp * ndf/2) ! Determine the Doppler ALLOCATE(dphdt(nc, n)) ALLOCATE(dphdt_sm(nc,n)) ALLOCATE(amp_sm(nc,n)) ALLOCATE(dph_sm(nc,n)) ALLOCATE(full_phase(nc,n)) ALLOCATE(teta(n)) ! Apply Fourier smoothing of the excess phase ! Reconstruct phase from smoothed Doppler ! Rising profile IF (ocd == 1) THEN CALL phs2dop(phase(:,n:1:-1), dphdt(:,:)) DO ic=1,nc CALL fsi_fourier_filter(n, ndf2, dphdt(ic, :), dphdt_sm(ic, :)) CALL fsi_fourier_filter(n, ndf, snr(ic, n:1:-1), amp_sm(ic, :)) ENDDO CALL dop2phs(phase(:,n:1:-1), dphdt_sm, dph_sm) ELSE ! Setting profile CALL phs2dop(phase, dphdt) DO ic=1,nc CALL fsi_fourier_filter(n, ndf2, dphdt(ic, :), dphdt_sm(ic, :)) CALL fsi_fourier_filter(n, ndf, snr(ic, :), amp_sm(ic, :)) ENDDO CALL dop2phs(phase, dphdt_sm, dph_sm) ENDIF ! adjustment phase for circular orbits IF (ocd == 1) THEN CALL gpsleofix(n, r_gns(n:1:-1, :), r_leo(n:1:-1, :), r_coc, & dph_sm, month, r10, r20, teta, full_phase, iflag) IF (iflag /= 0) THEN CALL message(msg_error,'Rising profile is not recognised. Check gpsleofix. ') ENDIF ELSE CALL gpsleofix(n, r_gns, r_leo, r_coc, dph_sm, month, r10, r20, teta, full_phase, iflag) IF (iflag /= 0) THEN CALL message(msg_error,'Setting profile is not recognised. Check gpsleofix. ') ENDIF ENDIF ! maximum impact height from L1 Pmax = MIN(hmax, MAXVAL(impact(1,:)) - roc) Pmin = PL !MIN(PL, MINVAL(impact(1,:)) - roc) ! Full-Spectrum-Inversion (FSI) (using open angle as a basic variable instead of time) ! teta should increasing with time (already treated as setting profile) ALLOCATE(ppL(nc,nmaxl1)) ALLOCATE(eeL(nc,nmaxl1)) ALLOCATE(aaL(nc,nmaxl1)) ALLOCATE(noutL(nc)) apply_SH = .TRUE. CALL fsi_leo(teta, full_phase, amp_sm, r10, r20, ppL, eeL, aaL, Pmin+roc, Pmax+roc, apply_SH, noutL) write(*,*)'min max impH after fsi_leo', ppL(1, 1)-roc, ppL(1, noutL(1))-roc, & ppL(2, 1)-roc, ppL(2, noutL(2))-roc ! calculate the LSW covariance nmaxlsw = MAXVAL(noutL) ALLOCATE(ppLSW(nc, nmaxlsw)) ALLOCATE(eeLSW(nc, nmaxlsw)) ALLOCATE(LSWMG(nc, nmaxlsw)) ALLOCATE(LSWSS(nc, nmaxlsw)) ALLOCATE(noutlws(nc)) DO ic=1,nc CALL fsi_lsw(roc, ppL(ic, 1:noutL(ic)), eeL(ic, 1:noutL(ic)), & ppLSW(ic, :), eeLSW(ic, :), LSWMG(ic, :), LSWSS(ic, :), noutlws(ic)) ENDDO DEALLOCATE(dphdt) DEALLOCATE(dphdt_sm) DEALLOCATE(amp_sm) DEALLOCATE(dph_sm) DEALLOCATE(full_phase) DEALLOCATE(teta) ALLOCATE(ppf(nc,nmaxl1)) ALLOCATE(eef(nc,nmaxl1)) ALLOCATE(aaf(nc,nmaxl1)) ALLOCATE(noutf(nc)) ! combine the data from Pmin to Pmax, with seperation at 10 km for L1 ! upper level from 10 km to 25 km using large smoothing window ! lower level below 10 km using small smmothing window ! for L2 using the large smoothing window pp, ee, and aa DO ic=1,nc IF(ic == 1) THEN imin = SUM(MAXLOC(ppL1(1:npts), ppL1(1:npts) <= hws+roc)) IF(imin > 0) THEN write(*, *) '10 km impact height imin =', imin, ppL1(imin+1), hws+roc ppf(ic, 1:imin) = ppL1(1:imin) eef(ic, 1:imin) = eeL1(1:imin) aaf(ic, 1:imin) = aaL1(1:imin) noutf(ic) = imin ELSE ! all impact heights are higher than 10 km ppf(ic, 1:npts) = ppL1(1:npts) eef(ic, 1:npts) = eeL1(1:npts) aaf(ic, 1:npts) = aaL1(1:npts) noutf(ic) = npts ENDIF imax = SUM(MAXLOC(ppL(ic, 1:noutL(ic)), ppL(ic, 1:noutL(ic)) <= hws+roc)) IF ( imax > 0) THEN !imapct height are lower than 10 km write(*, *) '10 km impact height imax =', imin, ppL(ic, imax-1), hws+roc ppf(ic, noutf(ic)+1:noutf(ic)+noutL(ic)-imax) = ppL(ic, imax+1:noutL(ic)) eef(ic, noutf(ic)+1:noutf(ic)+noutL(ic)-imax) = eeL(ic, imax+1:noutL(ic)) aaf(ic, noutf(ic)+1:noutf(ic)+noutL(ic)-imax) = aaL(ic, imax+1:noutL(ic)) noutf(ic) = noutf(ic)+noutL(ic)-imax ENDIF ELSE ! L2 using the large smoothing window pp, ee, and aa ppf(ic, 1:nout(ic)) = pp(ic, 1:nout(ic)) eef(ic, 1:nout(ic)) = ee(ic, 1:nout(ic)) aaf(ic, 1:nout(ic)) = aa(ic, 1:nout(ic)) noutf(ic) = nout(ic) ENDIF ENDDO DEALLOCATE(pp) DEALLOCATE(ee) DEALLOCATE(aa) DEALLOCATE(nout) DEALLOCATE(ppL) DEALLOCATE(eeL) DEALLOCATE(aaL) DEALLOCATE(noutL) DEALLOCATE(ppL1) DEALLOCATE(eeL1) DEALLOCATE(aaL1) DEALLOCATE(idxL1) ALLOCATE(PminC(nc)) ALLOCATE(spike_th(nc)) PminC(1) = Pmin PminC(2) = Pmin ! 250 m segment length ! set threshold at +- 75% (100%) of mean spike_th(1) = 0.75_wp spike_th(2) = 1.0_wp dx = 250.0_wp dx_smooth = 250.0_wp Channels: DO ic=1,nc ALLOCATE(pp_out(nmaxl1)) ALLOCATE(ee_out(nmaxl1)) ALLOCATE(aa_out(nmaxl1)) ALLOCATE(idx(nmaxl1)) ! remove outliers before applying smoothing dp = abs(ppf(ic, noutf(ic)) - ppf(ic, 1))/(noutf(ic) -1) nx = NINT(dx/dp) IF ( MOD(nx,2) /= 1) THEN nx = nx-1 ENDIF CALL remove_spikes(ppf(ic, 1:noutf(ic)), eef(ic, 1:noutf(ic)), & nx, spike_th(ic), pp_out, ee_out, idx, n_pe) write(*,*) 'n_pe', n_pe write(*,*) 'ee_out', ee_out(n_pe-10:n_pe) ALLOCATE(smooth_ba(n_pe)) DO i=1, n_pe aa_out(i) = aaf(ic, idx(i)) ENDDO pwmax = MAXVAL(pp_out(1:n_pe)) pwmin = MINVAL(pp_out(1:n_pe)) dpH = nr* ABS(pwmax - pwmin)/(n_pe - 1) nws = NINT(REAL(dx_smooth,wp) / dpH) write(*,*) 'dph, nws, nr', dph, nws, nr CALL ropp_pp_sliding_polynomial(pp_out(1:n_pe:nr), ee_out(1:n_pe:nr), & nws, np_smooth, smooth_ba(1:n_pe:nr)) IF (nr > 1) THEN CALL ropp_pp_interpol(pp_out(1:n_pe:nr), pp_out(1:n_pe), & smooth_ba(1:n_pe:nr), smooth_ba(1:n_pe)) ENDIF !write(*,*) 'smooth_ba', smooth_ba(n_pe-10:n_pe) ! Re-define Pmin IF (opt_SH) THEN PminC(ic) = min(pp_out(1),pp_out(n_pe)) - roc ELSE PminC(ic) = Pmin ENDIF IF(opt_DL2) THEN Pmin = PminC(ic) ELSE Pmin = PminC(1) ENDIF WRITE(outstr, '(2(A,F10.3))') 'Pmin = ', Pmin, ' Pmax= ', Pmax CALL message(msg_diag, outstr) WRITE(outstr, '(2(A,F18.3))') 'pp_out(1) = ', pp_out(1), ' pp_out(n_pe)= ', pp_out(n_pe) CALL message(msg_diag, outstr) ! 11.5 Redefenition of impact parameter grid imax = SUM(MAXLOC(impact(ic,:), impact(ic,:) < roc + Pmax)) IF (imax > n .OR. imax < 1) THEN CALL message(msg_warn, 'No data for FSI processing - will not process') RETURN ENDIF IF (impact(ic,1) < roc + Pmax) THEN imin = 1 ELSE imin = n ENDIF DO i=imin,imax,SIGN(1,imax-imin) impact(ic,i) = roc + (Pmin*REAL(imax-i,wp) + Pmax*REAL(i-imin,wp)) / & REAL(imax-imin,wp) ENDDO ! 11.6 Interpolation of bending angle and covariances onto impact grid ba_cov(ic,:) = 0.0_wp DO i=1,n IF (impact(ic,i) - roc < Pmax) THEN CALL ropp_pp_interpol(pp_out(1:n_pe), impact(ic,i), & smooth_ba(1:n_pe), bangle(ic,i)) IF ( noutlws(ic) > 0 ) THEN CALL ropp_pp_interpol(ppLSW(ic, 1:noutlws(ic)), impact(ic,i), & LSWMG(ic, 1:noutlws(ic)), ba_cov(ic,i)) ENDIF ENDIF ENDDO ! Interpolation amplititude to standard grid or output amplitude to diag DEALLOCATE(pp_out) DEALLOCATE(ee_out) DEALLOCATE(aa_out) DEALLOCATE(idx) DEALLOCATE(smooth_ba) ENDDO Channels DEALLOCATE(ppLSW) DEALLOCATE(eeLSW) DEALLOCATE(LSWMG) DEALLOCATE(LSWSS) DEALLOCATE(noutlws) DEALLOCATE(ppf) DEALLOCATE(eef) DEALLOCATE(aaf) DEALLOCATE(noutf) !------------------------------------------------------------------------------- ! 13. Synchronizing grids for all channels !------------------------------------------------------------------------------- IF (opt_DL2) THEN ! 13.1 Determination of Pmin from best channel icw = 1 Pmin = PminC(icw) WRITE(outstr, '(2X,A,I2,A,F10.3)') 'ICW = ', ICW, ' Pmin = ', Pmin CALL message(msg_diag, outstr) ! 13.2 Grid redefinition and interpolation ALLOCATE(EW(n)) ALLOCATE(ECW(n)) ALLOCATE(MW(n)) DO ic=1,nc CALL ropp_pp_interpol(impact(ic,:), impact(icw,:), bangle(ic,:), EW(:), & Cext=.true.) CALL ropp_pp_interpol(impact(ic,:), impact(icw,:), ba_cov(ic,:), ECW(:), & Cext=.true.) MW(:) = (impact(icw,:) >= roc + PminC(ic) + 1000.0 .AND. & impact(icw,:) <= roc + PminC(ic) + 6000.0) ndE = COUNT(Mask = MW(:)) IF (ndE > 20) THEN dE = SUM(EW(:) - bangle(icw,:), Mask = MW(:))/ndE ELSE dE = 0.0_wp ENDIF WHERE (impact(icw,:) < roc+PminC(ic)) bangle(ic,:) = bangle(icw,:) + dE ba_cov(ic,:) = ba_cov(icw,:) ELSEWHERE bangle(ic,:) = EW(:) ba_cov(ic,:) = ECW(:) END WHERE impact(ic,:) = impact(icw,:) ENDDO DEALLOCATE(EW) DEALLOCATE(ECW) DEALLOCATE(MW) ENDIF !------------------------------------------------------------------------------- ! 14. Clean up !------------------------------------------------------------------------------- DEALLOCATE(PminC) DEALLOCATE(k) DEALLOCATE(spike_th) END SUBROUTINE ropp_pp_fsi