fitp_surf1 Subroutine

  pure subroutine fitp_surf1 (m,n,x,y,z,iz,zx1,zxm,zy1,zyn,zxy11, &
       zxym1,zxy1n,zxymn,islpsw,zp,temp,sigma,ierr)
    implicit none
    integer, intent(in) :: m, n, iz, islpsw
    real, dimension(m), intent(in) :: x, zy1, zyn
    real, dimension(n), intent(in) :: y, zx1, zxm
    real, dimension(iz,n), intent(in) :: z
    real, intent(in) :: zxy11, zxym1, zxy1n, zxymn, sigma
    real, dimension(m,n,3), intent(out) :: zp
    real, dimension(n+n+m), intent(out) :: temp
    integer, intent(out) :: ierr
!
!                                 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 subroutine determines the parameters necessary to
! compute an interpolatory surface passing through a rect-
! angular grid of functional values. the surface determined
! can be represented as the tensor product of splines under
! tension. the x- and y-partial derivatives around the
! boundary and the x-y-partial derivatives at the four
! corners may be specified or omitted. for actual mapping
! of points onto the surface it is necessary to call the
! function surf2.
!
! on input--
!
!   m is the number of grid lines in the x-direction, i. e.
!   lines parallel to the y-axis (m .ge. 2).
!
!   n is the number of grid lines in the y-direction, i. e.
!   lines parallel to the x-axis (n .ge. 2).
!
!   x is an array of the m x-coordinates of the grid lines
!   in the x-direction. these should be strictly increasing.
!
!   y is an array of the n y-coordinates of the grid lines
!   in the y-direction. these should be strictly increasing.
!
!   z is an array of the m * n functional values at the grid
!   points, i. e. z(i,j) contains the functional value at
!   (x(i),y(j)) for i = 1,...,m and j = 1,...,n.
!
!   iz is the row dimension of the matrix z used in the
!   calling program (iz .ge. m).
!
!   zx1 and zxm are arrays of the m x-partial derivatives
!   of the function along the x(1) and x(m) grid lines,
!   respectively. thus zx1(j) and zxm(j) contain the x-part-
!   ial derivatives at the points (x(1),y(j)) and
!   (x(m),y(j)), respectively, for j = 1,...,n. either of
!   these parameters will be ignored (and approximations
!   supplied internally) if islpsw so indicates.
!
!   zy1 and zyn are arrays of the n y-partial derivatives
!   of the function along the y(1) and y(n) grid lines,
!   respectively. thus zy1(i) and zyn(i) contain the y-part-
!   ial derivatives at the points (x(i),y(1)) and
!   (x(i),y(n)), respectively, for i = 1,...,m. either of
!   these parameters will be ignored (and estimations
!   supplied internally) if islpsw so indicates.
!
!   zxy11, zxym1, zxy1n, and zxymn are the x-y-partial
!   derivatives of the function at the four corners,
!   (x(1),y(1)), (x(m),y(1)), (x(1),y(n)), and (x(m),y(n)),
!   respectively. any of the parameters will be ignored (and
!   estimations supplied internally) if islpsw so indicates.
!
!   islpsw contains a switch indicating which boundary
!   derivative information is user-supplied and which
!   should be estimated by this subroutine. to determine
!   islpsw, let
!        i1 = 0 if zx1 is user-supplied (and = 1 otherwise),
!        i2 = 0 if zxm is user-supplied (and = 1 otherwise),
!        i3 = 0 if zy1 is user-supplied (and = 1 otherwise),
!        i4 = 0 if zyn is user-supplied (and = 1 otherwise),
!        i5 = 0 if zxy11 is user-supplied
!                                       (and = 1 otherwise),
!        i6 = 0 if zxym1 is user-supplied
!                                       (and = 1 otherwise),
!        i7 = 0 if zxy1n is user-supplied
!                                       (and = 1 otherwise),
!        i8 = 0 if zxymn is user-supplied
!                                       (and = 1 otherwise),
!   then islpsw = i1 + 2*i2 + 4*i3 + 8*i4 + 16*i5 + 32*i6
!                   + 64*i7 + 128*i8
!   thus islpsw = 0 indicates all derivative information is
!   user-supplied and islpsw = 255 indicates no derivative
!   information is user-supplied. any value between these
!   limits is valid.
!
!   zp is an array of at least 3*m*n locations.
!
!   temp is an array of at least n+n+m locations which is
!   used for scratch storage.
!
! and
!
!   sigma contains the tension factor. this value indicates
!   the curviness desired. if abs(sigma) is nearly zero
!   (e. g. .001) the resulting surface is approximately the
!   tensor product of cubic splines. if abs(sigma) is large
!   (e. g. 50.) the resulting surface is approximately
!   bi-linear. if sigma equals zero tensor products of
!   cubic splines result. a standard value for sigma is
!   approximately 1. in absolute value.
!
! on output--
!
!   zp contains the values of the xx-, yy-, and xxyy-partial
!   derivatives of the surface at the given nodes.
!
!   ierr contains an error flag,
!        = 0 for normal return,
!        = 1 if n is less than 2 or m is less than 2,
!        = 2 if the x-values or y-values are not strictly
!            increasing.
!
! and
!
!   m, n, x, y, z, iz, zx1, zxm, zy1, zyn, zxy11, zxym1,
!   zxy1n, zxymn, islpsw, and sigma are unaltered.
!
! this subroutine references package modules ceez, terms,
! and snhcsh.
!
!-----------------------------------------------------------
    integer :: jbak, jbakp1, jm1, jp1, im1, ibakp1
    integer :: npibak, ip1, ibak, nm1, np1, mm1, mp1, npm, j
    integer :: npmpj, i, npi
    real :: diagi, sdiag1, del2, zxymns, delxmm, del1
    real :: diag1, deli, diag2, diagin, sdiag2, t
    real :: delxm, sigmay, dely1, c1, dely2, c2
    real :: delx1, delx2, zxy1ns, c3, delyn, sigmax, delynm

    mm1 = m-1
    mp1 = m+1
    nm1 = n-1
    np1 = n+1
    npm = n+m
    ierr = 0
    if (n .le. 1 .or. m .le. 1) go to 46
    if (y(n) .le. y(1)) go to 47
!
! denormalize tension factor in y-direction
!
    sigmay = abs(sigma) * (n - 1) / (y(n) - y(1))
!
! obtain y-partial derivatives along y = y(1)
!
    if ((islpsw/8)*2 .ne. (islpsw/4)) go to 2
    do i = 1,m
       zp(i,1,1) = zy1(i)
    end do
    go to 5
2   dely1 = y(2)-y(1)
    dely2 = dely1+dely1
    if (n .gt. 2) dely2 = y(3)-y(1)
    if (dely1 .le. 0. .or. dely2 .le. dely1) go to 47
    call fitp_ceez (dely1,dely2,sigmay,c1,c2,c3,n)
    do i = 1,m
       zp(i,1,1) = c1*z(i,1)+c2*z(i,2)
    end do
    if (n .eq. 2) go to 5
    do i = 1,m
       zp(i,1,1) = zp(i,1,1)+c3*z(i,3)
    end do
!
! obtain y-partial derivatives along y = y(n)
!
5   if ((islpsw/16)*2 .ne. (islpsw/8)) go to 7
    do i = 1,m
       npi = n+i
       temp(npi) = zyn(i)
    end do
    go to 10
7   delyn = y(n)-y(nm1)
    delynm = delyn+delyn
    if (n .gt. 2) delynm = y(n)-y(n-2)
    if (delyn .le. 0. .or. delynm .le. delyn) go to 47
    call fitp_ceez (-delyn,-delynm,sigmay,c1,c2,c3,n)
    do i = 1,m
       npi = n+i
       temp(npi) = c1*z(i,n)+c2*z(i,nm1)
    end do
    if (n .eq. 2) go to 10
    do i = 1,m
       npi = n+i
       temp(npi) = temp(npi)+c3*z(i,n-2)
    end do
10  if (x(m) .le. x(1)) go to 47
!
! denormalize tension factor in x-direction
!
    sigmax = abs(sigma) * (m - 1) / (x(m) - x(1))
!
! obtain x-partial derivatives along x = x(1)
!
    if ((islpsw/2)*2 .ne. islpsw) go to 12
    do j = 1,n
       zp(1,j,2) = zx1(j)
    end do
    if ((islpsw/32)*2 .eq. (islpsw/16) .and. (islpsw/128)*2  .eq. (islpsw/64)) go to 15
12  delx1 = x(2)-x(1)
    delx2 = delx1+delx1
    if (m .gt. 2) delx2 = x(3)-x(1)
    if (delx1 .le. 0. .or. delx2 .le. delx1) go to 47
    call fitp_ceez (delx1,delx2,sigmax,c1,c2,c3,m)
    if ((islpsw/2)*2 .eq. islpsw) go to 15
    do j = 1,n
       zp(1,j,2) = c1*z(1,j)+c2*z(2,j)
    end do
    if (m .eq. 2) go to 15
    do j = 1,n
       zp(1,j,2) = zp(1,j,2)+c3*z(3,j)
    end do
!
! obtain x-y-partial derivative at (x(1),y(1))
!
15  if ((islpsw/32)*2 .ne. (islpsw/16)) go to 16
    zp(1,1,3) = zxy11
    go to 17
16  zp(1,1,3) = c1*zp(1,1,1)+c2*zp(2,1,1)
    if (m .gt. 2) zp(1,1,3) = zp(1,1,3)+c3*zp(3,1,1)
!
! obtain x-y-partial derivative at (x(1),y(n))
!
17  if ((islpsw/128)*2 .ne. (islpsw/64)) go to 18
    zxy1ns = zxy1n
    go to 19
18  zxy1ns = c1*temp(n+1)+c2*temp(n+2)
    if (m .gt. 2) zxy1ns = zxy1ns+c3*temp(n+3)
!
! obtain x-partial derivative along x = x(m)
!
19  if ((islpsw/4)*2 .ne. (islpsw/2)) go to 21
    do j = 1,n
       npmpj = npm+j
       temp(npmpj) = zxm(j)
    end do
    if ((islpsw/64)*2 .eq. (islpsw/32) .and.(islpsw/256)*2 .eq. (islpsw/128)) go to 24
21  delxm = x(m)-x(mm1)
    delxmm = delxm+delxm
    if (m .gt. 2) delxmm = x(m)-x(m-2)
    if (delxm .le. 0. .or. delxmm .le. delxm) go to 47
    call fitp_ceez (-delxm,-delxmm,sigmax,c1,c2,c3,m)
    if ((islpsw/4)*2 .eq. (islpsw/2)) go to 24
    do j = 1,n
       npmpj = npm+j
       temp(npmpj) = c1*z(m,j)+c2*z(mm1,j)
    end do
    if (m .eq. 2) go to 24
    do j = 1,n
       npmpj = npm+j
       temp(npmpj) = temp(npmpj)+c3*z(m-2,j)
    end do
!
! obtain x-y-partial derivative at (x(m),y(1))
!
24  if ((islpsw/64)*2 .ne. (islpsw/32)) go to 25
    zp(m,1,3) = zxym1
    go to 26
25  zp(m,1,3) = c1*zp(m,1,1)+c2*zp(mm1,1,1)
    if (m .gt. 2) zp(m,1,3) = zp(m,1,3)+c3*zp(m-2,1,1)
!
! obtain x-y-partial derivative at (x(m),y(n))
!
26  if ((islpsw/256)*2 .ne. (islpsw/128)) go to 27
    zxymns = zxymn
    go to 28
27  zxymns = c1*temp(npm)+c2*temp(npm-1)
    if (m .gt. 2) zxymns = zxymns+c3*temp(npm-2)
!
! set up right hand sides and tridiagonal system for y-grid
! perform forward elimination
!
28  del1 = y(2)-y(1)
    if (del1 .le. 0.) go to 47
    deli = 1./del1
    do i = 1,m
       zp(i,2,1) = deli*(z(i,2)-z(i,1))
    end do
    zp(1,2,3) = deli*(zp(1,2,2)-zp(1,1,2))
    zp(m,2,3) = deli*(temp(npm+2)-temp(npm+1))
    call fitp_terms (diag1,sdiag1,sigmay,del1)
    diagi = 1./diag1
    do i = 1,m
       zp(i,1,1) = diagi*(zp(i,2,1)-zp(i,1,1))
    end do
    zp(1,1,3) = diagi*(zp(1,2,3)-zp(1,1,3))
    zp(m,1,3) = diagi*(zp(m,2,3)-zp(m,1,3))
    temp(1) = diagi*sdiag1
    if (n .eq. 2) go to 34
    do j = 2,nm1
       jm1 = j-1
       jp1 = j+1
       npmpj = npm+j
       del2 = y(jp1)-y(j)
       if (del2 .le. 0.) go to 47
       deli = 1./del2
       do i = 1,m
          zp(i,jp1,1) = deli*(z(i,jp1)-z(i,j))
       end do
       zp(1,jp1,3) = deli*(zp(1,jp1,2)-zp(1,j,2))
       zp(m,jp1,3) = deli*(temp(npmpj+1)-temp(npmpj))
       call fitp_terms (diag2,sdiag2,sigmay,del2)
       diagin = 1./(diag1+diag2-sdiag1*temp(jm1))
       do i = 1,m
          zp(i,j,1) = diagin*(zp(i,jp1,1)-zp(i,j,1)-sdiag1*zp(i,jm1,1))
       end do
       zp(1,j,3) = diagin*(zp(1,jp1,3)-zp(1,j,3)-sdiag1*zp(1,jm1,3))
       zp(m,j,3) = diagin*(zp(m,jp1,3)-zp(m,j,3)-sdiag1*zp(m,jm1,3))
       temp(j) = diagin*sdiag2
       diag1 = diag2
       sdiag1 = sdiag2
    end do
34  diagin = 1./(diag1-sdiag1*temp(nm1))
    do i = 1,m
       npi = n+i
       zp(i,n,1) = diagin*(temp(npi)-zp(i,n,1)-sdiag1*zp(i,nm1,1))
    end do
    zp(1,n,3) = diagin*(zxy1ns-zp(1,n,3)-sdiag1*zp(1,nm1,3))
    temp(n) = diagin*(zxymns-zp(m,n,3)-sdiag1*zp(m,nm1,3))
!
! perform back substitution
!
    do j = 2,n
       jbak = np1-j
       jbakp1 = jbak+1
       t = temp(jbak)
       do i = 1,m
          zp(i,jbak,1) = zp(i,jbak,1)-t*zp(i,jbakp1,1)
       end do
       zp(1,jbak,3) = zp(1,jbak,3)-t*zp(1,jbakp1,3)
       temp(jbak) = zp(m,jbak,3)-t*temp(jbakp1)
    end do
!
! set up right hand sides and tridiagonal system for x-grid
! perform forward elimination
!
    del1 = x(2)-x(1)
    if (del1 .le. 0.) go to 47
    deli = 1./del1
    do j = 1,n
       zp(2,j,2) = deli*(z(2,j)-z(1,j))
       zp(2,j,3) = deli*(zp(2,j,1)-zp(1,j,1))
    end do
    call fitp_terms (diag1,sdiag1,sigmax,del1)
    diagi = 1./diag1
    do j = 1,n
       zp(1,j,2) = diagi*(zp(2,j,2)-zp(1,j,2))
       zp(1,j,3) = diagi*(zp(2,j,3)-zp(1,j,3))
    end do
    temp(n+1) = diagi*sdiag1
    if (m  .eq. 2) go to 43
    do i = 2,mm1
       im1 = i-1
       ip1 = i+1
       npi = n+i
       del2 = x(ip1)-x(i)
       if (del2 .le. 0.) go to 47
       deli = 1./del2
       do j = 1,n
          zp(ip1,j,2) = deli*(z(ip1,j)-z(i,j))
          zp(ip1,j,3) = deli*(zp(ip1,j,1)-zp(i,j,1))
       end do
       call fitp_terms (diag2,sdiag2,sigmax,del2)
       diagin = 1./(diag1+diag2-sdiag1*temp(npi-1))
       do j = 1,n
          zp(i,j,2) = diagin*(zp(ip1,j,2)-zp(i,j,2)-sdiag1*zp(im1,j,2))
          zp(i,j,3) = diagin*(zp(ip1,j,3)-zp(i,j,3)-sdiag1*zp(im1,j,3))
       end do
       temp(npi) = diagin*sdiag2
       diag1 = diag2
       sdiag1 = sdiag2
    end do
43  diagin = 1./(diag1-sdiag1*temp(npm-1))
    do j = 1,n
       npmpj = npm+j
       zp(m,j,2) = diagin*(temp(npmpj)-zp(m,j,2)-sdiag1*zp(mm1,j,2))
       zp(m,j,3) = diagin*(temp(j)-zp(m,j,3)-sdiag1*zp(mm1,j,3))
    end do
!
! perform back substitution
!
    do i = 2,m
       ibak = mp1-i
       ibakp1 = ibak+1
       npibak = n+ibak
       t = temp(npibak)
       do j = 1,n
          zp(ibak,j,2) = zp(ibak,j,2)-t*zp(ibakp1,j,2)
          zp(ibak,j,3) = zp(ibak,j,3)-t*zp(ibakp1,j,3)
       end do
    end do
    return
!
! too few points
!
46  ierr = 1
    return
!
! points not strictly increasing
!
47  ierr = 2
    return
  end subroutine fitp_surf1