w3ft03.f

Go to the documentation of this file.

C> @file
C> @brief A point interpolater.
C> @author James Howcroft @date 1979-02-15
 
C> Do either bilinear or biquadratic interpolation for a
C> point within a two-dimensional data array.
C>
C> ### Program History Log:
C> Date | Programmer | Comment
C> -----|------------|--------
C> 1979-02-15 | James Howcroft | Initial.
C> 1989-01-25 | Ralph Jones | Change to microsoft fortran 4.10.
C> 1990-06-12 | Ralph Jones | Change to sun fortran 1.3.
C> 1991-03-30 | Ralph Jones | Convert to silicongraphics fortran.
C> 1993-03-29 | Ralph Jones | Add save statement.
C> 1996-07-01 | Ralph Jones | Compile on cray.
C>
C> @param[in] FL Real*4 two-dimensional cartesian array of data.
C> @param[in] MAXI Integer*4 i-dimension of fl.
C> @param[in] MAXJ Integer*4 j-dimension of fl.
C> @param[in] STI Real*4 i-coordinate to which a value is to be
C> interpolated.
C> @param[in] STJ Real*4 j-coordinate to which a value is to be
C> interpolated.
C> @param ITYPE
C> @param[out] HI Real*4 interpolated output value.
C>
C> @remark No error checks are made. it is left for the user to
C> determine that the point for which interpolation is desired
C> lies within the grid.
C>
C> @author James Howcroft @date 1979-02-15
      SUBROUTINE w3ft03(FL,HI,STI,STJ,MAXI,MAXJ,ITYPE)
C
       REAL            FL(MAXI,MAXJ)
       REAL            E (4)
C
       SAVE
C
       i  = sti
       j  = stj
       di = i
       dj = j
       di = sti - di
       dj = stj - dj
C
       hi = 0.
C                      TEST FOR POINT OFF GRID.
       IF (i.LT.1 .OR. i.GT.maxi) GO TO 300
       IF (j.LT.1 .OR. j.GT.maxj) GO TO 300
       IF (itype .NE. 2) GO TO 100
C                      DO BILINEAR IF POINT IS BETWEEN ULTIMATE AND
C                      PENULTIMATE ROWS, WHERE BIQUAD NOT POSSIBLE.
       IF (i.LT.2 .OR. i.GT.(maxi-1)) GO TO 100
       IF (j.LT.2 .OR. j.GT.(maxj-1)) GO TO 100
       GO TO 200
C
C        BILINEAR.
  100  CONTINUE
         hi = fl(i  ,j  )*(1.-di)*(1.-dj)  + fl(i+1,j  )*di*(1.-dj)
     &      + fl(i  ,j+1)*(1.-di)*dj       + fl(i+1,j+1)*di*dj
         GO TO 300
C
  200  CONTINUE
C        BIQUADRATIC.
       di2 = di*(di-1.)*.25
       dj2 = dj*(dj-1.)*.25
       j1  = j - 1
       DO 250  k=1,4
         e(k) = fl(i  ,j1)*(1.-di-di2) + fl(i+1,j1)*(di-di2)
     &        + (fl(i-1,j1) + fl(i+2,j1))*di2
         j1 = j1 + 1
  250  CONTINUE
       hi = e(2)*(1.-dj-dj2) + e(3)*(dj-dj2) + (e(1) + e(4))*dj2
C
  300  CONTINUE
       RETURN
       END