compute_bicubic_interp Subroutine

    subroutine compute_bicubic_interp(x, y, z, dzdx, dzdy, nx, ny, xgrd, ygrd, &
                                      v, vx, vy, vxy)
        !! Evaluate the z value at position (x, y) of a surface built as a 
        !! bicubic spline interpolating the points     
        !! (xgrd(\(x_i\), \(y_i\)), ygrd(\(x_i\), \(y_i\)),
        !! vxgrd(\(x_i\), \(y_i\))).     
        !! In order to do so, also the derivatives of the surface at the points
        !! of the grid are needed; in particular if we consider the surface
        !! points as
        !! \[(x_i, y_i, V(x_i, y_i))\]
        !! also value of \(\frac{\partial V}{\partial x}\),
        !! \(\frac{\partial V}{\partial y}\) and 
        !! \(\frac{\partial^2 V}{\partial x \partial y}\) at grid points
        !! are needed.
        
        use mod_memory, only: ip, rp

        implicit none

        real(rp), intent(in) :: x, y
        !! Coordinates at which the z value of the surface should be computed
        real(rp), intent(out) :: z
        !! The z value of the surface
        real(rp), intent(out) :: dzdx, dzdy
        !! Derivatives of z valure w.r.t. x and y coordinates
        integer(ip), intent(in) :: nx, ny
        !! Number of grid points along x and y direction
        real(rp), dimension(ny,nx), intent(in) :: xgrd
        !! X-coordinate value of points on the grid
        real(rp), dimension(ny,nx), intent(in) :: ygrd 
        !! Y-coordinate value of points on the grid
        real(rp), dimension(ny,nx), intent(in) :: v
        !! V(x, y) of points on the grid
        real(rp), dimension(ny,nx), intent(in) :: vx 
        !! \(\frac{\partial V}{\partial x}\) of points on the grid
        real(rp), dimension(ny,nx), intent(in) :: vy 
        !! \(\frac{\partial V}{\partial y}\) of points on the grid
        real(rp), dimension(ny,nx), intent(in) :: vxy
        !! \(\frac{\partial^2 V}{\partial x \partial y}\) of points on the grid

        integer(ip) :: ix, iy, ii, jj
        logical :: done

        real(rp), parameter :: A(4,4) = reshape([ 1.0,  0.0, -3.0,  2.0, &
                                                  0.0,  0.0,  3.0, -2.0, &
                                                  0.0,  1.0, -2.0,  1.0, &
                                                  0.0,  0.0, -1.0,  1.0],&
                                                 shape(A))

        real(rp) :: alpha(4,4), f(4,4), xx(4), yy(4), dxxdx(4), dyydy(4), &
                    deltax, deltay, dx, dy

        done = .false.
        ix = 0
        iy = 0
        do ix=1, nx-1
            do iy=1, ny-1
                if(x < xgrd(ix+1,iy) .and. x > xgrd(ix,iy) .and. &
                   y > ygrd(ix,iy) .and. y < ygrd(ix,iy+1)) then 
                    done = .true.
                end if
                
                if(done) exit
            end do
            if(done) exit
        end do
        
        if(.not. done) then
            write(*, *) "Cannot find the required point on the grid"
            stop
        end if
        
        dx = xgrd(ix+1,iy)-xgrd(ix,iy)
        dy = ygrd(ix,iy+1)-ygrd(ix,iy)
        do ii=0, 1
            do jj=0, 1
                f(jj+1,ii+1) = v(ix+jj,iy+ii)
                f(jj+1,ii+3) = vy(ix+jj,iy+ii) * dy
                f(jj+3,ii+1) = vx(ix+jj,iy+ii) * dx
                f(jj+3,ii+3) = vxy(ix+jj,iy+ii) * dx*dy
            end do
        end do

        alpha = matmul(A, matmul(f, transpose(A)))

        xx(1) = 1.0
        yy(1) = 1.0
        deltax = (x-xgrd(ix,iy)) / dx
        deltay = (y-ygrd(ix,iy)) / dy
        do ii=2, 4
            xx(ii) = xx(ii-1) * deltax
            yy(ii) = yy(ii-1) * deltay
        end do

        dxxdx(1) = 0.0
        dyydy(1) = 0.0
        do ii=2, 4
            dxxdx(ii) = (ii-1) * xx(ii-1) / dx
            dyydy(ii) = (ii-1) * yy(ii-1) / dy
        end do

        z = 0.0
        do ii=1,4
            do jj=1, 4
                z = z + alpha(ii, jj) * yy(jj) * xx(ii)
            end do
        end do
        
        dzdx = 0.0
        dzdy = 0.0
        do ii=1,4
            do jj=1, 4
                dzdx = dzdx + alpha(ii, jj) * yy(jj) * dxxdx(ii)
                dzdy = dzdy + alpha(ii, jj) * dyydy(jj) * xx(ii)
            end do
        end do
    end subroutine compute_bicubic_interp