smurray444 wrote:

> Hi there,

> I have a global dataset (text file) where for every 105 values

> (currently formatted as a continuous string), the 1st represents

> longitude, the 2nd = latitude, the 3rd is a 'spinup' value which can

> be ignored, and the 4th to 105th values are hydrological runoff values

> for the years 1901 to 2002.

> I've been trying to create a table where longitude values are listed

> along the top row (x-axis) and latitude values are listed in the far

> left column (y-axis), with the corresponding runoff values for each

> year located in a given grid - therefore there will be 102 grids

> created (one for each year). The longitude values are currently not

> ordered and will need to be sorted (low to high) into position along

> the x-axis.

> This would result in a table of data values (6509 in each dimension)

> which represents global point values at a resolution of 1.5-degrees...

> I want to perform a linear interpolation to change the resolution of

> this data to 5-minutes. This will require 18 new values to be created

> above, below and to both sides of each existing data value in the

> current table (except in instances where a value is not surrounded by

> four values, eg. at the edges of the table, where only values

> available should be used).

> I'm new to Fortran, so if anyone has any helpful tips, advice, or

> could create some Fortran code to achieve this, then I'd be very

> grateful.

Hello,

Sounds a bit like a homework problem, but I'll give you the benefit of the doubt. :o)

Here's a brain-dead simple 1-d interpolator:

SUBROUTINE interp_1D(y, xlp, & ! Input

y_int ) ! Output

! Arguments

REAL(fp), INTENT(IN) :: y(:), xlp(:)

REAL(fp), INTENT(OUT) :: y_int

! Local variables

INTEGER :: i

! Perform interpolation

y_int = ZERO

DO i = 1,NPTS

y_int = y_int + xlp(i)*y(i)

END DO

END SUBROUTINE interp_1D

and a corresponding 2-d interpolator that just reuses the above.

SUBROUTINE interp_2D(z, xlp, ylp, & ! Input

z_int ) ! Output

! Arguments

REAL(fp), INTENT(IN) :: z(:,:), xlp(:), ylp(:)

REAL(fp), INTENT(OUT) :: z_int

! Local variables

INTEGER :: i

REAL(fp) :: a(NPTS)

! Interpolate z in x dimension for all y

DO i = 1,NPTS

CALL interp_1D(z(:,i),xlp,a(i))

END DO

! Interpolate z in y dimension

CALL interp_1D(a,ylp,z_int)

END SUBROUTINE interp_2D

The NPTS parameter is the order of the interpolation you want to perform, plus 1. E.g. for

linear (order 1), npts is thus 2, etc.

The xlp and ylp inputs are the interpolating polynomials for your input x and y. I'll

leave that bit up to you since for linear interpolation it's only 2 lines of code. Since

you're redoing the interpolation for every year, you only have to compute them once

(assuming your input data are all on the same x- and y-axes) and then just reuse them.

But, I have to say, doing a linear interpolation from 1.5deg global grids to a 5' grid

doesn't sound too useful. Of course, it depends on your input data (which may be mostly

linear to start with) but at the very least try higher orders (or a spline) to see how

they munge your data.

cheers,

paulv