!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! RayleighBernard2D 
! advec.f90 module
! James Sandham 
! 11 Mar 2014
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


module advec
use def
use fft
use mpi

contains
  subroutine advection(uh,vh,uc,up,vc,vp,ik)
  complex, dimension(bs), intent(inout) :: uh,vh
  complex, dimension(bs), intent(inout) :: uc,up,vc,vp
  complex, dimension(bs), intent(in) :: ik

  integer :: status(MPI_STATUS_SIZE)
  integer :: i,j,k,tag=100
  complex, dimension(bs) :: ur=0.0,vr=0.0,uur=0.0,uvr=0.0,vvr=0.0
  complex, dimension(bs) :: uuh=0.0,uvh=0.0,vvh=0.0
  complex, dimension(n) :: lower_bc1=0.0,upper_bc1=0.0,lower_bc2=0.0,upper_bc2=0.0

  call ifft1(uh,ur)
  call ifft1(vh,vr)

  do i=1,bs
    uur(i) = ur(i)*ur(i)
    uvr(i) = ur(i)*vr(i)
    vvr(i) = vr(i)*vr(i)
  enddo

  call fft1(uur,uuh)
  call fft1(uvr,uvh)
  call fft1(vvr,vvh)

  if (mod(rank,2)==0) then
    if (rank==0) then
      call MPI_SEND(uvh((bs-n)+1:bs),n,MPI_COMPLEX,rank+1,tag,MPI_COMM_WORLD,ierr)
      call MPI_RECV(upper_bc1,n,MPI_COMPLEX,rank+1,tag,MPI_COMM_WORLD,status,ierr)

      call MPI_SEND(vvh((bs-n)+1:bs),n,MPI_COMPLEX,rank+1,tag,MPI_COMM_WORLD,ierr)
      call MPI_RECV(upper_bc2,n,MPI_COMPLEX,rank+1,tag,MPI_COMM_WORLD,status,ierr)
    else
      call MPI_SEND(uvh((bs-n)+1:bs),n,MPI_COMPLEX,rank+1,tag,MPI_COMM_WORLD,ierr)
      call MPI_SEND(uvh(1:n),n,MPI_COMPLEX,rank-1,tag,MPI_COMM_WORLD,ierr)
      call MPI_RECV(lower_bc1,n,MPI_COMPLEX,rank-1,tag,MPI_COMM_WORLD,status,ierr)
      call MPI_RECV(upper_bc1,n,MPI_COMPLEX,rank+1,tag,MPI_COMM_WORLD,status,ierr)

      call MPI_SEND(vvh((bs-n)+1:bs),n,MPI_COMPLEX,rank+1,tag,MPI_COMM_WORLD,ierr)
      call MPI_SEND(vvh(1:n),n,MPI_COMPLEX,rank-1,tag,MPI_COMM_WORLD,ierr)
      call MPI_RECV(lower_bc2,n,MPI_COMPLEX,rank-1,tag,MPI_COMM_WORLD,status,ierr)
      call MPI_RECV(upper_bc2,n,MPI_COMPLEX,rank+1,tag,MPI_COMM_WORLD,status,ierr)

    endif
  else
    if (rank==nproc-1) then
      call MPI_RECV(lower_bc1,n,MPI_COMPLEX,rank-1,tag,MPI_COMM_WORLD,status,ierr)
      call MPI_SEND(uvh(1:n),n,MPI_COMPLEX,rank-1,tag,MPI_COMM_WORLD,ierr)

      call MPI_RECV(lower_bc2,n,MPI_COMPLEX,rank-1,tag,MPI_COMM_WORLD,status,ierr)
      call MPI_SEND(vvh(1:n),n,MPI_COMPLEX,rank-1,tag,MPI_COMM_WORLD,ierr)
    else
      call MPI_RECV(lower_bc1,n,MPI_COMPLEX,rank-1,tag,MPI_COMM_WORLD,status,ierr)
      call MPI_RECV(upper_bc1,n,MPI_COMPLEX,rank+1,tag,MPI_COMM_WORLD,status,ierr)
      call MPI_SEND(uvh((bs-n)+1:bs),n,MPI_COMPLEX,rank+1,tag,MPI_COMM_WORLD,ierr)
      call MPI_SEND(uvh(1:n),n,MPI_COMPLEX,rank-1,tag,MPI_COMM_WORLD,ierr)

      call MPI_RECV(lower_bc2,n,MPI_COMPLEX,rank-1,tag,MPI_COMM_WORLD,status,ierr)
      call MPI_RECV(upper_bc2,n,MPI_COMPLEX,rank+1,tag,MPI_COMM_WORLD,status,ierr)
      call MPI_SEND(vvh((bs-n)+1:bs),n,MPI_COMPLEX,rank+1,tag,MPI_COMM_WORLD,ierr)
      call MPI_SEND(vvh(1:n),n,MPI_COMPLEX,rank-1,tag,MPI_COMM_WORLD,ierr)
    endif
  endif

  do i=1,bs
    j=i-n ; k=i+n
    if (j<=0) then
      if (rank==0) then
        uc(i) = 0.0
        vc(i) = 0.0 
      else
        uc(i) = -ik(i)*uuh(i)-(uvh(k)-lower_bc1(i))/(2*dy)
        vc(i) = -ik(i)*uvh(i)-(vvh(k)-lower_bc2(i))/(2*dy)
      endif
    elseif (k>bs) then
      if (rank==nproc-1) then
        uc(i) = 0.0
        vc(i) = 0.0
      else
        uc(i) = -ik(i)*uuh(i)-(upper_bc1(k-bs)-uvh(j))/(2*dy)
        vc(i) = -ik(i)*uvh(i)-(upper_bc2(k-bs)-vvh(j))/(2*dy)
      endif
    else
      uc(i) = -ik(i)*uuh(i)-(uvh(k)-uvh(j))/(2*dy)
      vc(i) = -ik(i)*uvh(i)-(vvh(k)-vvh(j))/(2*dy)
    endif
  enddo


  if (tt==1) then
    uh = uh + dt*uc
    vh = vh + dt*vc
  else
    uh = uh + (dtl*uc - dts*up)
    vh = vh + (dtl*vc - dts*vp)
  endif
  up = uc; vp = vc;

  end subroutine advection


end module advec

