w3ft01.f

Go to the documentation of this file.

C> @file
C> @brief Interpolate values in a data field.
C> @author James McDonell @date 1984-06-27
 
C> For a given grid coordinate in a data array, estimates
C> a data value for that point using either a linear or quadratic
C> interpolation method.
C>
C> ### Program History Log:
C> Date | Programmer | Commment
C> -----|------------|---------
C> 1984-06-27 | James McDonell | Initial
C> 1989-11-01 | Ralph Jones | Change to cray cft77 fortran
C>
C> @param[in] STI Real*4 i grid coordinate of the point for which
C> an interpolated value is desired.
C> @param[in] STJ Real*4 j grid coordinate of the point for which
C> an interpolated value is desired.
C> @param[in] FLD Real*4 size(ii,jj) data field.
C> @param[in] II Integer*4 number of columns in 'fld'.
C> @param[in] JJ Integer*4 number of rows in 'fld'.
C> @param[in] NCYCLK Integer*4 code to specify if grid is cyclic or
C> not:
C> - = 0 Non-cyclic in ii, non-cyclic in jj
C> - = 1 Cyclic in ii, non-cyclic in jj
C> - = 2 Cyclic in jj, non-cyclic in ii
C> - = 3 Cyclic in ii, cyclic in jj
C> @param[in] LIN Integer*4 code specifying interpolation method:
C> - = 1 Linear interpolation
C> - .NE.1 Quadratic interpolation
C> @param[out] HI Real*4 data field value at (sti,stj) obtained
C> by interpolation.
C>
C> @author James McDonell @date 1984-06-27
      SUBROUTINE w3ft01(STI,STJ,FLD,HI,II,JJ,NCYCLK,LIN)
C
      REAL    ERAS(4)
      REAL    FLD(II,JJ)
      REAL    JY(4)
C
      i     = sti
      j     = stj
      fi    = i
      fj    = j
      xdeli = sti - fi
      xdelj = stj - fj
      ip2   = i + 2
      im1   = i - 1
      ip1   = i + 1
      jy(4) = j + 2
      jy(1) = j - 1
      jy(3) = j + 1
      jy(2) = j
      xi2tm = 0.0
      xj2tm = 0.0
      IF (lin.NE.1) THEN
        xi2tm = xdeli * (xdeli - 1.0) * 0.25
        xj2tm = xdelj * (xdelj - 1.0) * 0.25
      ENDIF
      IF ((i.LT.2).OR.(j.LT.2))       GO TO 10
      IF ((i.GT.ii-3).OR.(j.GT.jj-3)) GO TO 10
C
C     QUADRATIC (LINEAR TOO) OK W/O FURTHER ADO SO GO TO 170
C
      GO TO 170
C
   10 CONTINUE
        icyclk = 0
        jcyclk = 0
        IF (ncyclk) 20,120,20
C
   20 CONTINUE
        IF (ncyclk / 2 .NE. 0) jcyclk = 1
        IF (ncyclk .NE. 2)     icyclk = 1
        IF (icyclk) 30,70,30
C
   30 CONTINUE
        IF (i.EQ.1)      GO TO 40
        IF (i.EQ.(ii-1)) GO TO 50
        ip2 = i + 2
        im1 = i - 1
        GO TO 60
C
   40 CONTINUE
        ip2 = 3
        im1 = ii - 1
        GO TO 60
C
   50 CONTINUE
        ip2 = 2
        im1 = ii - 2
C
   60 CONTINUE
        ip1 = i + 1
C
   70 CONTINUE
        IF (jcyclk) 80,120,80
C
   80 CONTINUE
        IF (j.EQ.1)      GO TO 90
        IF (j.EQ.(jj-1)) GO TO 100
        jy(4) = j + 2
        jy(1) = j - 1
        GO TO 110
C
   90 CONTINUE
        jy(4) = 3
        jy(1) = jj - 1
        GO TO 110
C
  100 CONTINUE
        jy(4) = 2
        jy(1) = jj - 2
C
  110 CONTINUE
        jy(3) = j + 1
        jy(2) = j
C
  120 CONTINUE
        IF (lin.EQ.1) GO TO 160
        IF (icyclk) 140,130,140
C
  130 CONTINUE
        IF ((i.LT.2).OR.(i.GE.(ii-1)))  xi2tm = 0.0
C
  140 CONTINUE
        IF (jcyclk) 160,150,160
C
  150 CONTINUE
        IF ((j.LT.2).OR.(j.GE.(jj-1)))  xj2tm = 0.0
C
  160 CONTINUE
C
C.....DO NOT ALLOW POINT OFF GRID,CYCLIC OR NOT
C
        IF (i.LT.1)   i   = 1
        IF (ip1.LT.1) ip1 = 1
        IF (ip2.LT.1) ip2 = 1
        IF (im1.LT.1) im1 = 1
C
C.....DO NOT ALLOW POINT OFF GRID,CYCLIC OR NOT
C
        IF (i.GT.ii)   i   = ii
        IF (ip1.GT.ii) ip1 = ii
        IF (ip2.GT.ii) ip2 = ii
        IF (im1.GT.ii) im1 = ii
C
  170 CONTINUE
      DO 180 k = 1,4
        j1 = jy(k)
C
C.....DO NOT ALLOW POINT OFF GRID,CYCLIC OR NOT
C
        IF (j1.LT.1)  j1 = 1
        IF (j1.GT.jj) j1 = jj
        eras(k) = (fld(ip1,j1) - fld(i,j1)) * xdeli + fld(i,j1) +
     &  (fld(im1,j1) - fld(i,j1) - fld(ip1,j1) + fld(ip2,j1)) * xi2tm
  180 CONTINUE
C
      hi = eras(2) + (eras(3) - eras(2)) * xdelj + (eras(1) -
     &     eras(2) -  eras(3) + eras(4)) * xj2tm
C
      RETURN
      END