module xbd3d_module
   
   implicit none
   
contains

   subroutine xbeach_in_delft3d(xd3d,yd3d,dps,s1,u1,v1,                  &
                                uorb_d3d,tp_d3d,teta_d3d,dis_d3d,        &
                                wsu_d3d,wsv_d3d,fxw_d3d,fyw_d3d,         &
                                grmasu_d3d,grmasv_d3d,                   &
                                hrms_d3d,wlen_d3d,                       &
                                mmax,nmax,nmaxus,timnow)
                                
      use libxbeach_module
      use initialize_module
      
      implicit none
   
      integer,                     intent(in)             :: mmax,nmax,nmaxus
      real*8,                      intent(in)             :: timnow
      real*8,dimension(nmax,mmax), intent(in)             :: xd3d,yd3d,dps,s1,u1,v1                  
      real*8,dimension(nmax,mmax), intent(out)            :: uorb_d3d,tp_d3d,teta_d3d,dis_d3d,        &
                                                             wsu_d3d,wsv_d3d,fxw_d3d,fyw_d3d,         &
                                                             grmasu_d3d,grmasv_d3d,                   &
                                                             hrms_d3d,wlen_d3d
   
      integer                                             :: rc
      logical, save                                       :: firsttimexbeach = .true.
      real*8, save                                        :: tend
      real*8                                              :: tcurrent
      
#ifdef USEMPI
      integer, dimension(12)                              :: info
      character(256)                                      :: line
      integer                                             :: rank,i
#endif
      
      if (firsttimexbeach) then
         rc = 0
         error   = 0
         n = 0
   
         ! setup of MPI
#ifdef USEMPI
         s=>slocal
         call xmpi_initialize
         call xmpi_barrier(toall)
         t0 = MPI_Wtime()
#endif
   
         ! create log files
         call start_logfiles(error)
   
         ! set starting time and date
         call cpu_time(tbegin)
         tend = tbegin
         
         ! show startup message
         call writelog_startup()
   
         !-----------------------------------------------------------------------------!
         ! Initialize simulation                                                       !
         !-----------------------------------------------------------------------------!
   
         ! initialize time counter
         it      = 0
   
         ! read input from params.txt
         params_inio = .false.
         call all_input(par)
         
         ! overwrite flags for safety
         par%gwflow = 0
         par%flow = 0
         par%nonh = 0
         par%ships = 0
         par%ndrifter = 0
         par%sedtrans = 0
         par%bchwiz = 0
         par%morphology = 0
         par%setbathy = 0
         par%xbd3d = 1
   
         ! allocate space scalars
         call space_alloc_scalars(sglobal)
         s => sglobal
   
         ! read grid and bathymetry
         call grid_bathy_in_d3d(xd3d,yd3d,dps,mmax,nmax,nmaxus,s,par)
   
         ! distribute grid over processors
#ifdef USEMPI
         call xmpi_determine_processor_grid(s%nx,s%ny,par%mpiboundary,par%mmpi,par%nmpi,par%cyclic,error)
         call writelog_mpi(par%mpiboundary,error)
#endif
   
         ! initialize timestep
         call timestep_init(par, tpar)
   
         if (xmaster) then
   
            call writelog('ls','','Initializing .....')
         endif
         ! initialize physics
         call readtide           (s,par)
         call readwind           (s,par)
         call flow_init          (s,par)
         call discharge_init     (s,par)
         call drifter_init       (s,par)
         call wave_init          (s,par)
         call gw_init            (s,par)
         ! TODO, fix ordening of arguments....
         call bwinit             (s)          ! works only on master process
   
         call sed_init           (s,par)
   
         call ship_init          (s,par,sh)   ! always need to call initialise in order
         ! to reserve memory on MPI subprocesses.
         ! Note: if par%ships==0 then don't allocate
         ! and read stuff for sh structures
         call veggie_init         (s,par,veg)
   
#ifdef USEMPI
         call distribute_par(par)
         s => slocal
         !
         ! here an hack to ensure that sglobal is populated, also on
         ! the not-(o)master processes, just to get valid addresses.
         !
         if (.not. xmaster .and. .not. xomaster) then
            call space_alloc_arrays_dummies(sglobal,par)
         endif
         call space_distribute_space (sglobal,s,par     )
#endif
   
         call ranges_init(s)
   
         ! nonh_init does not always need to be called
         if (par%nonh==1) call nonh_init(s,par)
   
         ! initialize output
         call means_init             (sglobal,s,par     )
   
         call output_init            (sglobal,s,par,tpar)
   
   
         ! store first timestep
         ! from this point on, xomaster will hang in subroutine output
         ! until a broadcast .true. is received
         call output(sglobal,s,par,tpar)
         firsttimexbeach = .false.
      end if
   
      ! Start simulation                                                            !
      if (par%instat==0.or.par%instat==40) then
         par%tnext = timnow
      end if
      par%tstop = timnow
      tcurrent = tend
      !
      ! Interpolate bathymetry
      !
      call delft3d_inputfields(s,par,dps,s1,u1,v1,mmax,nmax,nmaxus) 
      !
      do while (par%t<par%tstop)
#ifdef USEMPI
         if (execute_counter .eq. 1) then
            ! exclude first pass from time measurement
            call xmpi_barrier
            t01 = MPI_Wtime()
         endif
#endif
         execute_counter = execute_counter + 1
   
         call outputtimes_update(par, tpar)
   
         ! determine timestep
         if(xcompute) then
   
            ! determine this time step's wet points
            call compute_wetcells(s,par)
            !
            ! determine time step
            call timestep(s,par,tpar,it,ierr=error)
                 
            if (error==0) then
               
               ! boundary conditions
               call wave_bc        (sglobal,s,par)
   
               ! compute timestep
               if (par%swave==1)        call wave           (s,par)
               if (par%vegetation==1)   call vegatt         (s,par,veg)
            endif
         endif
   
         n = n + 1
   
         ! output
         rc = outputext()
      enddo
      
      call delft3d_outputfields(s,par,uorb_d3d,tp_d3d,teta_d3d,dis_d3d,  &
                                wsu_d3d,wsv_d3d,fxw_d3d,fyw_d3d,           &
                                grmasu_d3d,grmasv_d3d,                   &
                                hrms_d3d,wlen_d3d,                       &
                                mmax,nmax,nmaxus)
                                
#ifdef USEMPI
      end_program = .true.
      call xmpi_send_sleep(xmpi_imaster,xmpi_omaster) ! wake up omaster
      call xmpi_bcast(end_program,toall)
      call xmpi_barrier(toall)
      call writelog_finalize(tbegin,n,par%t,s%nx,s%ny,t0,t01)
      call xmpi_finalize
#else
      if (xmaster) then
         call cpu_time(tend)
         write(*,*)'XBeach calculation time for this timestep: ',tend-tcurrent,' seconds'
         write(*,*)'Total XBeach calculation time: ',tend-tbegin,' seconds'
         write(*,*)'End of XBeach part... Switching back to Delft3D.'
      end if
#endif
   end subroutine xbeach_in_delft3d
      
end module