Fortran Programming Language :: Fortran: creating a table and interpolating values


	Fortran: creating a table and interpolating values 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. Many thanks for your time, smurray444 smurray444 wrote: > 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. (snip) Are you asking how to read such data in Fortran, how to do linear interpolation in two dimensions, or how to write it out again? Most printers can't print 65090 columns, and finding paper that wide is also hard. What do you want to do with the table once you have computed it? -- glen - Hide quoted text - 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 -- Paul van Delst Ride lots. CIMSS @ NOAA/NCEP/EMC Eddy Merckx