fitp_surf2 Function

  real function fitp_surf2 (xx,yy,m,n,x,y,z,iz,zp,sigma)
    implicit none
    
    real, intent(in) :: xx, yy, sigma
    integer, intent(in) :: m, n, iz
    real, dimension(m), intent(in) :: x
    real, dimension(n), intent(in) :: y
    real, dimension(iz,n), intent(in) :: z
    real, dimension(m,n,3), intent(in) :: zp
!
!                                 coded by alan kaylor cline
!                           from fitpack -- january 26, 1987
!                        a curve and surface fitting package
!                      a product of pleasant valley software
!                  8603 altus cove, austin, texas 78759, usa
!
! this function interpolates a surface at a given coordinate
! pair using a bi-spline under tension. the subroutine surf1
! should be called earlier to determine certain necessary
! parameters.
!
! on input--
!
!   xx and yy contain the x- and y-coordinates of the point
!   to be mapped onto the interpolating surface.
!
!   m and n contain the number of grid lines in the x- and
!   y-directions, respectively, of the rectangular grid
!   which specified the surface.
!
!   x and y are arrays containing the x- and y-grid values,
!   respectively, each in increasing order.
!
!   z is a matrix containing the m * n functional values
!   corresponding to the grid values (i. e. z(i,j) is the
!   surface value at the point (x(i),y(j)) for i = 1,...,m
!   and j = 1,...,n).
!
!   iz contains the row dimension of the array z as declared
!   in the calling program.
!
!   zp is an array of 3*m*n locations stored with the
!   various surface derivative information determined by
!   surf1.
!
! and
!
!   sigma contains the tension factor (its sign is ignored).
!
! the parameters m, n, x, y, z, iz, zp, and sigma should be
! input unaltered from the output of surf1.
!
! on output--
!
!   surf2 contains the interpolated surface value.
!
! none of the input parameters are altered.
!
! this function references package modules intrvl and
! snhcsh.
!
!-----------------------------------------------------------
    integer :: im1, i, j, jm1
    real :: zxxi, zxxim1, zim1, zi
    real :: sigmax, del1, del2
    real :: sinhms, sinhm2, sinhm1, dels, sigmay
!
! inline one dimensional cubic spline interpolation
!
!/Moved to contained function
!    hermz (f1,f2,fp1,fp2) = (f2*del1+f1*del2)/dels-del1* &
!         del2*(fp2*(del1+dels)+ &
!         fp1*(del2+dels))/(6.*dels)
!
! inline one dimensional spline under tension interpolation
!
!/Moved to contained function
!    hermnz (f1,f2,fp1,fp2,sigmap) = (f2*del1+f1*del2)/dels &
!         +(fp2*del1*(sinhm1-sinhms) &
!         +fp1*del2*(sinhm2-sinhms) &
!         )/(sigmap*sigmap*dels*(1.+sinhms))

!
! denormalize tension factor in x and y direction
!
    sigmax = abs(sigma) * (m - 1) / (x(m) - x(1))
    sigmay = abs(sigma) * (n - 1) / (y(n) - y(1))
!
! determine y interval
!
    jm1 = fitp_intrvl (yy,y,n)
    j = jm1+1
!
! determine x interval
!
    im1 = fitp_intrvl (xx,x,m)
    i = im1+1
    del1 = yy-y(jm1)
    del2 = y(j)-yy
    dels = y(j)-y(jm1)
    if (is_zero(sigmay)) then
       ! perform four interpolations in y-direction
       zim1 = hermz(z(i-1,j-1),z(i-1,j),zp(i-1,j-1,1),zp(i-1,j,1))
       zi = hermz(z(i,j-1),z(i,j),zp(i,j-1,1),zp(i,j,1))
       zxxim1 = hermz(zp(i-1,j-1,2),zp(i-1,j,2),zp(i-1,j-1,3),zp(i-1,j,3))
       zxxi = hermz(zp(i,j-1,2),zp(i,j,2),zp(i,j-1,3),zp(i,j,3))
    else
       sinhm1 = sinhm_fun(sigmay * del1)
       sinhm2 = sinhm_fun(sigmay * del2)
       sinhms = sinhm_fun(sigmay * dels)
       zim1 = hermnz(z(i-1,j-1),z(i-1,j),zp(i-1,j-1,1),zp(i-1,j,1),sigmay)
       zi = hermnz(z(i,j-1),z(i,j),zp(i,j-1,1),zp(i,j,1),sigmay)
       zxxim1 = hermnz(zp(i-1,j-1,2),zp(i-1,j,2),zp(i-1,j-1,3),zp(i-1,j,3),sigmay)
       zxxi = hermnz(zp(i,j-1,2),zp(i,j,2),zp(i,j-1,3),zp(i,j,3),sigmay)
    end if

    ! perform final interpolation in x-direction

    del1 = xx-x(im1)
    del2 = x(i)-xx
    dels = x(i)-x(im1)

    if (is_zero(sigmax)) then
       fitp_surf2 = hermz(zim1,zi,zxxim1,zxxi)
    else
       sinhm1 = sinhm_fun(sigmax * del1)
       sinhm2 = sinhm_fun(sigmax * del2)
       sinhms = sinhm_fun(sigmax * dels)
       fitp_surf2 = hermnz(zim1,zi,zxxim1,zxxi,sigmax)
    end if
  contains
    !> FIXME : Add documentation
    elemental real function hermz(f1,f2,fp1,fp2)
      !Note del1,sinhm1 etc. are known because we are in subroutine scope
      implicit none
      real, intent(in) :: f1,f2,fp1,fp2
      hermz= (f2*del1+f1*del2)/dels-del1* &
           del2*(fp2*(del1+dels)+ &
           fp1*(del2+dels))/(6.*dels)
    end function hermz

    !> FIXME : Add documentation
    elemental real function hermnz(f1,f2,fp1,fp2,sigmap)
      !Note del1,sinhm1 etc. are known because we are in subroutine scope
      implicit none
      real, intent(in) :: f1,f2,fp1,fp2,sigmap
      hermnz= (f2*del1+f1*del2)/dels &
           +(fp2*del1*(sinhm1-sinhms) &
           +fp1*del2*(sinhm2-sinhms) &
           )/(sigmap*sigmap*dels*(1.+sinhms))
    end function hermnz
  end function fitp_surf2