!
! Derivative 1, MPI
!
subroutine FirstDeriv1dp(npts, dx, u, u_x, mynode, totalnodes)
  implicit none
  integer, intent(in) :: npts, mynode, totalnodes
  real*8, intent(in) :: dx
  real*8, dimension(:), intent(in) :: u
  real*8, dimension(:), intent(out) :: u_x
  real*8 two_invdx
  MPI_STATUS status
  two_invdx = 1d0/(2d0*dx)

  if(mynode == 0) then
     ! forward...
     u_x(1) = (-3.*u(1)+4.*u(2)-u(3))*two_invdx
  end if

  if(mynode == totalnodes) then
     ! backward
     u_x(npts) = (3.*u(npts)-4.*u(npts-1)+u(npts-2))*two_invdx
  end if

  do i=2,npts-1
     ! central...
     u_x(i) = (u(i+1)-u(i-1))*two_invdx
  end do

  if (mynode == 0) then
     mpitemp = u(npts)
     call MPI_SEND(mpitemp, 1, MPI_DOUBLE, 1, 1, MPI_COM_WORLD)
     call MPI_RECV(mpitemp, 1, MPI_DOUBLE, 1, 1, MPI_COM_WORLD, status)
     u_x(npts) = (mpitemp-u(npts-1))*two_invdx
  else if (mynode == totalnodes) then
     call MPI_RECV(mpitemp, 1, MPI_DOUBLE, mynode, 1, MPI_COM_WORLD, status)
     u_x(1) = (u(2)-mpitemp)*two_invdx
     mpitemp = u(1)
     call MPI_SEND(mpitemp, 1, MPI_DOUBLE, mynode, 1, MPI_COM_WORLD)
  else
     call MPI_RECV(mpitemp, 1, MPI_DOUBLE, mynode, 1, MPI_COM_WORLD, status)
     u_x(1) = (u(2)-mpitemp)*two_invdx
     mpitemp = u(1)
     call MPI_SEND(mpitemp, 1, MPI_DOUBLE, mynode-1, 1, MPI_COM_WORLD)
     mpitemp = u(npts)
     call MPI_SEND(mpitemp, 1, MPI_DOUBLE, mynode+1, 1, MPI_COM_WORLD)
     call MPI_RECV(mpitemp, 1, MPI_DOUBLE, mynode+1, 1, MPI_COM_WORLD, status)
     u_x(npts) = (mpitemp-u(npts-1))*two_invdx
  end if
end subroutine FirstDeriv1dp
