w3ft07.f

Go to the documentation of this file.

C> @file
C> @brief Transform gridpoint fld by interpolation.
C> @author McDonell & Howcroft @date 1974-09-01
 
C> Transforms data contained in a given grid array
C> by translation, rotation about a common point and dilatation
C> in order to create a new grid array according to specs.
C>
C> ### Program History Log:
C> Date | Programmer | Comment
C> -----|------------|--------
C> 1974-09-01 | J. McDonell, J.Howcroft | Initial.
C> 1984-06-27 | Ralph Jones | Change to ibm vs fortran
C> 1989-01-24 | Ralph Jones | Change to microsoft fortran 4.10
C> 1989-03-31 | Ralph Jones | Change to vax-11 fortran
C> 1993-03-16 | D. Shimomura | Renamed from w3ft00() to w3ft07()
C> in order to make minor mods while doing f77. Changes to call sequence;
C> changes to vrbl names; added comments.
C>
C> @param[in] FLDA Real*4 original source grid-point data field
C> @param[in] IA (Input for FLDA)
C> @param[in] JA (Input for FLDA)
C> @param[in] FLDB Real*4 original source grid-point data field
C> @param[in] IB (Input for FLDB)
C> @param[in] JB (Input for FLDB)
C> @param[in] AIPOLE Real*4 common point i-coordinates of the
C> original field, assuming a right-hand cartesian
C> coordinate system. the point need not be inside the bounds of either grid
C> @param[in] AJPOLE Real*4 common point j-coordinates of the
C> original field, assuming a right-hand cartesian
C> coordinate system. the point need not be inside the bounds of either grid
C> and can have fractional values. Common point about which to rotate the gridpoints.
C> @param[in] BIPOLE - Real*4 common point i-coordinates for
C> transformed destination grid
C> @param[in] BJPOLE - Real*4 common point j-coordinates for
C> transformed destination grid
C> @param[in] DSCALE - Real*4 scale-change (dilation) expressed as
C> a ratio of the transformed field to the original field
C> dscale = grdlenkm(destination) / grdlenkm(source)
C> @param[in] ANGLE  - Real*4 degree measure of the angle required to
C> rotate the j-row of the original grid into
C> coincidence with the new grid. (+ counter-
C> clockwise, - clockwise)
C> angle = vertlonw(source) - vertlonw(destination)
C>
C> @param[in] LINEAR - Logical*4 interpolation-method selection switch:
C> - .TRUE. Bi-linear interpolation.
C> - .FALSE. Bi-quadratic interpolation.
C> @param[in] LDEFQQ - Logical*4 default-value switch:
C> if .true. then
C>   use default-value for destination point
C>   out-of-bounds of given grid;
C> else
C>   extrapolate coarsely from nearby bndry point
C> @param[in] DEFALT - Real*4  the default-value to use if ldefqq = .true.
C>
C> @remark List caveats, other helpful hints or information
C> in general 'FLDA' and 'FLDB' cannot be equivalenced
C> although there are situations in which it would be safe to do
C> so. Care should be taken that all of the new grid points lie
C> within the original grid, no error checks are made.
C>
C> @author McDonell & Howcroft @date 1974-09-01
      SUBROUTINE w3ft07(FLDA,IA,JA,AIPOLE,AJPOLE,BIPOLE,BJPOLE,
     A                   DSCALE,ANGLE,LINEAR,LDEFQQ,DEFALT,FLDB,IB,JB)
C
      REAL      FLDA(IA,JA)
      REAL      AIPOLE,AJPOLE
      REAL      BIPOLE,BJPOLE
      REAL      DSCALE
      REAL      ANGLE
      REAL      DEFALT
      REAL      FLDB(IB,JB)
      REAL      ERAS(4)
      REAL      TINY
C
      LOGICAL   LINEAR
      LOGICAL   LDEFQQ
C
      SAVE
C
      DATA  tiny  / 0.001 /
C
C     ... WHERE TINY IS IN UNITS OF 1.0 = 1 GRID INTERVAL
C
C     . . . . .   S T A R T   . . . . . . . . . . . . . . . . . . .
C
      theta = angle * (3.14159/180.)
      sint  = sin(theta)
      cost  = cos(theta)
C
C     ... WE WILL SCAN ALONG THE J-ROW OF THE DESTINATION GRID ...
      DO 288 jn = 1,jb
        brelj  = float(jn) - bjpole
C
        DO 277 in = 1,ib
          breli = float(in) - bipole
          sti = aipole + dscale*(breli*cost - brelj*sint)
          stj = ajpole + dscale*(breli*sint + brelj*cost)
          im = sti
          jm = stj
C
C         ... THE PT(STI,STJ) IS THE LOCATION OF THE FLDB(IN,JN)
C         ... IN FLDA,S COORDINATE SYSTEM
C         ... IS THIS POINT LOCATED OUTSIDE FLDA?
C         ...           ON THE BOUNDARY LINE OF FLDA?
C         ...           ON THE FIRST INTERIOR GRIDPOINT OF FLDA?
C         ...           GOOD INSIDER, AT LEAST 2 INTERIOR GRIDS INSIDE?
          ioff = 0
          joff = 0
          kquad = 0
C
          IF (im .LT. 1) THEN
C           ... LOCATED OUTSIDE OF FLDA, OFF LEFT SIDE ...
            ii = 1
            ioff = 1
          ELSE IF (im .EQ. 1) THEN
C           ... LOCATED ON BOUNDARY OF FLDA, ON LEFT EDGE ...
            kquad = 5
          ELSE
C           ...( IM .GT. 1) ... LOCATED TO RIGHT OF LEFT-EDGE ...
            IF ((ia-im) .LT. 1) THEN
C             ... LOCATED OUTSIDE OF OR EXACTLY ON RIGHT EDGE OF FLDA ..
              ii = ia
              ioff = 1
            ELSE IF ((ia-im) .EQ. 1) THEN
C             ... LOCATED ON FIRST INTERIOR PT WITHIN RIGHT EDGE OF FLDA
              kquad = 5
            ELSE
C             ... (IA-IM) IS .GT. 1) ...GOOD INTERIOR, AT LEAST 2 INSIDE
            ENDIF
          ENDIF
C
C         . . . . . . . . . . . . . . .
C
          IF (jm .LT. 1) THEN
C           ... LOCATED OUTSIDE OF FLDA, OFF BOTTOM ...
            jj = 1
            joff = 1
          ELSE IF (jm .EQ. 1) THEN
C           ... LOCATED ON BOUNDARY OF FLDA, ON BOTTOM EDGE ...
            kquad = 5
          ELSE
C           ...( JM .GT. 1) ... LOCATED ABOVE BOTTOM EDGE ...
            IF ((ja-jm) .LT. 1) THEN
C             ... LOCATED OUTSIDE OF OR EXACTLY ON TOP EDGE OF FLDA ..
              jj = ja
              joff = 1
            ELSE IF ((ja-jm) .EQ. 1) THEN
C             ... LOCATED ON FIRST INTERIOR PT WITHIN TOP EDGE OF FLDA
              kquad = 5
            ELSE
C             ... ((JA-JM) .GT. 1) ...GOOD INTERIOR, AT LEAST 2 INSIDE
            ENDIF
          ENDIF
C
          IF ((ioff + joff) .EQ. 0) THEN
            GO TO 244
          ELSE IF ((ioff + joff) .EQ. 2) THEN
            GO TO 233
          ENDIF
C
          IF (ioff .EQ. 1) THEN
            jj = stj
          ENDIF
          IF (joff .EQ. 1) THEN
            ii = sti
          ENDIF
  233     CONTINUE
          IF (ldefqq) THEN
            fldb(in,jn) = defalt
          ELSE
            fldb(in,jn) = flda(ii,jj)
          ENDIF
          GO TO 277
C
C         . . . . . . . . . . . . .
C
  244     CONTINUE
          i = sti
          j = stj
          xdeli = sti - float(i)
          xdelj = stj - float(j)
C
          IF ((abs(xdeli) .LT. tiny) .AND. (abs(xdelj) .LT. tiny)) THEN
C           ... THIS POINT IS RIGHT AT A GRIDPOINT. NO INTERP NECESSARY
            fldb(in,jn) = flda(i,j)
            GO TO 277
          ENDIF
C
          IF ((kquad .EQ. 5) .OR. (linear)) THEN
C           ... PERFORM BI-LINEAR INTERP ...
            eras(1) = flda(i,j)
            eras(4) = flda(i,j+1)
            eras(2) = eras(1) + xdeli*(flda(i+1,j) - eras(1))
            eras(3) = eras(4) + xdeli*(flda(i+1,j+1) - eras(4))
            di = eras(2) + xdelj*(eras(3) - eras(2))
            GO TO 266
C
          ELSE
C           ... PERFORM BI-QUADRATIC INTERP ...
            xi2tm = xdeli * (xdeli-1.) * 0.25
            xj2tm = xdelj * (xdelj-1.) * 0.25
            j1 = j - 1
            DO 255 k=1,4
              eras(k)=(flda(i+1,j1)-flda(i,j1))*xdeli+flda(i,j1)+
     a        (flda(i-1,j1)-flda(i,j1)-flda(i+1,j1)+flda(i+2,j1))*xi2tm
              j1 = j1 + 1
  255       CONTINUE
C
            di = eras(2) +  xdelj*(eras(3)-eras(2)) +
     a                      xj2tm*(eras(4)-eras(3)-eras(2)+eras(1))
            GO TO 266
          ENDIF
C
  266     CONTINUE
          fldb(in,jn) = di
  277   CONTINUE
  288 CONTINUE
C
      RETURN
      END