interp_ts_bs_splin Subroutine

   subroutine interp_ts_bs_splin(ms_fe_f)
   implicit none
   type(MS_FE_FILM), intent(inout) :: ms_fe_f

      integer(kind=I4) :: i, j, inbvx, inbvy, iloy, iflag, ts_nx, ts_ny, e, ed, bs_nx, bs_ny, n_edg, n_nod, nn, node_number
      real(kind=R8) :: x_n, y_n, val
      real(kind=R8), dimension(:,:), allocatable :: log_press

      real(kind=R8), dimension(:,:), allocatable :: coeff
      real(kind=R8), dimension(:),   allocatable :: tx ! x knots
      real(kind=R8), dimension(:),   allocatable :: ty ! y knots
      real(kind=R8), dimension(:),   allocatable ::  x
      real(kind=R8), dimension(:),   allocatable ::  y

      integer(kind=I4), parameter :: deg = 2
      logical(kind=I4) :: mixture

      mixture = (ms_fe_f%ts_fe_f%data_f%fl%fluid_type == MIXT)

      ts_nx = ms_fe_f%ts_fe_f%m%nx
      ts_ny = ms_fe_f%ts_fe_f%m%ny

      allocate( coeff(1:ts_nx, 1:ts_ny) )
      allocate( tx(1:(ts_nx +deg)),  &
                ty(1:(ts_ny +deg)) )
      allocate(  x(1: ts_nx    ),  &
                 y(1: ts_ny    ) )

      allocate( log_press(1:ts_nx, 1:ts_ny) )
      log_press = reshape( ms_fe_f%ts_fe_f%vn(:, P_N), (/ts_nx, ts_ny/) )

      if (mixture) log_press = log(log_press)

      x(1:ts_nx) = ms_fe_f%ts_fe_f%m%x(1:ts_nx)
      i = 0
      do j = 1, ts_nx*(ts_ny -1) +1, ts_nx
         i = i +1
         y(i) = ms_fe_f%ts_fe_f%m%y(j)
      enddo

      iflag = 0
      call db2ink(   x = x(1:ts_nx),                    & ! Array of x abcissae. Must be strictly increasing.
                    nx = ts_nx,                         & ! Number of x abcissae
                     y = y(1:ts_ny),                    & ! Array of y abcissae. Must be strictly increasing.
                    ny = ts_ny,                         & ! Number of y abcissae
                   fcn = log_press(1:ts_nx, 1:ts_ny),   & ! Array of function values to interpolate. fcn(i,j) should
                                                          !    contain the function value at the point (x(i),y(j))
                    kx = deg,                           & ! The order of spline pieces in x (>= 2, < nx). (order = polynomial degree + 1)
                    ky = deg,                           & ! The order of spline pieces in y (>= 2, < ny). (order = polynomial degree + 1)
                    tx = tx(1:(ts_nx +deg)),            & ! The knots in the x direction for the spline interpolant.
                                                          !    If iflag=0 these are chosen by [[db2ink]].
                                                          !    If iflag=1 these are specified by the user.
                                                          !    Must be non-decreasing.
                    ty = ty(1:(ts_ny +deg)),            & ! The knots in the y direction for the spline interpolant.
                                                          !    If iflag=0 these are chosen by [[db2ink]].
                                                          !    If iflag=1 these are specified by the user.
                                                          !    Must be non-decreasing.
                 bcoef = coeff(1:ts_nx, 1:ts_ny),       & ! Array of coefficients of the b-spline interpolant.
                 iflag = iflag)                           ! **on input:**  0 = knot sequence chosen by [[db2ink]].
                                                          !                1 = knot sequence chosen by user.
                                                          ! **on output:** 1 = successful execution.
                                                          !                2 = iflag out of range.
                                                          !                3 = nx out of range.
                                                          !                4 = kx out of range.
                                                          !                5 = x not strictly increasing.
                                                          !                6 = tx not non-decreasing.
                                                          !                7 = ny out of range.
                                                          !                8 = ky out of range.
                                                          !                9 = y not strictly increasing.
                                                          !               10 = ty not non-decreasing.
      if (iflag/=1) error stop 'error calling db2ink'

      inbvx = 1
      inbvy = 1
      iloy  = 1
      do e = 1, ms_fe_f%ts_fe_f%m%ne
         bs_nx = ms_fe_f%bs_fe_f(e)%m%nx
         bs_ny = ms_fe_f%bs_fe_f(e)%m%ny
         n_edg = ms_fe_f%bs_fe_f(e)%m%ned
         do ed = 1, n_edg
            n_nod = ms_fe_f%bs_fe_f(e)%m%ed(ed)%n
            do nn = 1, n_nod
               node_number = ms_fe_f%bs_fe_f(e)%m%ed(ed)%nm(nn)
               x_n = ms_fe_f%bs_fe_f(e)%m%x(node_number)
               y_n = ms_fe_f%bs_fe_f(e)%m%y(node_number)

               call db2val(xval = x_n,                      &  ! xval     !! x coordinate of evaluation point.
                           yval = y_n,                      &  ! yval     !! y coordinate of evaluation point.
                            idx = 0,                        &  ! idx      !! x derivative of piecewise polynomial to evaluate.
                            idy = 0,                        &  ! idy      !! y derivative of piecewise polynomial to evaluate.
                             tx = tx(1:(ts_nx +deg)),       &  ! tx       !! sequence of knots defining the piecewise polynomial in the x direction. (same as in last call to [[db2ink]])
                             ty = ty(1:(ts_ny +deg)),       &  ! ty       !! sequence of knots defining the piecewise polynomial in the y direction. (same as in last call to [[db2ink]])
                             nx = ts_nx,                    &  ! nx       !! the number of interpolation points in x. (same as in last call to [[db2ink]])
                             ny = ts_ny,                    &  ! ny       !! the number of interpolation points in y. (same as in last call to [[db2ink]])
                             kx = deg,                      &  ! kx       !! order of polynomial pieces in x. (same as in last call to [[db2ink]])
                             ky = deg,                      &  ! ky       !! order of polynomial pieces in y. (same as in last call to [[db2ink]])
                          bcoef = coeff(1:ts_nx, 1:ts_ny),  &  ! bcoef    !! the b-spline coefficients computed by [[db2ink]].
                              f = val,                      &  ! f        !! interpolated value &
                          iflag = iflag,                    &  ! iflag    !! status flag: 0 : no errors, /=0 : error
                          inbvx = inbvx,                    &  ! inbvx    !! initialization parameter which must be set to 1 the first time this routine is called, and must not be changed by the user.
                          inbvy = inbvy,                    &  ! inbvy    !! initialization parameter which must be set to 1 the first time this routine is called, and must not be changed by the user.
                           iloy = iloy)                        ! iloy     !! initialization parameter which must be set to 1 the first time this routine is called, and must not be changed by the user.
               if (iflag/=0) error stop 'error'
               if (mixture) val = exp(val)

               ms_fe_f%bs_fe_f(e)%vn(node_number, P_N) = val
               ms_fe_f%bs_fe_f(e)%bc(node_number, REY) = 0
            enddo
         enddo
      enddo

      deallocate( log_press, coeff, tx, ty, x, y )

   return
   endsubroutine interp_ts_bs_splin