! $Id: ropp_pp_bending_angle_fsi.f90 2022 2022-07-05 10:49:04Z ychen $ SUBROUTINE ropp_pp_bending_angle_fsi(time, r_leo, r_gns, r_coc, roc, & phase_L1, phase_L2, snr_L1, snr_L2, month, & hmax, opt_DL2, & impact_L1, bangle_L1, ba_sigma_L1, & impact_L2, bangle_L2, ba_sigma_L2, diag) !****s* bangle/ropp_pp_bending_angle_fsi * ! ! NAME ! ropp_pp_bending_angle_fsi - Calculate L1 and L2 bending angle profiles from ! occultation data by Full Spectrum Inversion ! ! SYNOPSIS ! call ropp_pp_bending_angle_fsi(time, r_leo, r_gns, r_coc, phase_L1, ! phase_L2, snr_L1, snr_L2, month, ! hmax, opt_DL2, ! impact_L1, bangle_L1, ba_sigma_L1, ! impact_L2, bangle_L2, ba_sigma_L2, diag) ! ! DESCRIPTION ! This routine calculates L1 and L2 bending angles using the FSI algorithm. ! ! INPUTS ! real(wp), dimension(:) :: time ! Relative time of samples (s) ! 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 coords (m) ! real(wp) :: roc ! Radius curvature (m) ! integer :: month ! observation month ! integer :: w_ls ! Large-scale smoothing (points) ! integer :: w_smooth ! Smoothing window above 7km (point) ! integer :: w_low ! Smoothing window below 7km (point) ! real(wp) :: hmax ! Maximum height for WO (m) ! character(len=*) :: filter ! Filter method ('optest'/'slpoly') ! logical :: opt_DL2 ! Degraded L2 flag ! integer :: cff ! Complex filtering flag ! real(wp) :: dsh ! Shadow border width (m) ! real(wp), dimension(:) :: phase_L1 ! L1 excess phase (m) ! real(wp), dimension(:) :: phase_L2 ! L2 excess phase (m) ! real(wp), dimension(:) :: snr_L1 ! L1 amplitude ! real(wp), dimension(:) :: snr_L2 ! L2 amplitude ! real(wp), dimension(:) :: impact_L1 ! L1 impact parameters (m) ! real(wp), dimension(:) :: bangle_L1 ! L1 bending angles (rad) ! real(wp), dimension(:) :: impact_L2 ! L2 impact parameters (m) ! real(wp), dimension(:) :: bangle_L2 ! L2 bending angles (rad) ! ! OUTPUT ! real(wp), dimension(:) :: impact_L1 ! L1 impact parameters (m) ! real(wp), dimension(:) :: bangle_L1 ! L1 bending angles (rad) ! real(wp), dimension(:) :: ba_sigma_L1 ! L1 bending angle error (rad) ! real(wp), dimension(:) :: impact_L2 ! L2 impact parameters (m) ! real(wp), dimension(:) :: bangle_L2 ! L2 bending angles (rad) ! real(wp), dimension(:) :: ba_sigma_L2 ! L2 bending angle error (rad) ! type(PPdiag), optional :: diag ! Additional output diagnostic strt ! ! 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. ! !**** !------------------------------------------------------------------------------- ! 1. Declarations !------------------------------------------------------------------------------- USE typesizes, ONLY: wp => EightByteReal USE ropp_pp, not_this => ropp_pp_bending_angle_fsi USE ropp_pp_types, ONLY: PPDiag USE messages IMPLICIT NONE REAL(wp), DIMENSION(:), INTENT(in) :: time ! Time of samples (s) 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 curvature (m) REAL(wp), INTENT(in) :: roc ! Radius of curvature (m) REAL(wp), DIMENSION(:), INTENT(in) :: phase_L1 ! L1 excess phase (m) REAL(wp), DIMENSION(:), INTENT(in) :: phase_L2 ! L2 excess phase (m) REAL(wp), DIMENSION(:), INTENT(in) :: snr_L1 ! L1 amplitude REAL(wp), DIMENSION(:), INTENT(in) :: snr_L2 ! L2 amplitude INTEGER, INTENT(in) :: month ! Observation momth REAL(wp), INTENT(in) :: hmax ! Max height for WO (m) LOGICAL, INTENT(in) :: opt_DL2 ! Degraded L2 flag REAL(wp), DIMENSION(:), INTENT(inout) :: impact_L1 ! L1 impact parameter (m) REAL(wp), DIMENSION(:), INTENT(inout) :: bangle_L1 ! L1 bending angles (rad) REAL(wp), DIMENSION(:), INTENT(out) :: ba_sigma_L1 ! L1 bangle std dev REAL(wp), DIMENSION(:), INTENT(inout) :: impact_L2 ! L2 impact parameter (m) REAL(wp), DIMENSION(:), INTENT(inout) :: bangle_L2 ! L2 bending angles (rad) REAL(wp), DIMENSION(:), INTENT(out) :: ba_sigma_L2 ! L2 bangle std dev TYPE(PPdiag), OPTIONAL, INTENT(inout) :: diag ! Additional diagnostics REAL(wp), DIMENSION(:,:), ALLOCATABLE :: AP REAL(wp), DIMENSION(:), ALLOCATABLE :: A0 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: snr ! Amplitude array [ch,t] REAL(wp), DIMENSION(:,:), ALLOCATABLE :: phase ! Phase array [ch,t] REAL(wp), DIMENSION(:,:), ALLOCATABLE :: impact ! Impact array [ch, t] REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bangle ! Bangle array [ch, t] REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ba_cov ! Bangle covariance INTEGER :: ocd ! Occultation direction INTEGER :: n ! Number of data points CHARACTER(len = 256) :: routine CALL message_get_routine(routine) CALL message_set_routine('ropp_pp_bending_angle_fsi') !------------------------------------------------------------------------------- ! 2. Ensure monotonous impact parameter grid !------------------------------------------------------------------------------- n = SIZE(impact_L1) ocd = NINT(SIGN(1.0_wp, impact_L1(n)-impact_L1(1))) CALL ropp_pp_monotonous(impact_L1, -1) CALL ropp_pp_monotonous(impact_L2, -1) !------------------------------------------------------------------------------- ! 3. Create combined data arrays [channel, time] !------------------------------------------------------------------------------- ALLOCATE(snr(2,n)) ALLOCATE(phase(2,n)) ALLOCATE(impact(2,n)) ALLOCATE(bangle(2,n)) ALLOCATE(ba_cov(2,n)) ALLOCATE(AP(2,n)) ALLOCATE(A0(2)) snr(1,:) = snr_L1(:) snr(2,:) = snr_L2(:) phase(1,:) = phase_L1(:) phase(2,:) = phase_L2(:) impact(1,:) = impact_L1(:) impact(2,:) = impact_L2(:) bangle(1,:) = bangle_L1(:) bangle(2,:) = bangle_L2(:) !------------------------------------------------------------------------------- ! 4. Scale lowest amplitude data point !------------------------------------------------------------------------------- IF (ocd == 1) THEN snr(:,1) = 1e-5_wp*MAXVAL(snr(:,:)) ELSE snr(:,n) = 1e-5_wp*MAXVAL(snr(:,:)) END IF !------------------------------------------------------------------------------- ! 5. Perform full spectrum inversion (fsi) !------------------------------------------------------------------------------- IF (PRESENT(diag)) THEN CALL ropp_pp_fsi(time, snr, phase, r_leo, r_gns, r_coc, roc, month, & hmax, opt_DL2, impact, bangle, ba_cov, diag) ELSE CALL ropp_pp_fsi(time, snr, phase, r_leo, r_gns, r_coc, roc, month, & hmax, opt_DL2, impact, bangle, ba_cov) ENDIF !------------------------------------------------------------------------------- ! 6. Update output variables with computed bending angle and impact parameter !------------------------------------------------------------------------------- impact_L1(:) = impact(1,:) impact_L2(:) = impact(2,:) bangle_L1(:) = bangle(1,:) bangle_L2(:) = bangle(2,:) ba_sigma_L1(:) = SQRT(ba_cov(1,:)) ba_sigma_L2(:) = SQRT(ba_cov(2,:)) !------------------------------------------------------------------------------- ! 7. Clean up !------------------------------------------------------------------------------- DEALLOCATE(ba_cov) DEALLOCATE(bangle) DEALLOCATE(impact) DEALLOCATE(phase) DEALLOCATE(snr) DEALLOCATE(AP) DEALLOCATE(A0) CALL message_set_routine(routine) END SUBROUTINE ropp_pp_bending_angle_fsi