Replacing calls to E02DAF

E02DAF can be replaced with the following code.

*   Declare variables
      INTEGER ID                      ! Specifies use value not differential
      INTEGER IFAIL                   ! Was the surface successfully created?
      INTEGER MXY                     ! Size of the grid
      INTEGER ORD                     ! Order of polynomial used

      DOUBLE PRECISION DVALUE         ! Interpolated value returned 
      DOUBLE PRECISION FV1(8,8)       ! Grid Z values
      DOUBLE PRECISION X1(8), Y1(8)   ! X,Y grid locations
      DOUBLE PRECISION XD, YD         ! X,Y coord for interpolation
 
      DOUBLE PRECISION BCOEF(8,8)     ! Array used by PDA_DB2INK
      DOUBLE PRECISION TX(11),TY(11)  ! Array used by PDA_DB2INK
      DOUBLE PRECISION WORK(168)      ! Array used by PDA_DB2INK
 
*   Subroutine initial error values. 
      IFAIL=0
      STATUS=0

*   Order of polynomial used in this example.
      ORD=3
 
*   Use value rather than differential.
      ID=0

*   Size of data grid. 8x8 in this instance.
      MXY=8
     
*   Build the surface fit using the grid contents.
      CALL PDA_DB2INK(X1,MXY,Y1,MXY,FV1,MXY,
     :                ORD,ORD,TX,TY,BCOEF,
     :                WORK,IFAIL,STATUS)
 

*   If IFAIL=1 then okay to interpolate values.
      IF (IFAIL.EQ.1) THEN
   
*       Setup evaluation routine.
          IFAIL=0
          CALL PDA_DB2VAL(XD,YD,ID,ID,TX,TY,
     :                    MXY,MXY,ORD,ORD,BCOEF,WORK,
     :                    DVALUE,IFAIL,STATUS)
 
      END IF

In the example shown the value of the surface fitted is calculated at a point defined by XD and YD. The surface constructed is a 3rd order, polynomial and the value returned by PDA_DB2VAL is the surface value not its differential. The example assumes that the array FV1() already contains values for the data at each of the grid points defined in arrays X1 and Y1. The size of the work arrays is determined by the size of the grid required. The variable DVALUE contains the interpolated value.