!!$=head NAME
!!$
!!$ wei - wave equation migration
!!$
!!$=head1 STATUS PARAMETERS
!!$
!!$=over 4
!!$
!!$=item verb  -  logical  [.false.]  Whether or not to be verbose
!!$
!!$=item debug -  logical  [.false.]  Whether or not to be verbose
!!$
!!$=item depar -  logical  [.false.]  Whether or not to print processing inside 
!!$
!!$                               parallel loops
!!$
!!$=back
!!$
!!$=head1 TYPE OF OPERATION
!!$
!!$=over 4
!!$
!!$=item operator- char[5] [none]     What type of processing
!!$
!!$                              '.R_.D','.Rw.D' - migration
!!$
!!$
!!$=item operation- char[10] [dottest...]  What we are actually doing    
!!$
!!$                      dottest...   Do a dottest
!!$
!!$                      migration.   Do migration
!!$
!!$                      modeling..   Do modeling
!!$
!!$                      inversion.   Do inversion
!!$
!!$=item nfk       - integer              Number of fk operations (e.g. number
!!$
!!$                                   of reference velocities)
!!$
!!$=item adj       - logical [.false.]    Whether or not perform adjoint 
!!$
!!$=item add       - logical [.false.]    Whether or not to add to current model
!!$
!!$=back
!!$
!!$
!!$=head1 MIGRATION SCHEME
!!$
!!$=over 4
!!$
!!$=item distribute - char*
!!$
!!$  [reg] Distribution type 
!!$
!!$    reg -regular sampling
!!$
!!$    vs -variable sampling as a function of slowness
!!$
!!$    vsf-variable sampling as a function of slowness and frequency
!!$
!!$=item nz_mem  - integer (if vs/vsf)
!!$
!!$  [0] Number of z steps to keep in memory (more reduces IO
!!$
!!$      cost of lagrange interpolation, but takes more memory)
!!$
!!$=item nlangrange_pts - integer (if vs/vsf)
!!$
!!$  [5] Number of points to use in langrage interpolation of depth
!!$
!!$=item min_zsamp - real  (if vs/vsf)
!!$
!!$  [dz/3.] Minimum step size in depth
!!$
!!$=item max_zsamp - real  (if vs/vsf)
!!$
!!$  [10000] Maximum step size in depth
!!$
!!$=item cycles - integer  (if vs/vsf)
!!$
!!$  [100000] Cycles of the wavelet before dropping (because of ...)
!!$
!!$=item samp_per - integer  (if vs/vsf)
!!$
!!$  [8] Needed cycles of the wavelet for a given frequency
!!$
!!$=item base_freq -  real (vs)
!!$
!!$  [40.] Frequency to use when calculating needed sampling in vs case
!!$
!!$=item io_cost  - real (vs/vsf)
!!$
!!$  [20.] The cost (in units of number of steps) to do reading/writing
!!$
!!$   for given frequency block
!!$
!!$=item outer_freq - logical   [] Whether or not outer loop is frequency
!!$
!!$                        If nzs > 10  then defaults to .true.
!!$
!!$                        If nzs < 10  then defaults to .false.
!!$
!!$=back
!!$
!!$=head1 BLOCKING PARAMETERS
!!$
!!$=over 4
!!$
!!$=item nzs       - integer              Block size in z
!!$
!!$=item nws       - integer              Block size in w
!!$
!!$=back
!!$
!!$=head1 OTHER PARAMETERS`
!!$
!!$=over 4
!!$
!!$=item image_real  - logical
!!$
!!$ [.false.]  Whether or not image space is real data 
!!$
!!$=item velocity  - logical
!!$
!!$ [.false.]  Whether or not  S is actually velocity
!!$
!!$
!!$=item local_datapath - char*           If running mpi where to put
!!$
!!$                                   local data, image, and slowness
!!$=back
!!$
!!$=head1 LIBRARY
!!$
!!$ sepweif90
!!$
!!$=head1 SEE ALSO
!!$
!!$ L<weimig>
!!$
!!$=cut
!!$
!!$
module wei_process
  use sep
  use wei_util
  use mpi_sep
  use wei_type

  !---------------
  implicit none

  type  z_step
    real, pointer ::  dz(:)
    real, pointer ::  zloc(:)
    integer,pointer :: status(:)
    integer, pointer :: iloc(:,:)
    real, pointer :: frac(:,:)
  end type

  !----------------------------------------------------------------
  type weprocess
     character(len=128) :: operator
     character(len=128) :: operation,file_out

     logical            :: adj, add
     integer            :: opsign
     !! simultate conjg(phase) to avoid 'if'
     !! opsign=goB2T (e.g.  modeling)
     !! opsign=goT2B (e.g. migration)

     logical            :: verb    !! verbose
     logical            :: debug   !! Debug
     logical            :: depar   !! debug in the parallel regions
     logical            :: outer_freq !! wheter or not outer loop is frequency

     character(len=128) :: mod_
     character(len=128) :: dat_
     character(len=4)   :: tag_num

     character(len=128) ::  VS,VG, D, R, W,U !! file tags
     character(len=128) :: dS,   dR

     integer            :: nzs     !! number of depth     slices in a block
     integer            :: nws     !! number of frequency slices in a block

     integer            :: nzb     !! number of depth     blocks
     integer            :: nwb     !! local  number of frequency blocks
     integer            :: nwp     !! global  number of frequency block

     integer            :: nth     !! number of threads
     integer            :: nfk     !! number of f-k ops
     integer            :: mpi_ithread !! mpi thread number
     integer            :: mpi_nthread !! number of mpi threads
     integer            :: omp_nthread !! number of omp threads
     integer,pointer    :: iw_thread(:) !! which mpi thread  each
     ! iw_block belongs to
     integer,pointer    :: iw_first(:)  ! First actual iw in each block
     type(z_step),pointer:: z_samp(:)   ! Sampling of local frequencies
     integer            :: nw_local     !! number of iw owned by this process
     integer            :: distribute_pattern ! 0-reg,1-vel,2-vel+freq,
                                             ! 3-vel+freq+attenuation

     real,pointer       ::  dzb(:) !z sampling for the current block
     integer,pointer    ::  ihxb_beg(:) !beg loop for cur z block for hx axis
     integer,pointer    ::  ihyb_beg(:) !beg loop for cur z block for hy axis
     integer,pointer    ::  ihxb_end(:) !end loop for cur z block for hx axis
     integer,pointer    ::  ihyb_end(:) !end loop for cur z block for hy axis
     real,pointer       ::  max_slow(:) !maximum slowness for each z step
     complex,pointer    :: v_i_buf(:,:,:,:,:) !variable image buffer
     logical, pointer   :: filled(:) !whether or not output iz has been finished
     integer, pointer   :: disk_buf(:)!iz disk buffer element used
     integer, pointer   :: mem_buf(:)!iz memory buffer element used
     integer, pointer   :: write_freq(:) ! what block to start writing on
     integer            :: nmem_buf,ndisk_buf
     logical            :: use_disk_buf ! using disk buf
     logical            :: conv_wave,restart
     

     real               :: pi

     type(axis)         :: amx,amy !! midpoint

     logical            :: pCIG
     type(axis)         :: ahx,ahy !! offset
     type(axis)         :: apx,apy !! offset ray parameter

     type(axis)         :: az,bz   !! depth
     type(axis)         :: aw,bw   !! frequency
     type(weslow)   :: slow
     integer :: io_unit
  end type weprocess
  !----------------------------------------------------------------

  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!
  type(weprocess), public :: pro
  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!

  integer, public :: OK=1
  integer, public :: weiB2T=+1,weiT2B=-1
  integer, public :: goB2T=+1,goT2B=-1
  !! go Top2Bottom or Bottom2Top


  integer, external :: getch
  integer, external :: mpi_sep_tag_sum
  integer, external :: mpi_sep_tag_bcast
  integer, external :: mpi_sep_tag_combine
  integer, external :: mpi_sep_tag_distribute

#ifdef SEP_OMP
  integer, external :: omp_get_thread_num
  integer, external :: omp_get_num_threads
#endif
contains
  !----------------------------------------------------------------
  subroutine process_init()
    character(len=128) :: operator =""
    character(len=128) :: operation=""
    logical            :: adj,add
    integer            :: ierr
#ifdef DEBUG
    character(len=128), parameter :: unit='process_init'
#endif

    call from_param('verb', pro%verb, .false.)
    call from_param('debug',pro%debug,.false.)
    call from_param('restart',pro%restart,.false.)
    call from_param('depar',pro%depar,.false.)

    if(      pro%debug) pro%verb =.true.
    if(.not. pro%debug) pro%depar=.false.

#ifdef DEBUG
    if(pro%debug) call in(unit)
#endif

    call from_param('operator',operator,'.....')
    pro%operator = operator

    call from_param('adj', pro%adj, .false.)
    call from_param('add', pro%add, .false.)

    call from_param('operation',operation,'dottest...')
    pro%operation = operation
    !! inversion begins with adjoint, too
    if(operation(1:10) .eq. 'inversion.') then
       pro%adj=.true. 
    end if

    pro%mod_=''
    pro%dat_=''

    ierr=getch("S_S","d",ierr)+getch("S_G","d",ierr)

    if(2==ierr) then
        pro%VS='S_S'   !!       slowness
        pro%VG='S_G'   !!       slowness
        pro%conv_wave=.true.
    else if( 1==ierr) then
        call erexit("Must provide S_S and S_G")
    else
        pro%VS='S'
        pro%VG='S'
        pro%conv_wave=.false.
    end if


   call from_param("io_file",pro%file_out,"NONE")
   if(pro%file_out(1:4)=="NONE") then
     pro%io_unit=0
   else 
     pro%io_unit=6
     open(pro%io_unit,file=pro%file_out)
   end if


    pro%dS='dS' !! delta slowness
    pro%D='D'   !!       data
    pro%U='U'   !!       downward continued wavefield
    pro%R='R'   !!       image
    pro%dR='dR' !! delta image
    pro%W='W'   !!       wavefield

    call from_param("nzs",pro%nzs,1); pro%nzb=0;
    call from_param("nws",pro%nws,1); pro%nwb=0;
    allocate(pro%dzb(pro%nzs))

    !!  Decide how to set up outer loops
    pro%outer_freq=.false. ; if(pro%nzs > 10) pro%outer_freq=.true.
    call from_param("outer_freq",pro%outer_freq,pro%outer_freq)

    !! read in axes info from the parameter file
    call process_parameters()

    !! figure out the max number of threads
    call process_parallel()
    if(pro%nth>1) pro%depar = .false.

    pro%pi=acos(-1.)

    !! number of reference velocities
    call from_param("nfk",pro%nfk,1)

    !! CIG in offset or offset ray parameter
    call from_param("pCIG",pro%pCIG,.false.)

    if(pro%verb) call process_report()
#ifdef DEBUG
    if(pro%debug) call out(unit)
#endif
  end subroutine process_init
  !----------------------------------------------------------------
  subroutine process_parameters()

    pro%amx%l='mx'; call pinaxis(pro%amx)
    pro%amy%l='my'; call pinaxis(pro%amy)

    pro%ahx%l='hx'; call pinaxis(pro%ahx)
    pro%ahy%l='hy'; call pinaxis(pro%ahy)

    pro%apx%l='px'; call pinaxis(pro%apx)
    pro%apy%l='py'; call pinaxis(pro%apy)

    pro%aw%l ='w_'; call pinaxis(pro%aw )
    pro%az%l ='z_'; call pinaxis(pro%az )

    pro%bz = pro%az; pro%bz%n=pro%nzs
    pro%bw = pro%aw; pro%bw%n=pro%nws
  end subroutine process_parameters
  !----------------------------------------------------------------
  subroutine process_reset(adj,add,mod_,dat_)
    logical,            intent(in) :: adj,add
    character(len=128), intent(in) :: mod_,dat_

    pro%adj  = adj
    pro%add  = add
    if(pro%restart) pro%add=.true.
    pro%mod_ = mod_
    pro%dat_ = dat_

  end subroutine process_reset
  !----------------------------------------------------------------  
  subroutine process_report()
    call weiline()
    write(pro%io_unit,*)                    "        *       model=",pro%mod_
    write(pro%io_unit,*)                    "        *        data=",pro%dat_
    write(pro%io_unit,*)                    "        *    operator=",pro%operator (1: 5)
    write(pro%io_unit,*)                    "        *   operation=",pro%operation(1:10)
    call weiline()
    write(pro%io_unit,*)                    "        *         adj=",pro%adj
    write(pro%io_unit,*)                    "        *         add=",pro%add
    write(pro%io_unit,*)                    "        *        pCIG=",pro%pCIG
    call weiline()
    write(pro%io_unit,*)                    "        *        verb=",pro%verb
    write(pro%io_unit,*)                    "        *       debug=",pro%debug
    write(pro%io_unit,*)                    "        *       depar=",pro%depar
    call weiline()
    write(pro%io_unit,'(a,i5)')             "         *         nth=",pro%nth
    write(pro%io_unit,'(a,i5)')             "         *         nzs=",pro%nzs
    write(pro%io_unit,'(a,i5)')             "         *         nws=",pro%nws
    write(pro%io_unit,'(a,i5)')             "         *         nzb=",pro%nzb
    write(pro%io_unit,'(a,i5)')             "         *         nwb=",pro%nwb
    write(pro%io_unit,'(a,i5)')             "         *         nfk=",pro%nfk
    write(pro%io_unit,'(a,i5)')             "         *        nmpi=",pro%mpi_nthread
    call weiline()
    write(pro%io_unit,'(a,i4,f13.5,f13.5)') "         *          mx= ",pro%amx%n,  pro%amx%o,  pro%amx%d
    write(pro%io_unit,'(a,i4,f13.5,f13.5)') "         *          my= ",pro%amy%n,  pro%amy%o,  pro%amy%d
    write(pro%io_unit,'(a,i4,f13.5,f13.5)') "         *          hx= ",pro%ahx%n,  pro%ahx%o,  pro%ahx%d
    write(pro%io_unit,'(a,i4,f13.5,f13.5)') "         *          hy= ",pro%ahy%n,  pro%ahy%o,  pro%ahy%d
    write(pro%io_unit,'(a,i4,f13.5,f13.5)') "         *          px= ",pro%apx%n,  pro%apx%o,  pro%apx%d
    write(pro%io_unit,'(a,i4,f13.5,f13.5)') "         *          py= ",pro%apy%n,  pro%apy%o,  pro%apy%d
    call weiline()
  end subroutine process_report
  !----------------------------------------------------------------
  subroutine process_follow()

    call weiline()
    if(pro%adj) then
       if(pro%add) then
          write(pro%io_unit,*)'        *   ','[',trim(pro%operator),']  ',trim(pro%mod_),'<+<',trim(pro%dat_)
       else
          write(pro%io_unit,*)'        *   ','[',trim(pro%operator),']  ',trim(pro%mod_),'< <',trim(pro%dat_)
       end if
    else
       if(pro%add) then
          write(pro%io_unit,*)'        *   ','[',trim(pro%operator),']  ',trim(pro%mod_),'>+>',trim(pro%dat_)
       else
          write(pro%io_unit,*)'        *   ','[',trim(pro%operator),']  ',trim(pro%mod_),'> >',trim(pro%dat_)
       end if
    end if
    call weiline()
  end subroutine process_follow
  !----------------------------------------------------------------
  subroutine process_parallel()
    integer :: ierr,ilen
    character(len=1024) :: temp_ch,my_path,machine
#ifdef DEBUG
    character(len=128), parameter :: unit='process_parallel'
    if(pro%debug) call in(unit)
#endif

#ifdef SEP_MPI
    call MPI_Get_processor_name(machine,ilen,ierr)
    call MPI_COMM_SIZE(MPI_COMM_WORLD, pro%mpi_nthread,ierr)
    call MPI_COMM_RANK(MPI_COMM_WORLD, pro%mpi_ithread,ierr)
#else
    pro%mpi_nthread=1
    pro%mpi_ithread=0
#endif

    !TURN OFF VERB AND DEBUG FOR ALL BUT ZERO MPI THREAD
    if(pro%mpi_ithread /= 0) then
       pro%verb=.false.
       pro%debug=.false.
    end if

#ifdef SEP_MPI
    if(pro%mpi_nthread>1) then
       my_path="                                             &
       &                                                                       &
       &                    "
       if(0==getch("local_datapath","s",my_path)) &
       my_path="/net/"//machine(1:ilen)//"/scr1/temp/"


       if(pro%mpi_ithread<10) then
          write(pro%tag_num,"(a,i1)") "_",pro%mpi_ithread
       else if(pro%mpi_ithread <100) then
          write(pro%tag_num,"(a,i2)") "_",pro%mpi_ithread
       else
          write(pro%tag_num,"(a,i3)") "_",pro%mpi_ithread
       end if
       write(temp_ch,"(a,a,a,a)")"datapath=",trim(my_path),trim(pro%tag_num),&
       "DATA_"
       call getch_add_string(temp_ch)
    end if
#endif

    !$OMP PARALLEL
#ifdef
    pro%nth = omp_get_num_threads()
    pro%omp_nthread=omp_get_num_threads()
#else
    pro%nth=1
    pro%omp_nthread=1
#endif
    !$OMP END PARALLEL

    if(pro%nws/pro%nth * pro%nth /=pro%nws) then
       write(0,*) "nws=",pro%nws
       write(0,*) "nth=",pro%nth
       call seperr('nws /= k * nth')
    end if

    pro%nwp=pro%aw%n/pro%nws
    allocate(pro%iw_thread(pro%nwp))

#ifdef DEBUG
    if(pro%debug) call out(unit)
#endif
  end subroutine process_parallel
  !----------------------------------------------------------------
  subroutine wei_followup(izb,iwb,n_z)
    integer, intent(in) ::izb,iwb,n_z
    optional :: n_z
    integer :: i
    integer :: allblocks, doneblocks
    real    :: completed

    allblocks = pro%nzb*pro%nwb

    select case(pro%operator(1:5))
    case('.R_.D','.Rw.D')

       if(pro%adj) then !! migration
          if(pro%outer_freq) then
             doneblocks= (iwb-1)*pro%nzb + izb -1
          else
             doneblocks= (izb-1)*pro%nwb + iwb -1
          end if
       else             !! modeling
          if(pro%outer_freq) then
             doneblocks= (iwb-1)*pro%nzb + pro%nzb-izb
          else
             doneblocks= (pro%nzb-izb)*pro%nwb + iwb -1
          end if
       end if

    case('dS_dR','dSwdR')

       if(pro%adj) then !! backprojection
          if(pro%outer_freq) then
             doneblocks= (iwb-1)*pro%nzb + pro%nzb-izb
          else
             doneblocks= (pro%nzb-izb)*pro%nwb + iwb -1
          end if
       else             !! scattering
          if(pro%outer_freq) then
             doneblocks= (iwb-1)*pro%nzb + izb -1
          else
             doneblocks= (izb-1)*pro%nwb + iwb -1
         end if
       end if

    end select


    completed = 100.0 * doneblocks/allblocks !! percent completed

    write(pro%io_unit,'( "         *&
     & done ",f5.2,"%","    izb=",i4,"/",i4,"     iwb=",i4,"/",i4)') &
     completed,izb,pro%nzb,iwb,pro%nwb
  end subroutine wei_followup
  !----------------------------------------------------------------
  subroutine process_radians()
    pro%aw%o = 2*pro%pi * pro%aw%o
    pro%aw%d = 2*pro%pi * pro%aw%d
    pro%bw%o = 2*pro%pi * pro%bw%o
    pro%bw%d = 2*pro%pi * pro%bw%d
  end subroutine process_radians

  subroutine check_hp()
    if(pro%apx%n /= pro%ahx%n) then
       call seperr('apx.n /= ahx.n')
    end if
    if(pro%apy%n /= pro%ahy%n) then
       call seperr('apy.n /= ahy.n')
    end if
  end subroutine check_hp
  !----------------------------------------------------------------
end module wei_process
