module wei_datum_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_damp_mod
  use wei_dsr_mod
  use wei_work_mod
  use wei_file_mod
  use verb_mod
  use sep

  implicit none


contains

  subroutine datum_init(wave,datum,freq)
    type(wspace),pointer   :: wave
    type(wfile)   :: datum,freq
    type(sep3d) :: file
    integer     :: n
    integer     :: esize,i
    real        :: o,d
    logical     :: add,xyy

#ifdef DEBUG
    write(0,*) "IN DATUM INIT" 
#endif


   
   call set_alternate_putch("Ddown")
   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 init_damping(wave)



   !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

   
   call setup_sg_axes()   !We are doing shot receiver migration

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



   call from_param("add",add,.false.)
   if(add) then
     call init_file_tag(datum,wave,"Ddown",.false.,.true.,data_axes)
     datum%add=.true.
   else
     call init_file_struct(datum,wave,"Ddown",.false.,.true.,data_axes)
     datum%add=.false.
   end if

   !Initialize basic imaging
   call wei_sg_basic_init(wave)

    
   call from_param("nzb",wave%nzb,1) !Depth block
!   call init_one_slow(wave,"S","ref_slow")

   call dsr_init()  !Initialize lookup table
   call weikmap_init(wave)  !Initialize kmaps 

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



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

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


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

    call initial_work(wave,iws,izs)
    allocate(steps(size(wfld,5)))
    call from_param("nkxko",nkxko,2001)

    if(.not. dfile%add ) then
      if(my_verb()) write(0,*) "Zeroing datum"
      call wei_zero_i5(wfld,1)
      do i5=1,size(wave%ws)
        call write_i5(dfile,i5,1,wfld)
      end do
    end if

    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)
    call init_step(wave,step,ith,nkxko)
    call wem_ssf_init(wfld,step)

    call sep_begin_prog()


    do while(grab_next_ws_mig(step,iws,nwb)) !LOOP OVER FREQUENCY BLOCKS
      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

      !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 DATUM
        !-------------------
        if(step%ithread==1) then
          call read_slow(step%wave,izs(1,step%ithread))  !READ IN SLOWNESS
        end if
        call wei_barrier()

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

            !-------------------
            !DOWNWARD CONTINUE
            !-------------------
            call one_step(wfld,step,wave)

          end do  !LOOP OVER DEPTH STEPS

          if(izs(nzb,step%ithread)==size(step%wave%dz)) then
            endit=.true.
          end if
!      call wei_barrier()
!!if(step%ithread==1)  write(0,*) "CHECK THIS",nwb,nzb,iws(1,1),izs(1,1)
      call wei_barrier()
          if(step%ithread==1) call update_finished(step%wave,nwb*nzb,dfile%tag)
          if(step%ithread==1) call print_progress("Datumed ",nwb*nzb)

        end if    !IF WE HAVE A FREQUENCY


      end do !LOOP OVER Z

!      call wei_barrier()
!     write(0,*) "WHAT THE",nwb,nzb,iws(:,1),step%ithread
!      call wei_barrier()
       
      !----------------------
      !WRITE OUT DATUMED DATA
      !----------------------
!     write(0,*) "IN 1",step%ithread
!      call wei_barrier()
!     write(0,*) "OUT 1",step%ithread
      call wei_barrier()
      if(step%ithread==1)  then
        do i5=1,nwb
          call write_i5(dfile,iws(i5,step%ithread),i5,wfld)
        end do
      end if
      call wei_barrier()
!     write(0,*) "YO 1",step%ithread
!      call wei_barrier()
!  write(0,*) "WHAT2THE",nwb,nzb,iws(:,1),step%ithread,wave%finished
!      call wei_barrier()

    end do !LOOP OVER W

!    call wei_barrier() 

#ifdef OMP
!$OMP END PARALLEL
#endif


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


end module wei_datum_standard_mod
