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



module diff
use def
use fft
use solvers

contains
subroutine diffusion(uh,vh)
complex, dimension(bs), intent(inout) :: uh,vh

integer :: i,j,k,ii
real, dimension(n) :: a,b,c
complex, dimension(n) :: uht,vht
  

call mpitranspose(uh)
call mpitranspose(vh)

a=-dt*kappa/dy2 ; c=-dt*kappa/dy2
a(1)=0.0; a(n)=0.0; c(1)=0.0; c(n)=0.0

do i=1,bs/n
  uht = uh(n*(i-1)+1:n*i)
  vht = vh(n*(i-1)+1:n*i)

  j = rank*bs/n+i
  if (j<=n/2) then
     k=j-1
  else
     k=-(n-j+1)
  end if

  do ii=2,n-1
    uht(ii) = uht(ii)-0.5*kappa*dt*(k**2*uht(ii)-(uht(ii+1)-2*uht(ii)+uht(ii-1))/dy2)
    vht(ii) = vht(ii)-0.5*kappa*dt*(k**2*vht(ii)-(vht(ii+1)-2*vht(ii)+vht(ii-1))/dy2)
  enddo

  b=1.0+dt*kappa*(2.0/dy2+k**2); b(1)=1.0; b(n)=1.0
  call tridiag(uht,uh(n*(i-1)+1:n*i),a,b,c)
  call tridiag(vht,vh(n*(i-1)+1:n*i),a,b,c)
  
  !b=1.0+dt*kappa*(2.0/dy2+k**2) 
  !b(1)=1.0; b(n)=1.0
  !call tridiag(uh(n*(i-1)+1:n*i),uh(n*(i-1)+1:n*i),a,b,c)
  !call tridiag(vh(n*(i-1)+1:n*i),vh(n*(i-1)+1:n*i),a,b,c)
enddo

call mpitranspose(uh)
call mpitranspose(vh)
end subroutine diffusion


end module diff
