I'm trying to write some Fortran code to achieve a bilinear
interpolation upon a grid of dimension 240 values (x) by 120 values
(y) and where the data is currently at a spatial resolution of 1.5
degrees... I want the output grid to be at a resolution of 5 minutes.

How would I go about doing this, and what would the Fortran code look
like to achieve this?

subroutine grid_bilin(old,ns,new)
real,intent(in):: old(0:,0:) ! old data
integer,intent(in):: ns ! resolution factor (common for x and y)
real,intent(out):: new(0:,0:)
real:: im(0:ns,2) ! interpolation matrix (common for x and y)
integer:: i,j

! prepare interpolation matrix
forall(i=0:ns) im(i,2) = real(i,kind(im))/ns
im(:,1) = 1 - im(:,2)

do j=1,ubound(old,2)
do i=1,ubound(old,1)
new((i-1)*ns:i*ns,(j-1)*ns:j*ns) =
matmul(im,matmul(old(i-1:i,j-1:j),transpose(im)))
end do
end do
end subroutine