PROGRAM GSRM2p2_to_color_GRD
    !Written by Peter Bird, UCLA, 2014.04.10 (as "GSRM2p1_to_color_GRD.f90"); name changed 2019.08.19,
    !when I also changed the input data file name to "GSRM_average_strain.txt".
    !Another minor change was REAL --> REAL*8 throughout
    IMPLICIT NONE
    INTEGER :: i, ios, j, nColumns, nRows
    REAL*8 :: exx, exy, eyy, lat, lon
    !-----------------------------------------------------------------------------------------
    !Definitions specific to GSRM2.2 input file "GSRM_average_strain.txt" of Kreemer [p.c., 2019],
    !and to its internal (filled-out with zeros) representation GSRM2_grid:
    INTEGER :: GSRM2_rows, GSRM2_columns ! to be computed from data below:
    REAL*8, PARAMETER:: GSRM2_lat_min =  -89.800D0, GSRM2_lat_max =  89.800D0, GSRM2_d_lat = 0.200D0, &
                      & GSRM2_lon_min = -179.875D0, GSRM2_lon_max = 179.875D0, GSRM2_d_lon = 0.250D0, &
                      & GSRM2_central_lon = 0.0D0
                      ! in decimal degrees of East longitude & North latitude.
    REAL*8, PARAMETER:: strain_rate_units = 1.0D-9 / 3.15576D7 ! "units of 10^-9/year"
    REAL*8, DIMENSION(:, :, :), ALLOCATABLE :: GSRM2_grid
    !Note that longitude range in this dataset is from GSRM_lon_min to (GSRM_lon_min + 360.0D0).
    !Each strain-rate tensor in this grid (and each epicenter_rate_density point computed from it)
    !represents the average in a rectangular region centered on the grid point.
    !-----------------------------------------------------------------------------------------
    REAL*8, DIMENSION(:, :), ALLOCATABLE :: GSRM2_color
    DOUBLE PRECISION :: second_invariant

    nColumns = NINT(360.0D0 / GSRM2_d_lon)
    nRows = 1 + NINT((GSRM2_lat_max - GSRM2_lat_min) / GSRM2_d_lat)
    ALLOCATE ( GSRM2_grid(nRows, nColumns, 3 ) )
    GSRM2_grid = 0.0D0 ! whole 3-D array
    ALLOCATE ( GSRM2_color(nRows, nColumns) )
    GSRM2_color = 0.0D0 ! whole matrix

    !Read in Corne Kreemer's GSRM2 strain rates:
    OPEN (UNIT = 1, FILE = "GSRM_average_strain.txt", STATUS = "OLD", IOSTAT = ios)
    IF (ios /= 0) THEN
        WRITE (*, "(' ERROR: GSRM_average_strain.txt not found.  Fix, and try again...')")
        CALL Pause()
        STOP
    END IF
    DO i = 1, 25 ! lines of headers
        READ (1, *) ! read (& discard) header
    END DO
    scanning: DO
        READ (1, *, IOSTAT = ios) lat, lon, exx, eyy, exy
        IF (ios < 0) EXIT scanning ! hit the EOF
        !-------------------------------------------------------------
        !Line added specifically for this program:
        second_invariant = SQRT(exx**2 + eyy**2 + 2*(exy**2))
        !-------------------------------------------------------------
        exx = exx * strain_rate_units
        eyy = eyy * strain_rate_units
        exy = exy * strain_rate_units
        i = 1 + NINT((GSRM2_lat_max - lat) / GSRM2_d_lat)
        j = 1 + NINT((lon - GSRM2_lon_min) / GSRM2_d_lon)
        GSRM2_grid(i, j, 1) = exx ! = e_E-W
        GSRM2_grid(i, j, 2) = eyy ! = e_N-S
        GSRM2_grid(i, j, 3) = exy ! = e_NE-SW
        !-------------------------------------------------------------
        !Lines added specifically for this program:
        IF (second_invariant > 0.0D0) THEN
            GSRM2_color(i, j) = LOG10(second_invariant)
        ELSE
            GSRM2_color(i, j) = 0.0D0
        END IF
        !-------------------------------------------------------------
    END DO scanning
    CLOSE (1) ! GSRM_average_strain.txt, from Corne Kreemer

    !Output grid of colors for plotting with NeoKineMap (mosaic type 2; zeros left white):
    OPEN (UNIT = 2, FILE = "GSRM2.2_color.grd")
        WRITE (2, "(3F11.5,3X,'GSRM2.2: log_10(second-invariant(2-D strain-rate), in nano/year);')") GSRM2_lon_min, GSRM2_d_lon, GSRM2_lon_max
        WRITE (2, "(3F11.5,3X,'intended for plotting with NeoKineMap, as a type-2 mosaic (with zeros left white).')") GSRM2_lat_min, GSRM2_d_lat, GSRM2_lat_max
        WRITE (2, "(12F8.4)") ((GSRM2_color(i, j), j = 1, nColumns), i = 1, nRows)
    CLOSE (2)

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

CONTAINS

    SUBROUTINE Pause()
       IMPLICIT NONE
       WRITE (*,"(' Press [Enter]...'\)")
       READ  (*,*)
    END SUBROUTINE Pause

END PROGRAM GSRM2p2_to_color_GRD