module wei_mig_base
  use wei_basics_mod
  use wei_slow_mod
  use wei_imaging_mod
  use wei_data_mod
 use wei_cig
  use wei_image_mod
  use wei_wavespace_mod
  use wei_fft_mod
  use wei_wc_mod
  use wei_kmap_mod
  use wei_xmap_mod
  use wei_z_mod
  implicit none
  type(sep3d), private, save :: wsep,ssep
  complex, pointer :: wfld(:,:,:,:,:),image(:,:,:,:,:)


contains
!!$
!!$=head1 
!!$
!!$ weimig_init - Initialize migration operator
!!$
!!$=head1 SYNOPSIS
!!$
!!$  call weimig_init(SLin,WCin,FKin,FXin,IGin,in6,in7,in8,in9)
!!$
!!$=head1 PARAMETERS
!!$
!!$=over 4
!!$
!!$=item SLin - Slowness operator init  L<weislo1_init>
!!$
!!$=item WCin - WCop operator init  L<weikxmig1_init>
!!$
!!$=item FKin - FK operator init  L<weiwem_init>
!!$
!!$=item FXin - FX operator init  L<weissf_init>
!!$
!!$=item IGin - image gather creation operator init  L<weihcig_init>
!!$
!!$=item in6,..in9 -  integer,optional
!!$
!!$ Dummy operators for future expansion
!!$
!!$=back
!!$
!!$=head1 DESCRIPTION
!!$
!!$ Initialize migration/modeling operator
!!$
!!$=head1 SEE ALSO
!!$
!!$ L<weimig> L<wei>
!!$
!!$=head1 LIBRARY
!!$
!!$ sepweif90
!!$
!!$=cut
!!$
  !----------------------------------------------------------------
  subroutine weimigz_init(SLin,WCin,FKin,FXin,IGin,in7,in8,in9)
    integer,external,optional ::                     in7,in8,in9
    integer :: st,er,ierr
    type(sep3d) :: dsep
    type(sep3d) :: rsep
    character(len=128) :: stag
    logical  :: cR,cD
#ifdef DEBUG
    character(len=128), parameter :: unit='weimig_init'
#endif
    interface !!-------------------------
       integer function    SLin() result(st)
       end function SLin
       integer function WCin() result(st)
       end function WCin
       integer function    FKin() result(st)
       end function FKin
       integer function    FXin() result(st)
       end function FXin
       integer function    IGin() result(st)
       end function IGin
    end interface !!---------------------

    if(.not. init_wei_basics()) &
      call seperr("trouble initializing wei basics")

    !INITIALZE WAVEFIELD SPACE
    if(my_verb()) write(0,*) "BEFORE INITIALIZE WAVE SPACE"
    if(.not. init_reg_wave_space(wsep)) &
      call seperr("trouble initializing wave space")
    call store_current_wave_space(wsep)

    !INITIALIZE IMAGE
    if(my_verb()) write(0,*) "BEFORE INITIALIZE IMAGE"
    if(.not. init_image("R",rsep,CR)) &
      call seperr("trouble initializing image")

    !INITIALIZE DATA
    if(my_verb()) write(0,*) "BEFORE INITIALIZE DATA"
    if(.not. init_data("D",dsep,CD)) &
      call seperr("")

    if(.not. operation_check(cD,cR)) &
      call seperr("")
      

    !INITIALZE SLOWNESS
    if(my_verb()) write(0,*) "BEFORE INITIALIZE SLOWNESS"
    if(.not. init_slow()) &
      call seperr("")

    !INITIALZE SLOWNESS FROM FILE
    if(my_verb()) write(0,*) "BEFORE INITIALIZE SLOWNESS FILE"
    if(.not. init_slow_tag("S",stag,ssep)) &
      call seperr("")

    call store_slowness_masters((/ssep/),(/stag/))

    !INITIALIZE SLOWNESS OPERATOR
    if(my_verb()) write(0,*) "BEFORE INITIALIZE SLOWNESS OPERATOR"
    st=SLin()     !! slownesses init

    !INITIALIZE CIG OPERATOR
    if(my_verb()) write(0,*) "BEFORE INITIALIZE IMAGING OPERATOR"
    st=IGin()     !! cig op init

    !INITIALIZE KMAPS
    if(my_verb()) write(0,*) "BEFORE INITIALIZE MAPS"
    call weikmap_init(wsep)
    call weixmap_init(wsep)
    call show_mem("map")

    if(my_verb()) write(0,*) "BEFORE INITIALIZE FK"
    st=FKin()     !! FK op init
    if(my_verb()) write(0,*) "BEFORE INITIALIZE FX"
    st=FXin()     !! FX op init
    if(my_verb()) write(0,*) "BEFORE INITIALIZE WC"
    st=WCin()     !! WCop init
    if(my_verb()) write(0,*) "BEFORE INITIALIZE PATH"
    st=Path_z_init()   !! Path init

    !INTIALIZE IMAGING OPERATOR
    if(my_verb()) write(0,*) "BEFORE IMAGING SUB INIT"
    call weiimage_init() 

    if(my_verb()) write(0,*) "BEFORE SPACE LOGIC "
    if(.not. check_space_logic(wsep,dsep,rsep,(/ssep/))) &
    	call seperr("spaces don't conform" )
     
    if(my_verb()) write(0,*) "BEFORE MEMORY ALLOCATION"
    if(.not. wem_allocate_memory(wsep,rsep,my_nomp()))&
      call seperr("trouble allocating memory")


    !INITIALIZE FFTS after memory allocation
    if(my_verb()) write(0,*) "BEFORE FFT INIT"
    call assign_wave_space(wfld)  !get memory for wfld
    call weifft_init(wfld)  
    if(my_verb()) write(0,*) "BEFORE STRUCTURE DELETE"
    call sep3df_delete(dsep,"D")
    call sep3df_delete(rsep,"R")

#ifdef DEBUG
    if(my_debug()) call out(unit)
#endif
  end subroutine

  integer function weimigz(adj,add,mod_,dat_,SLop,  WCop,FKop,FXop,IGop,op7,op8,op9)
  logical,            intent(in) :: adj,add
  character(len=128), intent(in) :: mod_,dat_
  character(len=128) :: mod_big
  integer,external,optional :: WCop,FKop,FXop,IGop,op7,op8,op9
  real , pointer :: dstep(:),omega(:)
  integer :: i,path,ierr,iz,iw,ith,iws,izs
  character(len=60) :: string
  complex, pointer :: image(:,:,:,:,:)
  type(sep3d) :: dsep
  logical :: ipath,logic
  type(sep3d) :: rsmall,rbig
#ifdef DEBUG
    character(len=128), parameter :: unit='weimigs'
#endif
    interface !!----------------------------------------------
       integer function SLop() result(st)
       end function SLop
    end interface !!------------------------------------------
                                                                                           
#ifdef DEBUG
    if(run%debug) call in(unit)
#endif
  
   call init_sep3d(dat_,dsep,"SCRATCH")
   call init_sep3d(mod_,rsmall,"SCRATCH")
   call store_current_data_space(dsep)
   call store_current_data_tag(dat_)
   call store_current_add(add)
   call store_current_adj(adj)



   if(.not. distribute_image(mod_,rsmall,mod_big,rbig)) &
     call seperr("trouble spreading image")
   call store_current_image_space(rbig)
   call store_current_image_tag(mod_big)
    
   if(.not. my_restart()) then
     if(.not. add)  then
       if(.not. adj) then
          if(.not. zero_data(dat_,dsep)) &
           call seperr("trouble zeroing data")
       else
         if(.not. zero_image(mod_big,rbig)) &
           call seperr("trouble zeroing image")
       end if
     end if
   end if


     
   call assign_wave_space(wfld)  !get memory for wfld
   call assign_image_space(image) !get memory for image



   nullify(dstep,omega)
   if(.not. path_z_begin()) call seperr("trouble initializing path")


   string="modeling"
   if(adj) string="migration"


logic= timer_start("PATH BEGIN")
     if(.not. path_z_op(mode_begin,wfld,image,dstep,omega)) &!Grab the first block
       call seperr("trouble grabbing first path block")

logic= timer_stop("PATH BEGIN")
     
     do while(associated(omega))  !{Go until we hit an error or are finished
        call store_current_s_r(.true.)
        ierr=SLop()
        do iw=1,size(omega) !{ loop over the frequencies in the block
          call store_current_omega(omega(iw))
          ith=1
          do iz=1,size(dstep) !{ loop over the the depths in this block
            call store_current_dstep(dstep(iz))
logic= timer_start("W2R")
            if(adj) call weiimgW2R(iz,iw,ith,image,wfld,IGop)
logic= timer_start("W2R")
logic= timer_start("KXMIG")
            call weikxmig(wfld,iw,iz,ith,WCop,FKop,FXop)
logic= timer_stop("KXMIG")
logic= timer_start("R2W")
            if(.not. adj) call weiimgR2W(iz,iw,ith,image,wfld,IGop)
logic= timer_stop("R2W")
            call print_progress(string)
          end do ! }
        end do! }
logic= timer_start("PATH END")
        if(.not.path_z_op(mode_end,wfld,image,dstep,omega)) &
            call seperr("Trouble finishing block")
logic= timer_stop("PATH END")
logic= timer_start("PATH BEGIN")
       if(.not.path_z_op(mode_begin,wfld,image,dstep,omega)) then
          call seperr("Trouble grabbing next block")
       end if
logic= timer_stop("PATH BEGIN")
     end do !}

    if(.not. combine_image(mod_big,rbig,mod_,rsmall)) &
     call seperr("trouble combining image files")
    
call sep_mpi_stop()
if(sep_thread_num()==0) call print_timers()
    weimigz=ok
#ifdef DEBUG
    if(run%debug) call out(unit)
#endif
  end function
end module wei_mig_base
