PROGRAM XML_2_GRD

!  Utility program to read a CSEP-format seismicity forecast file
!  (3.7 GB, with 31 magnitude bins from 6.00+-0.05 to 8.95+),
!  and convert it to a single-magnitude-bin (all above threshold)
!  forecast in Peter Bird's GRD format.
!  This allows plotting with NeoKineMap, and also allows
!  testing with Kagan_2009_GJI_I_scores and/or PseudoCSEP.
!  Note that the spatial grid is unchanged; however,
!  the forecast is also transformed from specific earthquake counts
!  in a specific time window to earthquake rates densities
! (/m**2 /s) for an indefinite or open window.
!  By Peter Bird, UCLA, 2014.05.12.

    IMPLICIT NONE
    CHARACTER*12 :: duration_days_c12
    CHARACTER*79 :: GRD_file, line, XML_file
    INTEGER, PARAMETER :: nMagnitudeBins = 31
    INTEGER :: i, ios, iRow, j, jCol, k1, k2, n, nColumns, nRows, whole_steps
    REAL, PARAMETER :: Earth_radius_m = 6371000.
    REAL :: bottom_lat, &
          & d_lat, d_lon, duration_days_real, duration_s, &
          & EQ_count, EQ_rate_density, &
          & highest_threshold, &
          & lat, lat_max, lat_min, lon, lon_max, lon_min, lowest_threshold, &
          & magnitudeBinDimension, magnitudeMinValue, &
          & R2, &
          & s_per_year, strip_area_steradian, &
          & threshold_magnitude, top_lat
    REAL, DIMENSION(nMagnitudeBins) :: magnitude_lower_limit
    REAL, DIMENSION(:), ALLOCATABLE :: cell_area_m2 ! with one entry per row of seismicity matrix (N to S)
    DOUBLE PRECISION, PARAMETER :: unity = 1.0D0
    DOUBLE PRECISION, PARAMETER :: Pi = 3.141592654D0
    DOUBLE PRECISION :: dx, dy, exponent, fk1, fk2, fraction, integral, integral_converted, radians_per_degree, total_area_m2
    DOUBLE PRECISION, DIMENSION(nMagnitudeBins) :: earthquake_count ! for one cell, over duration_s, in magnitude bin #i = 1, ... , nMagnitudeBins
    DOUBLE PRECISION, DIMENSION(nMagnitudeBins) :: log10_earthquake_rate_density ! in one cell, for indefinite duration, in magnitude bin #i = 1, ... , nMagnitudeBins
    DOUBLE PRECISION, DIMENSION(:, :), ALLOCATABLE :: seismicity

    radians_per_degree = Pi / 180.0D0
    s_per_year = 365.25 * 24. * 60. * 60.

    WRITE (*, "(' PROGRAM XML_2_GRD')")
    WRITE (*, "(' ')")
    WRITE (*, "('   Utility program to read a CSEP-format seismicity forecast file')")
    WRITE (*, "('   (3.7 GB, with ',I2,' magnitude bins from 6.00+-0.05 to 8.95+),')") nMagnitudeBins
    WRITE (*, "('   and convert it to a single-magnitude-bin (all above threshold)')")
    WRITE (*, "('   forecast in Peter Bird''s GRD format.')")
    WRITE (*, "('   This allows plotting with NeoKineMap, and also allows')")
    WRITE (*, "('   testing with Kagan_2009_GJI_I_scores and/or PseudoCSEP.')")
    WRITE (*, "('   Note that the spatial grid is unchanged; however,')")
    WRITE (*, "('   the forecast is also transformed from specific earthquake counts')")
    WRITE (*, "('   in a specific time window to earthquake rates densities')")
    WRITE (*, "('  (/m**2 /s) for an indefinite or open window.')")
    WRITE (*, "('   By Peter Bird, UCLA, 2014.05.12.')")

 10 WRITE (*, *)
    WRITE (*, "(' Enter name of existing XML file: '\)")
    READ (*, "(A)") XML_file
    WRITE (*, "(' Note that each OPEN takes a minute or two, for buffering...')")
    OPEN (UNIT = 1, FILE = TRIM(XML_file), STATUS = "OLD", PAD = "YES", IOSTAT = ios)
    IF (ios /= 0) THEN
        WRITE (*, "(' ERROR: This file could not be found/opened (in current directory).')")
        CALL Pause()
        GO TO 10
    END IF

    WRITE (*, *)
    WRITE (*, "(' HEADER LINES in this XML file:')")
    WRITE (*, "(' -----------------------------------------------------------------------')")
    DO i = 1, 16 ! expecting 16 lines of header information...
        READ (1, "(A)") line
        WRITE (*, "(' ',A)") TRIM(line)

        IF (i == 11) THEN ! try to get duration_days_c12
            IF (line(1:22) == "    <forecastDuration>") THEN
                duration_days_c12 = line(23:34)
                n = SCAN(duration_days_c12, '<') ! check for intrusion of </forecastDuration>
                IF (n > 0) duration_days_c12(n:n) = ' ' ! create a space to permit READ(duration_days_c12, *)
                READ (duration_days_c12, *, IOSTAT = ios) duration_days_real
                IF (ios /= 0) THEN
                    WRITE (*, *)
                    WRITE (*, "(' ERROR: Could not get duration_days_real from: ',A)") duration_days_c12
                    CALL Pause()
                    STOP
                END IF
                duration_s = duration_days_real * 24. * 60. * 60.
            ELSE
                WRITE (*, *)
                WRITE (*, "(' ERROR in file format: Expecting <forecastDuration>, but found:')")
                WRITE (*, "(A)") TRIM(line)
                CALL Pause()
                STOP
            END IF
        END IF ! line #11: duration_day_c1

        IF (i == 12) THEN ! try to get d_lat, d_lon
            IF (line(1:28) == "    <cellDimension latRange=") THEN
                line = line(30:79) ! delete prefix, including leading "
                n = SCAN(line, '"') ! find close-quote after d_lat
                line(n:n) = ' '     ! replace with a space
                READ (line, *, IOSTAT = ios) d_lat
                IF (ios /= 0) THEN
                    WRITE (*, *)
                    WRITE (*, "(' ERROR: Could not get d_lat from start of: ',A)") TRIM(line)
                    CALL Pause()
                    STOP
                END IF
                n = SCAN(line, '=')
                IF (n == 0) THEN
                    WRITE (*, *)
                    WRITE (*, "(' ERROR: Could not find second ""="", to get d_lon from interior of: ',A)") TRIM(line)
                    CALL Pause()
                    STOP
                END IF
                line = line(n+2:79) ! cut out the string introducting d_lon
                n = SCAN(line, '"') ! find close-quote after d_lon
                line(n:n) = ' '     ! replace with a space
                READ (line, *, IOSTAT = ios) d_lon
                IF (ios /= 0) THEN
                    WRITE (*, *)
                    WRITE (*, "(' ERROR: Could not get d_lon from start of edited string: ',A)") TRIM(line)
                    CALL Pause()
                    STOP
                END IF
            ELSE
                WRITE (*, *)
                WRITE (*, "(' ERROR in file format: Expecting <cellDimension latRange=, but found:')")
                WRITE (*, "(A)") TRIM(line)
                CALL Pause()
                STOP
            END IF
        END IF ! line #12: d_lat, d_lon

        IF (i == 13) THEN ! line #13: magnitudeBinDimension
            IF (line(1:27) == "    <magnitudeBinDimension>") THEN
                line = line(28:79) ! cut out the prefix
                n = SCAN(line, '<') !find start of </magnitudeBinDimension>
                line(n:n) = ' '     !convert it to a space
                READ (line, *, IOSTAT = ios) magnitudeBinDimension
                IF (ios /= 0) THEN 
                    WRITE (*, *)
                    WRITE (*, "(' ERROR: Could not get magnitudeBinDimension from start of string: '/' ', A)") TRIM(line)
                    CALL Pause()
                    STOP
                END IF
            ELSE
                WRITE (*, *)
                WRITE (*, "(' ERROR: Expecting <magnitudeBinDimension> but found: '/' ', A)") TRIM(line)
                CALL Pause()
                STOP
            END IF
        END IF ! line #13: magnitudeBinDimension

        IF (i == 14) THEN ! line #14: magnitudeminValue
            IF (line(1:23) == "    <magnitudeMinValue>") THEN
                line = line(24:79) ! cut out the prefix
                n = SCAN(line, '<') !find start of </magnitudeMinValue>
                line(n:n) = ' '     !convert it to a space
                READ (line, *, IOSTAT = ios) magnitudeMinValue
                IF (ios /= 0) THEN 
                    WRITE (*, *)
                    WRITE (*, "(' ERROR: Could not get magnitudeMinValue from start of string: '/' ', A)") TRIM(line)
                    CALL Pause()
                    STOP
                END IF
            ELSE
                WRITE (*, *)
                WRITE (*, "(' ERROR: Expecting <magnitudeMinValue> but found: '/' ', A)") TRIM(line)
                CALL Pause()
                STOP
            END IF
        END IF ! line #14: magnitudeMinValue

    END DO
    WRITE (*, "(' -----------------------------------------------------------------------')")
    WRITE (*, "(' Inferred forecast duration = ',ES12.5,' seconds.')") duration_s
    WRITE (*, "(' Inferred d_lat = ',ES12.5,', d_lon = ',ES12.5)") d_lat, d_lon
    WRITE (*, "(' Inferred magnitudeBinDimension = ',F6.3,', magnitudeMinValue = ',F6.3)") magnitudeBinDimension, magnitudeMinValue
    magnitude_lower_limit(1) = magnitudeMinValue
    DO i = 2, nMagnitudeBins
        magnitude_lower_limit(i) = magnitude_lower_limit(i-1) + magnitudeBinDimension
    END DO
    CALL Pause()

    !scan file for ranges of lon and lat:
    WRITE (*, *)
    WRITE (*, "(' Scanning XML file for extreme values of lon and lat (SLOW)...')") 
    lat_max = -9999.
    lat_min = +9999.
    lon_max = -9999.
    lon_min = +9999.
    scanning: DO
        READ (1, "(A)", IOSTAT = ios) line
        IF (ios /= 0) EXIT scanning
        IF (line(1:16) == "      <cell lon=") THEN
            line = line(18:79) ! cutting out prefix, and also the leading "
            n = SCAN(line, '"')
            line(n:n) = ' '    ! replace closing " with a space
            READ (line, *, IOSTAT = ios) lon 
            IF (ios /= 0) THEN
                WRITE (*, *)
                WRITE (*, "(' ERROR: Could not get lon from start of altered line: '/' ', A)") TRIM(line)
                CALL Pause()
                STOP
            END IF
            n = SCAN(line, '=')   ! find end of prefix for lat
            line = line((n+2):79) ! cut out prefix, and also prior "
            n = SCAN(line, '"')   ! find closing "
            line(n:n) = ' '       ! convert it to a space
            READ (line, *, IOSTAT = ios) lat
            IF (ios /= 0) THEN
                WRITE (*, *)
                WRITE (*, "(' ERROR: Could not get lat from start of altered line: '/' ', A)") TRIM(line)
                CALL Pause()
                STOP
            END IF
            lat_max = MAX(lat_max, lat)
            lat_min = MIN(lat_min, lat)
            lon_max = MAX(lon_max, lon)
            lon_min = MIN(lon_min, lon)
        END IF ! lon/lat line
    END DO scanning
    WRITE (*, "(' Inferred lon_min = ',F10.4,', lon_max = ',F10.4)") lon_min, lon_max
    WRITE (*, "(' Inferred lat_min = ',F10.4,', lat_max = ',F10.4)") lat_min, lat_max
    CALL Pause()
    CLOSE (1)

    !create grid to store eventual GRD seismicity rate densities:
    nRows = 1 + NINT((lat_max - lat_min) / d_lat)
    nColumns = 1 + NINT((lon_max - lon_min) / d_lon)
    WRITE (*, *)
    WRITE (*, "(' Allocating seismicity grid (REAL*8) of ',I6,' rows and ',I6,' columns...')") nRows, nColumns
    ALLOCATE ( seismicity(nRows, nColumns) )
    seismicity = 0.0D0 ! whole matrix

    !compute areas of cells:
    ALLOCATE ( cell_area_m2(nRows) )
    R2 = Earth_radius_m**2
    total_area_m2 = 0.0D0
    DO i = 1, nRows ! N to S
        top_lat = lat_max + 0.5 * d_lat - (i-1) * d_lat
        bottom_lat = top_lat - d_lat
        strip_area_steradian = 2 * Pi * (unity - DCOS((90.0D0 - bottom_lat) * radians_per_degree)) &
                           & - 2 * Pi * (unity - DCOS((90.0D0 - top_lat   ) * radians_per_degree))
        cell_area_m2(i) = strip_area_steradian * R2 * ((lon_max + d_lon - lon_min) / 360.0) / nColumns
        total_area_m2 = total_area_m2 + nColumns * cell_area_m2(i)
    END DO
    WRITE (*, "(' ')") 
    WRITE (*, "(' Total grid area = ',ES12.5,' square meters,')") total_area_m2

    !re-open XML file and skip the headers this time:
    WRITE (*, *)
    WRITE (*, "(' Closing and then re-opening the big XML file.')")
    WRITE (*, "(' Note that each OPEN takes a minute or two, for buffering...')")
    OPEN (UNIT = 1, FILE = TRIM(XML_file), STATUS = "OLD", PAD = "YES", IOSTAT = ios)
    DO i = 1, 16
        READ (1, "(A)") line
    END DO

    !select threshold magnitude for GRD file:
    lowest_threshold = magnitudeMinValue - 2.0 * magnitudeBinDimension
    highest_threshold = magnitude_lower_limit(nMagnitudeBins) + 2.0 * magnitudeBinDimension
 90 WRITE (*, *)
    WRITE (*, "(' Enter threshold magnitude for GRD file: '\)") 
    READ (*, *) threshold_magnitude
    IF ((threshold_magnitude < lowest_threshold).OR.(threshold_magnitude > highest_threshold)) THEN
        WRITE (*, *)
        WRITE (*, "(' ERROR: Please enter a threshold_magnitude value from ',F6.3,' to ',F6.3,';')") lowest_threshold, highest_threshold 
        WRITE (*, "('        because off-end extrapolation is limited to 2 bin-widths.')") 
        CALL Pause()
        GO TO 90
    END IF
    whole_steps = (threshold_magnitude - magnitudeMinValue) / magnitudeBinDimension ! truncates to: whole-steps (ignoring suplus fraction) to right from lower-limit of lowest bin
    k1 = whole_steps + 1     ! may be negative, at this point....
    k1 = MAX(1, k1)
    k1 = MIN((nMagnitudeBins-1), k1)  ! within range from 1 ~ (nMagnitudeBins-1)
    k2 = k1 + 1                       ! within range from 2 ~ nMagnitudeBins
    fraction = ((threshold_magnitude - magnitudeMinValue) / magnitudeBinDimension) - (k1 - 1.0) ! usually in 0.0 ~ 1.0, but not always!
    fk1 = (1.0D0 - fraction) ! usually in 0.0 ~ 1.0, but not always!
    fk2 = fraction           ! usually in 0.0 ~ 1.0, but not always!
    WRITE (*, *)
    WRITE (*, "(' Log_10(seismicity) at threshold ',F6.3,' will get weight of ',F7.4)") magnitude_lower_limit(k1), fk1
    WRITE (*, "(' Log_10(seismicity) at threshold ',F6.3,' will get weight of ',F7.4)") magnitude_lower_limit(k2), fk2

    !select name for output GRD file:
    WRITE (*, *)
    WRITE (*, "(' Enter filename for new, output GRD file: '\)")
    READ (*, "(A)") GRD_file
    OPEN (UNIT = 2, FILE = TRIM(GRD_file)) ! unconditional OPEN; overwrites any existing file!

    WRITE (*, *)
    WRITE (*, "(' Starting primary working loop; SLOW ...')")

    rescanning: DO

        READ (1, "(A)", IOSTAT = ios) line
        IF (ios /= 0) EXIT rescanning

        IF (line(1:16) == "      <cell lon=") THEN ! Note: This should always be true, I think.

            !extract the (lon, lat) information:
            line = line(18:79) ! cutting out prefix, and also the leading "
            n = SCAN(line, '"')
            line(n:n) = ' '    ! replace closing " with a space
            READ (line, *, IOSTAT = ios) lon 
            IF (ios /= 0) THEN
                WRITE (*, *)
                WRITE (*, "(' ERROR: Could not get lon from start of altered line: '/' ', A)") TRIM(line)
                CALL Pause()
                STOP
            END IF
            n = SCAN(line, '=')   ! find end of prefix for lat
            line = line((n+2):79) ! cut out prefix, and also prior "
            n = SCAN(line, '"')   ! find closing "
            line(n:n) = ' '       ! convert it to a space
            READ (line, *, IOSTAT = ios) lat
            IF (ios /= 0) THEN
                WRITE (*, *)
                WRITE (*, "(' ERROR: Could not get lat from start of altered line: '/' ', A)") TRIM(line)
                CALL Pause()
                STOP
            END IF

            !determine target cell of matrix:
            iRow = 1 + NINT((lat_max - lat) / d_lat)
            jCol = 1 + NINT((lon - lon_min) / d_lon)

            !read earthquake counts in this cell, binned, and per forecast window (e.g., per year):
            DO i = 1, nMagnitudeBins
                READ (1, *, IOSTAT = ios) earthquake_count(i)
                IF (ios /= 0) THEN
                    WRITE (*, *)
                    WRITE (*, "(' ERROR: Could not read earthquake_count(',I2,') of row ',I6,', column ',I6,' at lat = ',F10.4,', lon = ',F10.4)") i, iRow, jCol, lat, lon
                    CALL Pause()
                    STOP
                END IF
            END DO
            READ (1, *) ! </cell>

            !recast earthquake numbers as cumulative:
            DO i = (nMagnitudeBins-1), 1, -1
                earthquake_count(i) = earthquake_count(i) + earthquake_count(i + 1)
            END DO

            !compute log10(cumulative earthquake rate densities) for this one cell:
            DO i = 1, nMagnitudeBins
                log10_earthquake_rate_density(i) = LOG10(earthquake_count(i) / (duration_s * 1.0D0 * cell_area_m2(iRow)))
            END DO

            !interpolate, using straight-line on a Gutenberg-Richter plot (of log rate vs. magnitude):
            exponent = fk1 * log10_earthquake_rate_density(k1) + &
                     & fk2 * log10_earthquake_rate_density(k2)
            seismicity(iRow, jCol) = 10.0D0**exponent

        END IF ! line had (lon, lat)

    END DO rescanning
    CLOSE(1) ! input XML file

    !Output the seismicity (EQs/m^2/s) .grd file:
    WRITE (*, *)
    WRITE (*, "(' Writing seismicity (EQs/m^2/s) .grd file and computing area-integral...')")
    WRITE (2, "(3F11.5,' = lon_min, d_lon, lon_max  Gridded long-term seismicity above magnitude ',F6.3,' in events/square-meter/s')") &
         & lon_min, d_lon, lon_max, threshold_magnitude
    WRITE (2, "(3F11.5,' = lat_min, d_lat, lat_max  from program XML_2_GRD, P. Bird, UCLA, 2014')") &
         & lat_min, d_lat, lat_max
    WRITE (2, "(1P,10D10.2)") ((seismicity(i, j), j = 1, nColumns), i = 1, nRows)
    !Note: Line breaks in output file have no relation to the logical shape of this array.
    !      Programs that read .grd files must take this into account, by reading the whole
    !      array with a single READ (as it was written), and not row-by-row.
    !-----------------------------------------------------------------------------------
   
    !Integrate global shallow seismicity: 

    !Integrate seismicity over grid area:
    integral = 0.0D0 ! (just initializing)
    DO i = 1, nRows
        DO j = 1, nColumns
            integral = integral + seismicity(i, j) * cell_area_m2(i)
        END DO
    END DO
    integral_converted = integral * 100.0 * s_per_year ! now, earthquakes per century
    WRITE (*, "(/' Area integral of seismicity =' &
              & /' ',ES11.4,' earthquakes/s = ',ES11.4,' earthquakes/century.')") integral, integral_converted
    WRITE (2, "(  'Area integral of seismicity =' &
              &     /ES11.4,' earthquakes/s = ',ES11.4,' earthquakes/century.')") integral, integral_converted
    CLOSE (UNIT = 2, DISP = "KEEP") ! seismicity.grd file

    !------------------------------------------------------------    

    WRITE (*, *)
    WRITE (*, "(' Job completed.')")
    CALL Pause()

CONTAINS

    SUBROUTINE Pause ()
        IMPLICIT NONE
        CHARACTER*1 c1
        WRITE (*,"(' Press [Enter] to continue...'\)")
        READ (*,"(A)") c1
    END SUBROUTINE Pause

END PROGRAM XML_2_GRD