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

module init
use def

  complex, dimension(bs) :: u=0.0, v=0.0, rho=0.0, p=0.0, omega=0.0
  complex, dimension(bs) :: ascr1=0.0, ascr2=0.0, ascr3=0.0, ascr4=0.0
  complex, dimension(bs) :: sscr1=0.0, sscr2=0.0
  complex, dimension(bs) :: ik=0.0
  real, dimension(bs) :: x=0.0, y=0.0

 

contains
  subroutine initialization(u,v,rho,ik,x,y)
  complex, dimension(bs), intent(inout) :: u, v, rho
  complex, dimension(bs), intent(inout) :: ik
  real, dimension(bs), intent(inout) :: x, y

  integer :: i,j
  integer, dimension(12) :: seed
  real :: r

  real, dimension(n) :: kk
  real, dimension(n) :: xx
  real, dimension(n) :: yy
  real, dimension(n) :: ones=1.0
  real, dimension(bs) :: k

  u=0.0; v=0.0; rho=1.0

  seed(:) = int(86345*rank+23456)
  call random_seed(put=seed)

  kk = wave_num(n)

  do i=1,bs/n
    k(1+n*(i-1):n*i)=kk
  enddo

  ik = ii*k

  xx = xgrid(n)
  yy = ygrid(n)

  do i=1,bs/n
    x(1+n*(i-1):n*i)=xx
    y(1+n*(i-1):n*i)=yy(bs*rank/n+i)*ones
  enddo

  if (rank==0) then
    do i=1,n
      call random_number(r)
      rho(i) = Ta-0.5*r
      !rho(i) = -0.5*exp(-10*(x(i)-pi)**2)+1.0
    enddo
  endif

  if (rank==nproc-1) then
    do i=bs-n+1,bs
      call random_number(r)
      rho(i) = Tb+0.5*r
    enddo
  endif


  !if (rank==0) then
  !  do i=1,bs
  !    u(i) = exp(-10*((x(i)-pi)**2+(y(i)-pi/2)**2))
  !    v(i) = exp(-10*((x(i)-pi)**2+(y(i)-pi/2)**2))
  !  enddo
  !endif
  end subroutine initialization

  
  function wave_num(s)
  integer :: s,i
  real, dimension(s) :: wave_num
  do i=1,s
     if (i<=s/2) then
        wave_num(i)=i-1
     else
        wave_num(i)=-(s-i+1)
     end if
  end do
  end function wave_num


  function xgrid(s)
  integer :: s,i
  real, dimension(s) :: xgrid
  do i=1,s
     xgrid(i)=i*dx
  end do
  end function xgrid


  function ygrid(s)
  integer :: s,i
  real, dimension(s) :: ygrid
  do i=1,s
     ygrid(i)=i*dy
  end do
  end function ygrid

end module init


