sw_init Subroutine

private subroutine sw_init(me, filename, status)

Initialize the space weather module by reading a CSSI file

Type Bound

sw_data_type

Arguments

Type IntentOptional Attributes Name
class(sw_data_type), intent(inout) :: me
character(len=*), intent(in) :: filename

Path to CSSI space weather file

integer(kind=ip), intent(out) :: status

Output status (0=success, non-zero=error)


Calls

proc~~sw_init~~CallsGraph proc~sw_init sw_data_type%sw_init proc~date_to_mjd date_to_mjd proc~sw_init->proc~date_to_mjd proc~sw_cleanup sw_data_type%sw_cleanup proc~sw_init->proc~sw_cleanup proc~julian_day julian_day proc~date_to_mjd->proc~julian_day

Called by

proc~~sw_init~~CalledByGraph proc~sw_init sw_data_type%sw_init proc~jr_init jacchia_roberts_type%jr_init proc~jr_init->proc~sw_init

Source Code

   subroutine sw_init(me, filename, status)
      class(sw_data_type), intent(inout) :: me
      character(len=*), intent(in) :: filename !! Path to CSSI space weather file
      integer(ip), intent(out) :: status !! Output status (0=success, non-zero=error)

      integer(ip) :: unit, io_stat, i
      character(len=500) :: line
      integer(ip) :: year, month, day, bsrn, nd, kp_sum, kp_int(8), ap_int(8), ap_avg_int
      integer(ip) :: c9, isn, qual
      real(dp) :: cp, f107_adj, f107_adj_ctr81, f107_adj_lst81
      real(dp) :: f107_obs, f107_obs_ctr81
      real(dp) :: mjd_val
      integer(ip) :: n_lines
      logical :: in_data_section

      character(len=*),parameter :: fmt = '(I4,I3,I3,I5,I3,8I3,I4,8I4,I4,F4.1,I2,I4,F6.1,I2,5F6.1)' !! Format string for reading data lines

      status = 0

      ! Deallocate if already allocated
      call me%destroy()

      ! Open file
      open(newunit=unit, file=filename, status='old', action='read', iostat=io_stat)
      if (io_stat /= 0) then
         write(*,'(A,A)') 'ERROR: Cannot open space weather file: ', trim(filename)
         status = 1
         return
      end if

      ! First pass: count valid data records
      in_data_section = .false.
      n_lines = 0

      do
         read(unit, '(A)', iostat=io_stat) line
         if (io_stat /= 0) exit

         if (len_trim(line) == 0) cycle
         if (line(1:1) == '#') cycle

         ! Check for BEGIN sections (OBSERVED, DAILY_PREDICTED, MONTHLY_PREDICTED)
         if (index(line, 'BEGIN OBSERVED') > 0 .or. &
             index(line, 'BEGIN DAILY_PREDICTED') > 0 .or. &
             index(line, 'BEGIN MONTHLY_PREDICTED') > 0) then
            in_data_section = .true.
            cycle
         end if

         ! Check for END sections
         if (index(line, 'END OBSERVED') > 0 .or. &
             index(line, 'END DAILY_PREDICTED') > 0 .or. &
             index(line, 'END MONTHLY_PREDICTED') > 0) then
            in_data_section = .false.
            cycle
         end if

         ! Stop at monthly fit section
         if (index(line, 'BEGIN MONTHLY_FIT') > 0) exit

         if (in_data_section .and. len_trim(line) > 80) then
            n_lines = n_lines + 1
         end if
      end do

      if (n_lines == 0) then
         write(*,'(A)') 'ERROR: No data records found in file'
         status = 4
         close(unit)
         return
      end if

      ! Allocate arrays with exact size needed
      allocate(me%mjd(n_lines))
      allocate(me%f107_obs(n_lines))
      allocate(me%f107_adj(n_lines))
      allocate(me%f107a_obs_ctr(n_lines))
      allocate(me%f107a_adj_ctr(n_lines))
      allocate(me%kp(8, n_lines))
      allocate(me%ap_avg(n_lines))

      ! Rewind file for second pass
      rewind(unit)

      ! Second pass: read the data
      in_data_section = .false.
      me%n_records = 0

      do
         read(unit, '(A)', iostat=io_stat) line
         if (io_stat /= 0) exit

         ! Skip empty lines and comments
         if (len_trim(line) == 0) cycle
         if (line(1:1) == '#') cycle

         ! Check for BEGIN sections (OBSERVED, DAILY_PREDICTED, MONTHLY_PREDICTED)
         if (index(line, 'BEGIN OBSERVED') > 0 .or. &
             index(line, 'BEGIN DAILY_PREDICTED') > 0 .or. &
             index(line, 'BEGIN MONTHLY_PREDICTED') > 0) then
            in_data_section = .true.
            cycle
         end if

         ! Check for END sections
         if (index(line, 'END OBSERVED') > 0 .or. &
             index(line, 'END DAILY_PREDICTED') > 0 .or. &
             index(line, 'END MONTHLY_PREDICTED') > 0) then
            in_data_section = .false.
            cycle
         end if

         ! Stop at monthly fit section
         if (index(line, 'BEGIN MONTHLY_FIT') > 0) exit

         ! Parse data lines using fixed-width format
         if (in_data_section .and. len_trim(line) >= 130) then
            ! Column positions per CelesTrak documentation
            read(line, fmt, iostat=io_stat) &
                 year, month, day, bsrn, nd, &
                 kp_int(1), kp_int(2), kp_int(3), kp_int(4), &
                 kp_int(5), kp_int(6), kp_int(7), kp_int(8), kp_sum, &
                 ap_int(1), ap_int(2), ap_int(3), ap_int(4), &
                 ap_int(5), ap_int(6), ap_int(7), ap_int(8), ap_avg_int, &
                 cp, c9, isn, f107_adj, qual, f107_adj_ctr81, f107_adj_lst81, &
                 f107_obs, f107_obs_ctr81

            if (io_stat /= 0) cycle

            ! Convert date to MJD
            call date_to_mjd(year, month, day, mjd_val)

            me%n_records = me%n_records + 1

            ! Store data (Kp values are stored as Kp*10, convert back)
            me%mjd(me%n_records) = mjd_val
            do i = 1, 8
               me%kp(i, me%n_records) = real(kp_int(i), dp) / 10.0_dp
            end do
            me%ap_avg(me%n_records) = real(ap_avg_int, dp)
            me%f107_adj(me%n_records) = f107_adj
            me%f107a_adj_ctr(me%n_records) = f107_adj_ctr81
            me%f107_obs(me%n_records) = f107_obs
            me%f107a_obs_ctr(me%n_records) = f107_obs_ctr81
         end if
      end do

      close(unit)

      if (me%n_records == 0) then
         write(*,'(A)') 'ERROR: No valid data records found'
         status = 4
         return
      end if

      ! Set epoch range tracking
      me%historic_start = me%mjd(1)
      ! Note: add 1.0 to last epoch to reach end of day
      me%historic_end = me%mjd(me%n_records) + 1.0_dp

      ! Find where daily data ends and monthly data begins (gaps > 1 day)
      me%historic_daily_end = me%historic_end
      do i = 2, me%n_records
         if (me%mjd(i) - me%mjd(i-1) > 1.01_dp) then
            me%historic_daily_end = me%mjd(i-1)
            exit
         end if
      end do

      me%initialized = .true.

      write(*,'(A,I0,A)') 'Space weather data loaded: ', me%n_records, ' records'
      write(*,'(A,I0,A,I0)') '  (Allocated for ', n_lines, ' lines, successfully parsed ', me%n_records, ')'
      write(*,'(A,F12.2,A,F12.2)') '  Date range (MJD): ', me%mjd(1), ' to ', &
                                    me%mjd(me%n_records)

   end subroutine sw_init