Initialize the space weather module by reading a CSSI file
| Type | Intent | Optional | 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) |
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