Format of PB2002 tectonic zones.grd

From BemlarGroup

Reading PB2002_tectonic_zones.grd with Fortran (for example):

The first two lines of the file are literally:

    0.000     0.100   360.000

  -90.000     0.100    90.000

which means that X (East longitude) goes from 0 to 360 in steps of 0.1, while Y (North latitude) goes from -90 to 90 in steps of 0.1. The convention is that X and Y ranges are always given from low to high.

However, data points that follow are listed as if English text were printed over a map.

The order of grid points in my PB2002_tectonic_zones.grd file is: Organized into rows, with West->East progression in each row, from 0.0 to 360.0 degrees East. Therefore the last value in each row is equal to the first value in each row. The first row is at +90.0N, and all values should be the same. The second row is at +89.9N. The last row is at -90.0N = 90 S, and all values should be the same.

Since each line of data values has 40 entries (followed by CR and LF bytes, according to DOS convention), the line breaks in the file do not correspond to rows. Therefore, it is easiest to use a single Fortran READ (with 2 implied DO loops) to input the whole array.

However, it is not required that you store the data in memory in the same order that it is found in the input file.

If you want your
INTEGER*2, DIMENSION (1801,3601):: array
to have rows numbered from 1 = N pole to 1801 = S pole
(as in the input file) then you could use:
READ (1, *)((array(i, j), j = 1, 3601), i = 1, 1801).

Or, if you want row number to be proportional to latitude, and also column number to be proportional to longitude, you could set up:
INTEGER*2, DIMENSION (-900:900, 0:3600) :: array
and then use
READ (1, *)((array(i, j), j = 0, 3600), i = 900, -900, -1)
In this case, you choose the correct grid point by simple logic:
IF (longitude < 0.0) longitude = longitude + 360.0
IF (longitude < 0.0) longitude = longitude + 360.0
IF (longitude > 360.0) longitude = longitude - 360.0
IF (longitude > 360.0) longitude = longitude - 360.0
jColumn = NINT(10 * longitude)
latitude = MIN(MAX(latitude, -90.0), 90.0)
iRow = NINT(10 * latitude)

Another usage that is possible in Fortran 90 is:
INTEGER*2, DIMENSION(1801, 3601):: array
READ (1, *) array
but I don't recommend this because it is too confusing! (Someone reading this program years from now would not be able to figure out the meaning of the row numbers.)

N.B. You can change INTEGER*2 to CHARACTER*1 in the above, but then you have to use a format of "(40(1X,A1))" instead of "*" in any READ statement.

Best wishes,
Peter Bird
UCLA
2007.10.08