    PROGRAM nag_spline_2d_ex01

      ! Example Program Text for nag_spline_2d
      ! NAG fl90, Release 4. NAG Copyright 2000.

      ! .. Use Statements ..
      USE nag_examples_io, ONLY : nag_std_in, nag_std_out
      USE nag_spline_2d, ONLY : nag_spline_2d_comm_wp => nag_spline_2d_comm_dp &
       , nag_spline_2d_auto_fit, nag_spline_2d_eval, nag_deallocate
      ! .. Implicit None Statement ..
      IMPLICIT NONE
      ! .. Intrinsic Functions ..
      INTRINSIC KIND, REAL
      ! .. Parameters ..
      INTEGER, PARAMETER :: wp = KIND(1.0D0)
      ! .. Local Scalars ..
      INTEGER :: i, j, mx, my, nx, ny, p, q
      REAL (wp) :: dx, dy, smooth, theta, xmax, xmin, ymax, ymin
      CHARACTER (1) :: start
      TYPE (nag_spline_2d_comm_wp) :: spline
      ! .. Local Arrays ..
      REAL (wp), ALLOCATABLE :: f(:,:), s(:,:), u(:), v(:), x(:), y(:)
      ! .. Executable Statements ..
      WRITE (nag_std_out,*) 'Example Program Results for nag_spline_2d_ex01'

      READ (nag_std_in,*)          ! Skip heading in data file
      READ (nag_std_in,*) mx, my
      READ (nag_std_in,*) nx, xmin, xmax
      READ (nag_std_in,*) ny, ymin, ymax

      ALLOCATE (x(mx),y(my),f(mx,my),u(nx),v(ny),s(nx,ny)) ! Allocate storage

      READ (nag_std_in,*) x
      READ (nag_std_in,*) y
      READ (nag_std_in,*) f

      start = 'cold start'

      DO

        ! Read in successive values of smooth and generate spline for
        ! each.
        READ (nag_std_in,*,end=20) smooth

        ! Determine the spline approximation.

        CALL nag_spline_2d_auto_fit(start,x,y,f,smooth,spline,p=p,q=q, &
         theta=theta)

        WRITE (nag_std_out,'(//1X,A,1P,E13.4)') &
         'Calling with smoothing factor =', smooth
        WRITE (nag_std_out,'(1X,A,1P,E13.4)') 'Sum of squared residuals =', &
         theta
        WRITE (nag_std_out,'(1X,A,I8)') &
         'Total number of knots in x-direction =', p
        WRITE (nag_std_out,'(1X,A,I8)') &
         'Total number of knots in y-direction =', q

        ! Evaluate the spline on a rectangular grid at nx*ny points
        ! over the domain [xmin, xmax] x [ymin, ymax].
        dx = (xmax-xmin)/REAL(nx-1,kind=wp)
        dy = (ymax-ymin)/REAL(ny-1,kind=wp)
        u = (/ (xmin+(i-1)*dx,i=1,nx-1), xmax/)
        v = (/ (ymin+(j-1)*dy,j=1,ny-1), ymax/)

        CALL nag_spline_2d_eval(spline,u,v,s)

        WRITE (nag_std_out,'(/1X,A)') 'Values of computed spline:'
        WRITE (nag_std_out,'(/11X,A,7F8.2)') 'x', u
        WRITE (nag_std_out,*) '     y'
        DO j = ny, 1, -1
          WRITE (nag_std_out,'(1X,F8.2,3X,7F8.2)') v(j), s(:,j)
        END DO

        start = 'warm start'

      END DO
20    CONTINUE

      DEALLOCATE (x,y,f,u,v,s)     ! Deallocate storage

      CALL nag_deallocate(spline)  ! Free structure allocated by NAG fl90

    END PROGRAM nag_spline_2d_ex01
