SUBROUTINE ropp_fm_roprof2obs1dbangle(ro_data, y) ! 1.1 Declarations ! ---------------- USE typesizes, ONLY: wp => EightByteReal USE ropp_utils USE ropp_io USE ropp_io_types, ONLY: ROprof USE ropp_fm USE ropp_fm_types, ONLY: Obs1dBangle USE ropp_fm_copy, not_this => ropp_fm_roprof2obs1dbangle USE geodesy IMPLICIT NONE TYPE(ROprof), INTENT(in) :: ro_data ! RO data structure TYPE(Obs1dBangle), INTENT(inout) :: y ! Bending angle structure INTEGER :: i, n, ipn INTEGER, DIMENSION(8) :: DT8 CHARACTER(len = 256) :: routine CHARACTER(len = 10) :: err_val CHARACTER(LEN=6) :: sr_leo ! 1.2 Error handling ! ------------------ CALL message_get_routine(routine) CALL message_set_routine('ropp_fm_roprof2obs (1D bending angles)') y%obs_ok = .TRUE. ! 1.3 Check and copy geolocation and time ! --------------------------------------- IF (isinrange(ro_data%georef%lon, ro_data%georef%range%lon)) THEN y%lon = ro_data%georef%lon ELSE WRITE(err_val, '(e8.1)') ro_data%georef%lon CALL message(msg_warn, & "Longitude data for observations out of range or missing. " // & "(longitude value = " // TRIM(err_val) // ")") CALL message(msg_warn, "Check input data and valid_range attributes") CALL message(msg_warn, "Set status flag obs%ok to FALSE") y%obs_ok = .FALSE. ENDIF IF (isinrange(ro_data%georef%lat, ro_data%georef%range%lat)) THEN y%lat = ro_data%georef%lat ELSE WRITE(err_val, '(e8.1)') ro_data%georef%lat CALL message(msg_warn, & "Latitude data for observations out of range or missing. " // & "(latitude value = " // TRIM(err_val) // ")") CALL message(msg_warn, "Check input data and valid_range attributes") CALL message(msg_warn, "Set status flag obs%ok to FALSE") y%obs_ok = .FALSE. ENDIF IF (isinrange(ro_data%georef%azimuth, ro_data%georef%range%azimuth)) THEN y%azimuth = ro_data%georef%azimuth ELSE WRITE(err_val, '(e8.1)') ro_data%georef%azimuth CALL message(msg_warn, "Azimuth (of " // TRIM(err_val) // " deg) out of range or missing.") CALL message(msg_warn, "Check input data and valid_range attributes.") CALL message(msg_warn, "Continuing, as azimuth not needed for 1D bending angles.") ENDIF IF ( isinrange(ro_data%DTocc%year, ro_data%DTocc%range%year) .AND. & isinrange(ro_data%DTocc%month, ro_data%DTocc%range%month) .AND. & isinrange(ro_data%DTocc%day, ro_data%DTocc%range%day) .AND. & isinrange(ro_data%DTocc%hour, ro_data%DTocc%range%hour) .AND. & isinrange(ro_data%DTocc%minute, ro_data%DTocc%range%minute) .AND. & isinrange(ro_data%DTocc%second, ro_data%DTocc%range%second) ) THEN DT8 = (/ro_data%DTocc%year, ro_data%DTocc%month, & ro_data%DTocc%day, 0, & ro_data%DTocc%hour, ro_data%DTocc%minute, & ro_data%DTocc%second, ro_data%DTocc%msec/) CALL TimeSince ( DT8, y%time, 1, Base="JS2000" ) ELSE CALL message(msg_warn, & "Time data for observations out of range or missing.") CALL message(msg_warn, "Check input data and valid_range attributes") CALL message(msg_warn, "Set status flag obs%ok to FALSE") y%obs_ok = .FALSE. ENDIF ! 1.4 Check that profiles are increasing in height - 1st element towards surface ! ------------------------------------------------------------------------------ CALL ropp_io_ascend(ro_data) ! 1.5 Check and copy observation data ! ----------------------------------- IF (ro_data%Lev1b%Npoints == 0) THEN CALL message(msg_warn, & "RO data has no level 1b (bending angle) parameters.") CALL message(msg_warn, "Check input data and valid_range attributes") CALL message(msg_warn, "Set status flag obs%ok to FALSE") y%obs_ok = .FALSE. RETURN ENDIF n = ro_data%Lev1b%Npoints IF (.NOT. ro_data%Lev2c%direct_ion) THEN ALLOCATE(y%impact(n), y%bangle(n), y%weights(n)) y%impact = ro_data%Lev1b%impact(1:n) y%bangle = ro_data%Lev1b%bangle(1:n) y%weights = 1.0_wp y%n_L1 = n ELSE ALLOCATE(y%impact(2*n), y%bangle(2*n), y%weights(2*n)) IF (ALL(ro_data%Lev1b%impact_L1 < ropp_MDTV)) THEN ! Use neutral IP if no L1 IPs y%impact(1:n) = ro_data%Lev1b%impact(1:n) ELSE y%impact(1:n) = ro_data%Lev1b%impact_L1(1:n) END IF IF (ALL(ro_data%Lev1b%impact_L2 < ropp_MDTV)) THEN ! Use neutral IP if no L2 IPs y%impact(n+1:2*n) = ro_data%Lev1b%impact(1:n) ELSE y%impact(n+1:2*n) = ro_data%Lev1b%impact_L2(1:n) END IF y%bangle(1:n) = ro_data%Lev1b%bangle_L1(1:n) y%bangle(n+1:2*n) = ro_data%Lev1b%bangle_L2(1:n) y%weights = 1.0_wp y%n_L1 = n END IF y%g_sfc = gravity(ro_data%georef%lat, 0.0_wp) y%r_earth = R_eff(ro_data%georef%lat) IF (isinrange(ro_data%georef%roc, ro_data%georef%range%roc)) THEN y%r_curve = ro_data%georef%roc ELSE CALL message(msg_warn, "Radius of curvature out of range or missing.") CALL message(msg_warn, "Check input data and valid_range attributes") CALL message(msg_warn, "Set status flag obs%ok to FALSE") y%obs_ok = .FALSE. ENDIF IF ( ro_data%Lev2c%direct_ion ) THEN ! Need the radius of the LEO orbit for single frequency work IF ( ANY ( ABS(ro_data%georef%leo_pod%pos - ropp_MDFV) < 1.0_wp ) .OR. & ANY ( ABS(ro_data%georef%r_coc - ropp_MDFV) < 1.0_wp ) ) THEN CALL message (msg_warn, 'Missing LEO position or centre of curvature ... ' // & 'will try to estimate r_LEO from satellite ID.') SELECT CASE (ro_data%leo_id(1:2)) ! Data from https://www.wmo-sat.info/oscar/ CASE ('C0') CALL message (msg_info, 'Assuming LEO height of 800 km (COSMIC-1 profile)') y%r_leo = y%r_curve + 800.E3_wp CASE ('C2') CALL message (msg_info, 'Assuming LEO height of 520 km (COSMIC-2a profile)') y%r_leo = y%r_curve + 520.E3_wp CASE ('CH') CALL message (msg_info, 'Assuming LEO height of 470 km (CHAMP profile)') y%r_leo = y%r_curve + 470.E3_wp CASE ('FY') CALL message (msg_info, 'Assuming LEO height of 840 km (FY-3C/D/... profile)') y%r_leo = y%r_curve + 840.E3_wp CASE ('GR') SELECT CASE (ro_data%leo_id(1:4)) CASE ('GRAA', 'GRAB') CALL message (msg_info, 'Assuming LEO height of 340 km (GRACE-A/B profile)') y%r_leo = y%r_curve + 340.E3_wp CASE ('GRAC', 'GRAD') CALL message (msg_info, 'Assuming LEO height of 490 km (GRACE-FO profile)') y%r_leo = y%r_curve + 490.E3_wp CASE DEFAULT CALL message(msg_warn, 'Unrecognised LEO satellite ID (' // ro_data%leo_id // & ') ... r_LEO cannot be estimated.') CALL message(msg_warn, 'Check that LEO POD and CoC positions, or LEO ID, are valid.') CALL message(msg_warn, 'Setting status flag obs%ok to FALSE') y%r_leo = ropp_MDFV y%obs_ok = .FALSE. END SELECT CASE ('KO') CALL message (msg_info, 'Assuming LEO height of 550 km (KOMPSAT-5 profile)') y%r_leo = y%r_curve + 550.E3_wp CASE ('ME') CALL message (msg_info, 'Assuming LEO height of 830 km (Metop-A/B/C profile)') y%r_leo = y%r_curve + 830.E3_wp CASE ('MG') CALL message (msg_info, 'Assuming LEO height of 870 km (Megha-Tropiques profile)') y%r_leo = y%r_curve + 870.E3_wp CASE ('PA') CALL message (msg_info, 'Assuming LEO height of 510 km (PAZ profile)') y%r_leo = y%r_curve + 510.E3_wp CASE DEFAULT CALL message(msg_warn, 'Unrecognised LEO satellite ID (' // ro_data%leo_id // & ') ... r_LEO cannot be estimated.') CALL message(msg_warn, 'Check that LEO POD and CoC positions, or LEO ID, are valid.') CALL message(msg_warn, 'Setting status flag obs%ok to FALSE') y%r_leo = ropp_MDFV y%obs_ok = .FALSE. END SELECT ELSE y%r_leo = SQRT ( SUM ( (ro_data%georef%leo_pod%pos - ro_data%georef%r_coc)**2 ) ) WRITE (sr_leo, FMT='(F6.1)') (y%r_leo - y%r_curve) * 1.0E-3_wp CALL message (msg_info, 'Using a LEO altitude of ' // sr_leo // ' km.') END IF ELSE y%r_leo = ropp_MDFV ! For safety's sake END IF ! Single frequency work IF (isinrange(ro_data%georef%undulation, ro_data%georef%range%undulation)) THEN y%undulation = ro_data%georef%undulation ELSE y%undulation = ropp_MDFV CALL message(msg_warn, "Undulation out of range or missing.") CALL message(msg_warn, "Check input data and valid_range attributes") CALL message(msg_warn, "Set status flag obs%ok to FALSE") y%obs_ok = .FALSE. END IF ! 1.6 Check and copy sigmas to diagonal error covariance ! ------------------------------------------------------ y%cov_ok = .TRUE. IF (ASSOCIATED(y%cov%d)) DEALLOCATE(y%cov%d) IF (.NOT. ro_data%Lev2c%direct_ion) THEN CALL callocate(y%cov%d, n*(n+1)/2) DO i = 1, n IF (ro_data%Lev1b%bangle(i) > 0.0_wp) THEN IF (ro_data%Lev1b%bangle_sigma(i) > 0.0_wp) THEN ! matrix_pp type, uplo = 'U' y%cov%d(i + i*(i-1)/2) = ro_data%Lev1b%bangle_sigma(i)**2 ELSE y%cov_ok = .FALSE. END IF ELSE y%cov%d(i + i*(i-1)/2) = 0.0003_wp ! = 1 deg**2 ENDIF END DO ELSE CALL callocate(y%cov%d, n*(2*n+1)) DO i = 1, n IF (ro_data%Lev1b%bangle_L1(i) > 0.0_wp) THEN IF (ro_data%Lev1b%bangle_L1_sigma(i) > 0.0_wp) THEN ! matrix_pp type, uplo = 'U' y%cov%d(i + i*(i-1)/2) = ro_data%Lev1b%bangle_L1_sigma(i)**2 ELSE y%cov_ok = .FALSE. END IF ELSE y%cov%d(i + i*(i-1)/2) = 0.0003_wp ! = 1 deg**2 END IF ipn = i + n IF (ro_data%Lev1b%bangle_L2(i) > 0.0_wp) THEN IF (ro_data%Lev1b%bangle_L2_sigma(i) > 0.0_wp) THEN ! matrix_pp type, uplo = 'U' y%cov%d(ipn + ipn*(ipn-1)/2) = ro_data%Lev1b%bangle_L2_sigma(i)**2 ELSE y%cov_ok = .FALSE. END IF ELSE y%cov%d(ipn + ipn*(ipn-1)/2) = 0.0003_wp ! = 1 deg**2 END IF END DO END IF IF (ASSOCIATED(y%cov%e)) DEALLOCATE(y%cov%e) IF (ASSOCIATED(y%cov%f)) DEALLOCATE(y%cov%f) IF (ASSOCIATED(y%cov%s)) DEALLOCATE(y%cov%s) y%cov%fact_chol = .FALSE. y%cov%equi_chol = 'N' ! 1.7 Clean up ! ------------ CALL message_set_routine(routine) END SUBROUTINE ropp_fm_roprof2obs1dbangle