prepare_flux_data Subroutine

   subroutine prepare_flux_data(sw_data, flux_data, utc_mjd, kp_out, f107_out, f107a_out)
      type(sw_data_type), intent(inout) :: sw_data !! Space weather data manager
      type(flux_data_type), intent(in) :: flux_data !! Current day's flux data
      real(dp), intent(in) :: utc_mjd !! Current epoch (MJD)
      real(dp), intent(out) :: kp_out !! Selected Kp index for current epoch
      real(dp), intent(out) :: f107_out !! Selected F10.7 daily value for current epoch
      real(dp), intent(out) :: f107a_out !! Selected F10.7a average for current epoch

      real(dp) :: frac_epoch !! Fractional day from midnight of flux_data%mjd
      real(dp) :: frac_epoch_kp !! Fractional day for Kp selection (accounts for 6.7 hour lag)
      real(dp) :: f107_offset !! Time of day offset for F10.7 validity boundary
      integer(ip) :: sub_index !! Index for 3-hour Kp period (0-7 for 8 periods per day)
      integer(ip) :: f107_index !! Index for F10.7 selection (0 for current day, -1 for previous day)
      logical :: sw_status !! Status flag for SW data retrieval
      type(flux_data_type) :: flux_data_prev !! Previous day's flux data (for Kp and F10.7 fallback)
      type(flux_data_type) :: flux_data_2days_ago !! Flux data from 2 days ago (for F10.7 fallback)

      real(dp), parameter :: F107_REF_EPOCH = 48407.5_dp  !! MJD for 5/31/91 noon (matches C++ GSFC MJD 18408.0)

      logical,parameter :: use_gmat_bug = .false. !! Set to .false. to use corrected logic
                                                  !! for F10.7a selection instead of GMAT's
                                                  !! apparent bug

      ! Calculate fractional day from midnight
      frac_epoch = utc_mjd - flux_data%mjd

      ! Apply 6.7 hour lag for Kp per Vallado & Finkleman
      frac_epoch_kp = frac_epoch - 6.7_dp / 24.0_dp

      ! Determine which 3-hour period (0-7 for 8 periods per day)
      sub_index = min(floor(frac_epoch_kp * 8.0_dp), 7)

      ! F10.7 is measured at 8pm (5pm before 5/31/91), but the data row is
      ! valid from 8am on the current day to 8am the next day (5am before ref epoch)
      ! So f107_offset represents the DATA VALIDITY BOUNDARY, not measurement time
      if (utc_mjd < F107_REF_EPOCH) then
         f107_offset = 5.0_dp / 24.0_dp  ! 5am validity boundary
      else
         f107_offset = 8.0_dp / 24.0_dp  ! 8am validity boundary
      end if

      ! Determine f107_index based on time of day
      if (frac_epoch < f107_offset) then
         f107_index = -1  ! Before validity boundary - use previous day's index
      else
         f107_index = 0   ! After validity boundary - use current day's index
      end if

      ! ===== KP SELECTION =====
      if (sub_index >= 0) then
         ! Use Kp from current day (Fortran uses 1-based indexing)
         kp_out = flux_data%kp(sub_index + 1)
      else
         ! Need previous day's data
         call sw_data%get_flux_data(flux_data%mjd - 1.0_dp, flux_data_prev, sw_status)

         if (sw_status) then
            ! Use appropriate Kp from previous day
            ! sub_index is negative, so 8 + sub_index gives the right period
            kp_out = flux_data_prev%kp(8 + sub_index + 1)
         else
            ! If we can't get previous day, use first period of current day
            kp_out = flux_data%kp(1)
         end if
      end if

      ! ===== F10.7 DAILY VALUE (always uses previous day) =====
      if (f107_index == 0) then
         ! Current time is after measurement time - use yesterday's F10.7
         call sw_data%get_flux_data(flux_data%mjd - 1.0_dp, flux_data_prev, sw_status)
         if (sw_status) then
            f107_out = flux_data_prev%f107_obs
         else
            f107_out = flux_data%f107_obs  ! Fallback to current
         end if
      else
         ! Current time is before measurement time - use 2 days ago F10.7
         call sw_data%get_flux_data(flux_data%mjd - 2.0_dp, flux_data_2days_ago, sw_status)
         if (sw_status) then
            f107_out = flux_data_2days_ago%f107_obs
         else
            ! Try 1 day ago as fallback
            call sw_data%get_flux_data(flux_data%mjd - 1.0_dp, flux_data_prev, sw_status)
            if (sw_status) then
               f107_out = flux_data_prev%f107_obs
            else
               f107_out = flux_data%f107_obs  ! Last resort
            end if
         end if
      end if

      ! ===== F10.7a AVERAGE (centered 81-day average) =====
      ! if (f107_index == 0) then
      if (f107_index == 0) then
         if (use_gmat_bug) then
            !-----------------------------------------------------
            ! GMAT bug?: uses previous day's F10.7a even after 8am
            ! see: https://github.com/nasa/GMAT/issues/4
            call sw_data%get_flux_data(flux_data%mjd - 1.0_dp, flux_data_prev, sw_status)
            if (sw_status) then
               f107a_out = flux_data_prev%f107a_obs_ctr
            else
               f107a_out = flux_data%f107a_obs_ctr  ! Fallback
            end if
            !-----------------------------------------------------
         else
            ! Current time is after measurement time - use today's F10.7a
            f107a_out = flux_data%f107a_obs_ctr
         end if
      else
         ! Current time is before measurement time - use yesterday's F10.7a
         call sw_data%get_flux_data(flux_data%mjd - 1.0_dp, flux_data_prev, sw_status)
         if (sw_status) then
            f107a_out = flux_data_prev%f107a_obs_ctr
         else
            f107a_out = flux_data%f107a_obs_ctr  ! Fallback to current
         end if
      end if

   end subroutine prepare_flux_data