module wei_mig_standard_mod
  use wei_wavespace_mod
  use wei_slow_mod
  use wei_omp_mod
  use wei_buffer_mod
  use sep3d_struct_mod
  use  wei_mig_shared_mod
  use wei_sg_basic_mod
  use wei_wem_ssf_mod
  use wei_fft_mod
  use wei_dsr_mod
  use wei_work_mod
  use wei_file_mod
  use verb_mod
  use sep

  implicit none
  character(len=128), private, save :: wem_type
  
  



contains

  subroutine mig_init(wave,image,zero,freq)
    type(wspace),pointer   :: wave
    type(wfile)   :: image,freq,zero
    type(sep3d) :: file
    integer     :: n
    integer     :: esize,i,nsz
    real        :: o,d
    logical     :: image_real,add,xyy,verb
    character(len=10024) :: tt
#ifdef DEBUG
    write(0,*) "IN MIG INIT" 
#endif


   call set_no_putch()
   call from_param("verb",verb,.true.)
   call from_param("add",add,.false.)

   if(.not. add) then

      call auxin("Data")

     call grab_history("Data",tt,10024,n)

     call auxclose("R")

     call auxputlin("R",tt)

     call srite("R",4,4)
     call auxclose("R")
     call set_alternate_putch("R")
     call set_yes_putch()
   end if
  
   call auxinout("R")

      
   if(.not. init_verbose()) call seperr("Trouble initializing verbosity")

 

   
   !SETUP DATA
   call init_sep3d("Data",file,"INPUT")             ! Read in data description
   call auxclose("Data")
  
   call sep3d_to_wavespace(file,wave)                     ! Set first four axes
   call set_reg_freq(wave,file%n(5),file%o(5),file%d(5))   ! Set frequency axis


 
   call from_param("mem_save",wave%mem_save,0)  !How important is to save memory

   call from_param("wem_type",wem_type,"WEM")
  
   !Read in regular depth axis description
   call from_param("az__n",n);
   call from_param("az__o",o);
   call from_param("az__d",d)
 
   call set_reg_dz(wave,n,o,d)  !Setup regular depth sampling

   
   !Setup data file
   call init_file_tag(freq,wave,"Data",.false.,.true.,data_axes)

   !Initialize basic imaging
   call wei_sg_basic_init(wave)
   
  
   call check_setup_mig(wave,image,zero)




   wave%adj=.true. !migration
   wave%opsign=-1  !migration

   if(.not. init_sep_timers()) call seperr("trouble initializing timers") 


#ifdef DEBUG
    write(0,*) "OUT MIG INIT" 
#endif
  end subroutine


  subroutine mig_of(wave,iofile,izfile,freq,wfld,image)
    type(wspace),pointer  :: wave
    type(wfile)   :: izfile,freq,iofile
    complex        :: wfld(:,:,:,:,:)
    complex        :: image(:,:,:,:,:)
    type(westep),pointer   :: steps(:)
    type(westep),pointer   :: step
    integer              :: ith,i5,i1,i2,i3,i4,nkxko,i
    integer              :: nzs,nfreqs,nwblock,nzblock,nwb,nzb,iz
    integer,pointer      :: iws(:,:),izs(:,:)
    logical              :: wblock,zblock,endit,logit
    integer              :: t1,t2,t3,t4,t5
#ifdef OMP
    integer, external    :: omp_get_thread_num
#endif
#ifdef DEBUG
    write(0,*) "IN MIG OF" 
#endif

    call initial_work(wave,iws,izs)
    allocate(steps(size(wfld,5)))
    logit= setup_next_timer("IO READ",t1)
    logit= setup_next_timer("IO WRITE",t4)
    logit= setup_next_timer("IO SLOW",t5)
    logit= setup_next_timer("step",t2)
    logit= setup_next_timer("image",t3)
    call from_param("nkxko",nkxko,2001)


   
    if(.not. init_print_progress(0,size(wave%dz)*size(wave%ws))) &
      call seperr("trouble initializing pring progress")

   !INITIALIZE FFT
   if(my_verb()) write(0,*) "Initialize FFTs"
   call wei_fft_init_ser("FFTW",wfld)




    #ifdef OMP
    !$OMP PARALLEL default(shared) private(ith,step,i5,nwb,nzb,iz)
    #endif

    call wei_ith(ith)
    step=>steps(ith)
    step%ithread=ith
    call init_step(wave,step,ith,nkxko)
    call wem_ssf_init(wfld,step,wem_type)

    if(ith==1) then
      if(izfile%doit) izfile%wave=wave
      iofile%doit=.true.
      freq%wave=wave
      if(.not. iofile%add ) then
        if(my_verb()) write(0,*) "Zeroing image"
        call wei_zero_i5(image,1)
        call sep_open_lock(iofile%tag)
        do i5=1,size(wave%dz)
          call write_i5(iofile,i5,1,image)
          if(izfile%doit) then
             write(0,*) "DO IZ WRITE",i5
             call write_i5(izfile,i5,1,image)
          end if
        end do
        call sep_close_lock(iofile%tag)
      end if
    end if

    call sep_begin_prog()



    do while(grab_next_ws_mig(step,iws,nwb)) !LOOP OVER FREQUENCY BLOCKS
      if(step%ithread==1) call start_timer_num(t1)
      if(step%ithread==1) then         !MASTER THREAD READS DATA

    do i5=1,nwb          !LOOP OVER FREQUENCIES
          call read_i5(freq,iws(i5,step%ithread),i5,wfld) !READ IN FREQUNCY
        end do
      end if
      if(step%ithread==1) call stop_timer_num(t1)

      !ASSIGN THREAD ITS FREQUENCY
      if(ith <=nwb) step%w=step%wave%ws(iws(ith,step%ithread)) 

      endit=.false.


      !LOOP OVER DEPTH BLOCKS
      do while(.not. endit .and. grab_next_dzs_mig(step,izs,nzb)) 

        !-------------------
        !ZERO IMAGE
        !-------------------
        if(step%ithread==1) then
          if(step%ithread==1) call start_timer_num(t5)
          call read_slow(step%wave,izs(1,step%ithread))  !READ IN SLOWNESS
          if(step%ithread==1) call stop_timer_num(t5)
          do i5=1,nzb
            call wei_zero_i5(image,i5)
          end do
        end if
        call wei_barrier()


        step%izb=izs(1,step%ithread)-1
        if(ith<=nwb) then !DO MIGRATION IF THIS THREAD HAS WORK
          do iz=1,nzb
            step%iz=iz 
            step%dstep=step%wave%dz(izs(iz,step%ithread))



            !-------------------
            !PERFORM THE IMAGING;
            !-------------------;
            #ifdef OMP
             !$OMP CRITICAL   !SYNC ALL THREADS
            #endif
            if(step%ithread==1) call start_timer_num(t3)
            call s_g_basic_image(step,wfld,image)
            if(step%ithread==1) call stop_timer_num(t3)
            #ifdef OMP
            !$OMP END CRITICAL   !SYNC ALL THREADS
            #endif

            !-------------------
            !DOWNWARD CONTINUE
            !-------------------
            if(step%ithread==1) call start_timer_num(t2)
            call one_step(wfld,step,wave)
            if(step%ithread==1) call stop_timer_num(t2)

          end do  !LOOP OVER DEPTH STEPS
        end if    !IF WE HAVE A FREQUENCY

       
        
        !-------------------
        !WRITE OUT IMAGE
        !-------------------
        call wei_barrier()
        if(step%ithread==1)  then
          call sep_open_lock(iofile%tag)
          do i5=1,nzb
            if(step%ithread==1) call start_timer_num(t4)
            call write_i5(iofile,izs(i5,step%ithread),i5,image)
            if(izfile%doit) call write_i5(izfile,izs(i5,step%ithread),i5,image)
            if(step%ithread==1) call stop_timer_num(t4)
            call update_finished(step%wave,size(iws,1),iofile%tag)
          end do
          call sep_close_lock(iofile%tag)
        end if
        call wei_barrier()


         
        if(step%ithread==1) call print_progress("Migration ",nwb*nzb)
        if(izs(nzb,step%ithread)==size(step%wave%dz)) then
          endit=.true.
          iofile%add=.true.
          izfile%add=.true.
        end if
      end do
      call wei_barrier() 

    end do
#ifdef OMP
!$OMP END PARALLEL
#endif


  call print_timers()
#ifdef DEBUG
    write(0,*) "OUT MIG OF" 
#endif
  end subroutine


end module
