!*******************************************************************
!  
!   2D turbulent flow/ binary fluids/ surfactant
!
!   2D, pseudo-spectral for flow (with phase-field modeling)
!   2D, WENO 5th/TVD RK3/particles for level set
!
!--------------------------------------------------------------------
!   This part of the code is mostly initialization
!   
!  1) Initializes MPI stuff
!  2) Processor 0 reads the content of the file control.in 
!     which contains all the control parameters
!  3) Processor 0 broadcasts the content of the file to all the other
!     processors
!  4) Each processor allocates dynamically the space for all the arrays
!
!
!-----------------------------------------------------------------------
      program Surf
!-----------------------------------------------------------------------
       use communicators
       use physics_stuff
       use grid_stuff			

       include "mpif.h"
       integer n,n1,n2,np1,lsize
       parameter(n=100)
       parameter(n1=0256,n2=0512,np1=04)
       parameter(lsize=n1*n2/np1)
       parameter(lpadding = 128)
       real*8 params(n)
       real clok
       integer ierr,lsize,color,key

       character(LEN=80) ::  f_name='NN_BJ'
! 
!  Define the array variables and work arrays
!   
       real*8, dimension(lsize+lpadding):: &
                                           pu1,pu2,cu1,cu2,&
                                           ptt,ctt,ftt,ttt,&
					   paa,caa,faa,taa,&
					   pbb,cbb,fbb,tbb,&
					   pcc,ccc,fcc,tcc,&
					   ucc,vcc,wcc,xcc,ycc,zcc,&
					   ww1,ww2,ww3,ww4,ww5,ww6,ww7
       real*8 pfactor
	
!  Initialise MPI stuff and get processor id number
!
       gcomm = MPI_COMM_WORLD
       call MPI_INIT(ierr)
       call MPI_COMM_RANK(gcomm,nloc,ierr)
       nproc = np1

       if(nloc .eq. 0) then  
	params = 0.0
 
         open(10,file="control.in")
          do i = 1,n
            read(10,*)params(i)
          end do
         close(10)
       end if  
!
!  Content of file control.in is broadcasted to all the processors
!
       params(1) = n1   ! size along the 1st dimension in cs (x)
       params(2) = n2   ! size along the 1st dimension in cs (y)
       params(4) = np1  ! Number of procs along first proc dim

      call bcast(params,n,gcomm)
!
!  Begin the initialization of the MODULES
!

       call setup_confg(n1,n2)  ! module domain_confg
       call setup_phase(n1,n2)  ! module domain_phase
       call setup_fft()                    ! module fft_stuff
       call setup_wave()                   ! module wave_stuff
       call IO_setup_phase()               ! module used by Romio
       call IO_setup_confg()               ! module used by Romio
       call setup_physics(params,n)        ! module physics_stuff
       call setup_grid(params,n)           ! module grid_stuff 
       call setup_steps(params,n)          ! module steps_stuff
       call setup_deriv()                  ! module deriv_stuff
       call setup_diffuse()                ! module diffuse_stuff
       call setup_statistics()		   ! module statisitics_stuff
       call setup_alias(params,n)          ! module alias_stuff
       call setup_output(params,n,f_name)      ! module output_stuff
       call setup_input(params,n,f_name,nloc)  ! module input_stuff

	alpha = alpha * (dx)**1.
	gc    = gc * dx
	pfactor = log(2./0.01-1.)/2./gc
	beta  = 2.* alpha * pfactor**2.
!	beta = 2.*20.*dx*pfactor**2.
!
!  Ready to start
!
!---------measure time, get the initial time------------------
      clok = MPI_WTIME()
!YN
	if(nloc.eq.0) print *,'time = ',clok
!YN
!--------Enable streams for speed--------------------------------
!      call set_d_stream(1)        

      call drive(pu1,pu2,ptt,paa,pbb,pcc, &
	         cu1,cu2,ctt,caa,cbb,ccc, &
                         ftt,faa,fbb,fcc, &
                 	 ttt,taa,tbb,tcc, &
	         ucc,vcc,     &
		 wcc,xcc,ycc,zcc,&	
                 ww1,ww2,ww3,ww4,ww5,ww6,ww7)

!--------Disable streams ------------------------------------
!      call set_d_stream(0)    

      clok = MPI_WTIME() - clok

      if(nloc.eq.0) then
	write(*,*)'execution time ',clok
      end if
!
!  Close shop and get out
!     

       call MPI_FINALIZE(ierr)
 
      end program Surf