# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
module SPmig_mod !! WEI kerneled migration, Shot Profile style

  !! mosteverything
  use down_mod
  !!file stuff
  use S_mod
  use R_mod
  use R_write_mod
  use D_mod
  use W_mod
  !!transforms
  use kmap_mod
  use xmap_mod
  use fft_mod
  !! extrapolation
  use kxmig_mod
  !! imaging condition
  use image_mod
  use taper_mod
  use xmap_mod
  implicit none

  !! files
  type(dat) ,public :: D,D_d,HD,HD_d
  type(imag),public :: R,M,N_u,N_d
  type(twfld),public:: W,W_d
  !! timers
  integer:: time_array(8),jwork,jstep,donez,donew
  real::startd,work(10),steps
  logical:: inim=.false.
contains
!----------------------------------------------------------------
  subroutine SPmig_init(SLin,WCin,FKin,FXin,IGin)
    integer :: st,er
    real::a



    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 !!---------------------

    startd=0.
    !! initialize the main running variables
    if(.not. run_init()) &
    call erexit("trouble initializing run_init")
    if(.not. down_init()) &
    call erexit("trouble initializing down_init")




    !!bring in the files
    call file_io()

    !! convert frequency to radians
    call process_radians()

    !! init FFT
    write(run%io,*)'        ************************************************'
    write(run%io,*)"        *    Generating FFT plans "
    if(down%sourcefn.or.down%passive) then
       call fft_init_s(D_d%bin)
    end if
    if(down%gathers)then
       call fft_init_r(D%bin)
    end if
    !!linear maps of k values
    call kmap_init()

! !! init imaging
! call imaging_init()

    !! allocate slowness arrays for FK/FX ops
    call slow_allocate()

    !! init operators
    st=SLin() !! FK slownesses init
    st=WCin() !! WCop init
    st=FKin() !! FKop init
    st=FXin() !! FXop init
    st=IGin() !! Imaging condition init

    !! reset tag variables
    D_d%tag= trim(down%dat_d_)
    D%tag = trim(down%dat_ )
    R%tag = trim(down%mod_ )
    M%tag = trim(down%m_ )
    N_d%tag= trim(down%n_d_ )
    N_u%tag= trim(down%n_u_ )
    W_d%tag= trim(down%w_d_ )
    W%tag = trim(down%w_ )
    HD%tag = trim(down%h_d_ )
    HD_d%tag=trim(down%h_d_d_)

    steps=(down%az%n-down%awtr%j+down%awtr%n)*down%aw%n
    do st=1,size(work)
       work(st)=nint(steps*0.1*real(st))
    end do

    write(run%io,'(a,f6.2,a)')"         *    Allocated ",down%mem/1024./1024," Mbytes of arrays. "
    write(run%io,*)'        ************************************************'




  end subroutine SPmig_init
  !----------------------------------------------------------------
  integer function SPmig(adj,add,SLop,WCop,FKop,FXop,IGop,MIGop,NDIGop)
    logical,intent(in):: adj,add
    integer,external :: WCop,FKop,FXop,IGop,MIGop,NDIGop
    integer :: izb,iwb,isx,is
    integer :: st
    real::startSHOT,endSHOT



    interface !!----------------------------------------------
       integer function SLop() result(st)
       end function SLop
    end interface !!------------------------------------------



    if(run%verb)write(run%io,*) '        * Start parallel'
    call sep_begin_prog()
    if(run%verb)write(run%io,*) '        * '
    down%shotnum=0
    jstep=0;jwork=1
    shot: do is=1,down%evnt%n
       isx=is+down%evnt%b-1
       call start_timer(startSHOT)

       call line()
       write(run%io,'(a,i5,a,i5)')"         *  Shot number",&
         down%shotnum+1,' of',down%evnt%n

       call data_position(D,isx)
       call source_pos(D_d,isx)
       call imagespace()
       if(down%headD) call data_header_read(HD)

       SPmig=SPmigof(isx,SLop,WCop,FKop,FXop,IGop,MIGop,NDIGop)

       call end_timer(endSHOT,startSHOT)
       if(run%verb)write(run%io,'(a,i4,a,f6.2,a)') &
       '         *       Finished migrating shot',down%shotnum +1,&
       ' in',endSHOT/60.,'(m)'

       down%shotnum=down%shotnum+1
       if (down%sourcefn.and.D_d%asx%n.le.1) then
          st=sseek(D_d%tag,0,0) !! rewind to beginning of source
       end if
       inim=.false.
    end do shot

    if(down%passive)call aux_unlink('D_d')
    call writeparams()
    call message("Done and Done")



    SPmig=OK
  end function SPmig
!--------------------------------------------------------------------
  integer function SPmigof(isx,SLop,WCop,FKop,FXop,IGop,MIGop,NDIGop)
    integer,external :: WCop,FKop,FXop,IGop,MIGop,NDIGop
    integer :: izb,iwb,position
    integer :: st,ierr,isx,nz
    real:: dz
    real::startREAD,endREAD,startRITE,endRITE
    logical::readvel=.true.



    interface !!----------------------------------------------
       integer function SLop() result(st)
       end function SLop
    end interface !!------------------------------------------

    donew=0;iwb=1;jstep=0;jwork=1
    down%az%b=1
    if(down%adj) then !! MIGRATION




       call start_timer(startd)
       freqblock:do while (donew.lt.down%aw%n)
          down%bw%o = down%aw%o + (iwb-1)*down%bw%d
          if(run%verb)write(run%io,'(a,i3,a,f5.2)') &
          '         * Frequency block ',iwb,' of ',real(down%aw%n)/real(down%wblock)
          call read_data_source(iwb)

          donez=0;izb=1
          if(down%watercolumn) then
             if(run%verb)write(run%io,*)'        ************************************************'
             if(run%verb)write(run%io,'(a,i3,a)')'       **  Datuming through watercolumn in',down%awtr%n,' steps'

             call read_slowness(1) !! read a slowness block
             st = SLop() !! setup sS, sR, s(z) for working space

             !!save parameters
             dz=down%az%d
             nz=down%nzs

             !!modify parameters
             down%az%d=down%awtr%d
             down%nzs=down%awtr%n

             !! do the datuming deed
             call migrate(izb,iwb,WCop,FKop,FXop,IGop,MIGOP,NDIGop)

             !! restore parameters
             down%nzs=nz
             down%az%d=dz

             !izb and izs to continue subrtn migrate()
             izb=1+down%awtr%j/down%nzs
             if (izb.eq.1) readvel=.false.
             down%az%b=1.5+down%awtr%j -(izb-1)*down%nzs-1
             down%awtr%b=izb
             if(run%verb)write(run%io,'(a,i3,a)')'       **  Finished datuming '
             if(run%verb)write(run%io,*)'        ************************************************'
             donez=(izb-1)*down%nzs
          end if

          depthblock:do while (donez.lt.down%az%n)
          startREAD=0.;endREAD=0.;startRITE=0.;endRITE=0.
       call start_timer(startREAD)
             call read_image(izb)
       call end_timer(endREAD,startREAD)
             if(readvel) then
       call start_timer(startREAD)
                call read_slowness(izb) !! only if you didnt for datuming thru wtrcolumn
       call end_timer(endREAD,startREAD)
                st = SLop() !! setup sS, sR, s(z) for working space
             end if

             call migrate(izb,iwb,WCop,FKop,FXop,IGop,MIGop,NDIGop)
       call start_timer(startRITE)
             call write_image(izb)
       call end_timer(endRITE,startRITE)
             donez=donez+down%nzs
             down%nzs=min(down%nzs,down%az%n-donez)
             izb=izb+1
             down%az%b=1
             readvel=.true.
             down%awtr%b=1
             if(run%verb)write(run%io,'(a,f9.3,a,f9.3,a)') &
               '         *       IO time read:',endREAD,'(s). time write:',endRITE,'(s).'
             if(run%verb)write(run%io,'(a,f9.3,a)') &
               '         *       IO time total:',(endREAD+ endRITE),'(s).'
          end do depthblock
          down%nzs=down%zblock
          st=1+down%az%n-(izb-2)*down%zblock
          call writewavefields(st,izb-1,iwb) !! last depth
          donew=donew+down%nws
          down%nws=min(down%nws,down%aw%n-donew)
          iwb=iwb+1
       end do freqblock
       down%nws=down%wblock



    else !! MODELING



       write(run%io,*) "There is no such beast as the Adjoint of ShotMig"



    end if



    SPmigof=OK
  end function SPmigof
!---------------------------------------------------------------
  subroutine migrate( izb,iwb,WCop,FKop,FXop,IGop,MIGOP,NDIGop)
    integer,intent(in)::izb,iwb
    integer,external :: WCop,FKop,FXop,IGop,MIGOP,NDIGop
    integer :: izs,iws,a,currentz,ith=1
    real :: startIC,finishIC,startEX,finishEX




    startIC=0.;finishIC=0.;startEX=0.;finishEX=0.
    if (down%image_freqs) a=sseek(R%tag,0,0)
    finishIC=0.
    finishEX=0.
    write(run%io,*) "        *       Migration begins"

    zloop: do izs=down%az%b,down%nzs !LOOP OVER z STEPS
       call writewavefields(izs,izb,iwb)
!#ifdef 1
!!$OMP PARALLEL PRIVATE(ith) FIRSTPRIVATE(finishIC,startIC,startEX,finishEX)
! ith=1+omp_get_thread_num()
!!$OMP DO PRIVATE(iws)
!#endif
! ws: do iws=1,down%nws !LOOP OVER w STEPS

          call start_timer(startIC)

          call imagingcondition(izb,izs,iws,ith,IGop,MIGOP,NDIGop)

          call end_timer(finishIC,startIC)
          call start_timer(startEX)

          call kxmig(izs,iws,ith,D%bin,D_d%bin,WCop,FKop,FXop)


          call end_timer(finishEX,startEX)

! end do ws
!!$OMP END DO
!!$OMP END PARALLEL
       jstep=jstep+down%nws
       if (work(jwork).le.jstep) call depth_report(jwork,jstep)
    end do zloop

    finishIC=finishIC/60.
    finishEX=finishEX/60.
    write(run%io,'(8x,a,f6.2,a)')' *       Time imaging       :', &
    finishIC,' (m)'
    write(run%io,'(8x,a,f6.2,a)')' *       Time extrapolating :', &
         finishEX,' (m)'





  end subroutine migrate
!----------------------------------------------------------------
  subroutine read_data_source(iwb)
    integer:: iwb,position

    if(.not.down%headD) then !! read data
       position=down%shotnum*down%bw%n + (iwb-1)

       if(down%gathers) then
          if(.not. data_read(D,position)) &
            call seperr("trouble reading D data") !! read data at block top
       end if

       if(down%sourcefn) then !! read the source function
          if(.not. source_read(D_d,iwb-1)) & !!read the source file (conjugate at same time)
          call seperr("trouble reading from D_d.")
       end if
    else !! read data with header locations
       position=(iwb-1)

       if(.not. data_sort_read(D,D_d,position,HD)) &
         call seperr("trouble reading from D or HD.") !! read data at block top

       if(down%sourcefn) then
          if(.not. source_read(D_d,iwb-1)) & !!read the source file (conjugate at same time)
            call seperr("trouble reading from D_d.")
       end if

    end if

    if(down%passive) then!!use conj of data as source wavefield
       call data_copy_conjg(D,D_d)
    end if

  end subroutine read_data_source
!----------------------------------------------------------------
  subroutine writewavefields(izs,izb,iwb)
    integer, intent(in) :: izs,izb,iwb
    integer::z

    z=(izb-1)*down%zblock+izs
    !! write wavefields if you want them
    if (down%wIOu)then
       if(z.ge.down%firstzW) then !! limit W file size when min. target depth known
          if (.not. wavefield_write(W ,D%bin ,iwb,izb,izs-down%firstzW+1)) &
          call seperr('wavefield_write(W ==badness')
       end if
    end if
    if (down%wIOd)then
       if(z.ge.down%firstzW) then !! limit W file size when min. target depth known
          if (.not. wavefield_write(W_d,D_d%bin,iwb,izb,izs-down%firstzW+1)) &
          call seperr('wavefield_write(W_d ==badness')
       end if
    endif
  end subroutine writewavefields
!----------------------------------------------------------------
  subroutine imagingcondition(izb,izs,iws,ith,IGop,MIGOP,NDIGop)
    integer, intent(in) :: izb,izs,iws,ith
    integer,external :: IGop,MIGOP,NDIGop
    integer::idepth,iw,dd,zz
    !! taper the edges of the wavefields

    if(down%taperme) then
!!$OMP PARALLEL DO PRIVATE( iw)
! do iw=1,down%nws
          if(down%sourcefn.or.down%passive)call taperwfld(D_d,iw,ith)
! end do
! do iw=1,down%nws
          if(down%gathers) call taperwfld(D,iw,ith)
! end do
!!$OMP END PARALLEL DO
    end if

    if(inim) then
       zz=1
    else
       zz=down%firstz
    end if
    dd=(izb-1)*down%nzs+izs
    if(dd.lt.down%firstz) return

    idepth=izs-zz+1 + (down%awtr%b-1)*down%nzs

    !! all the various imaging conditions
    if(down%normd)call SPimg(idepth,iws,ith,D%bin,D_d%bin,N_d,NDIGop)!DD*
    if(down%mult) call SPimg(idepth,iws,ith,D%bin,D_d%bin,M, MIGop)!UU
    if(down%image)call SPimg(idepth,iws,ith,D%bin,D_d%bin,R, IGop) !UD*
  end subroutine imagingcondition
!----------------------------------------------------------------
  subroutine addimages(izs)
    integer, intent(in) :: izs
    integer::idepth
    if(izs.lt.down%firstz) return
    idepth=izs-down%firstz+1
    if(down%normd)call addimagethrds(idepth,N_d)
    if(down%mult) call addimagethrds(idepth,M)
    if(down%image)call addimagethrds(idepth,R)
  end subroutine addimages
!----------------------------------------------------------------
  subroutine writeparams()
    if(down%normd)call param_write(N_d%tag)
    if(down%image)call param_write(R%tag)
    if(down%mult) call param_write(M%tag)
    if(down%wIOu) call param_write(W%tag)
    if(down%wIOd) call param_write(W_d%tag)
  end subroutine writeparams
!----------------------------------------------------------------
  subroutine depth_report(jw,js)
    integer :: jw,js
    real:: stopd
    jw=jw+1
    stopd=0.
    call end_timer(stopd,startd)
    stopd=stopd/60.
    if(run%verb) then
       write(run%io,'(a,i5,a,i5,a,i2,a,f6.2,a)')'         *  Finished ',&
         js/1000,'k of',int(steps/1000.),'k w*z steps (',int(100.*js/steps),'%):',stopd,'(m)'
       call start_timer(startd)
    end if
  end subroutine depth_report
!----------------------------------------------------------------
  subroutine file_io()
    integer:: err
    logical:: output
    character(len=1024):: present
    output=F
    !! . . save the file tags in case they change
    down%dat_d_ = down%D_d !"/homes/sep/claudio/research.ln/Prst_ExpRefl/Src_Brad/SPmig.F90" DATA IS FREQUENCY SLICES
    down%dat_ = down%D !DATA IS FREQUENCY SLICES
    down%mod_ = down%R !IMAGE SPACE (REFLECTIVITY)
    down%m_ = down%M !image space SRME model
    down%n_u_ = down%N_u !image space SRME model
    down%n_d_ = down%N_d !image space SRME model
    down%w_d_ = down%W_d !downgoing wfld
    down%w_ = down%W !upcomping wfld
    down%h_d_ = down%HD !headers
    down%h_d_d_ = down%HD_d !headers

    !! initialize objects/files
    if (down%gathers)then
       if(down%datum) then !!only situation that doesnt need source wflds
          if(down%headD) then !! receiver headers
             if(run%verb)write(run%io,*)'        ************************************************'
             if(run%verb)write(run%io,*)'        *  D/= sep77.  Workspace sized from parameters'
             if(.not.data_sort_init(D,down%D,HD,down%HD)) &
               call erexit("trouble initializing D,HD for sort_read ")
          else !! regular cube data init
             if(.not. data_init(D,down%D)) &
               call erexit("trouble initializing D")
          end if
       else !! source and data wflds used
          if(down%headD_d)then !! source headers data init
             if(.not.source_sort_init(HD_d,down%HD_d))&
               call erexit("trouble initializing source headers HD_d")
          end if
          if(down%headD) then !! receiver headers
             write(run%io,*)'        ************************************************'
             write(run%io,*)'        *  D/= sep77.  Workspace sized from parameters'
             if(.not.data_sort_init(D,down%D,HD,down%HD)) &
               call erexit("trouble initializing D,HD for sort_read ")
          else !! regular cube data init
             if(.not. data_init(D,down%D)) &
             call erexit("trouble initializing D")
          end if
          if (down%passive)then
             if(.not. aux_data_init(D_d,down%D_d)) &
             call erexit("trouble initializing passive D_d")
          else
             if (down%sourcefn.and..not.allocated(D_d%bin)) then
                if(.not. source_init(D_d,down%D_d)) &
                call erexit("trouble initializing D_d")
! if(D%asx%d/down%amx%d ) then
! write(run%io,*)"Source locations dont fall exactly on wfld grid (X) "
! if(D%asy%d/down%amy%d ) then
! write(run%io,*)"Source locations dont fall exactly on wfld grid (Y) "
             else
                call seperr("Need source (D_d=file.H) you moron.")
             end if
          end if
       end if
    else !! no data, only have D_d
       if(down%headD_d.and.down%sourcefn)then
          if(.not.source_sort_init(HD_d,down%HD_d))&
          call erexit("trouble initializing source headers HD_d")
          if(.not.workspace_params(D)) call seperr("Badness in workspace_params")
       elseif(down%sourcefn) then
          if(.not.workspace_params(D_d)) call seperr("Badness in workspace_params")
          call smap_init(D_d)
       else
          call seperr("Need source (D_d=file.H) you moron.")
       end if

       if(run%verb)write(run%io,*)'        ************************************************'
       if(run%verb)write(run%io,*)'        *  No gathers.  Workspace sized from parameters'
       if(.not. source_init(D_d,down%D_d)) call erexit("trouble initializing D_d")
    end if


    !! cosine taper for the wavefields
    if(down%taperme) call xytaper_init()

    call image_size_calc() !! image domain coordinate calculations
    down%ax=down%amx;down%ay=down%amy

    if (down%gathers) then
       if (down%image.and.(.not.down%datum)) then
          if((.not.down%sourcefn).and.(.not.down%passive)) call seperr("Need D and D_d for R")
          if(.not. image_init(R,down%R)) &
          call erexit("trouble initializing image R")
          if(run%verb) write(run%io,*) "        *      Image Space     (UD*)          R=",down%image
          output=.true.
       end if
       if (down%mult.and.(.not.down%datum) ) then
          if (.not.output) then
             if(.not. image_init(M,down%M)) &
             call erexit("trouble initializing image M")
             if(run%verb) write(run%io,*) "        *      Image Space SRMP (UU)          M=",down%mult
             output=.true.
          else
             if(.not. aux_image_init(M,R,down%M))call seperr('badness at image_init(M)')
             if(run%verb) write(run%io,*) "        *      Image Space SRMP (UU)          M=",down%mult
          end if
       end if
    end if


    if(down%sourcefn) then
       if (down%normd )then
          if(output) then
             if(down%image) then
                if(.not. aux_image_init(N_d,R,down%N_d))call seperr('badness at bwfld_init(N_d)')
             elseif(down%mult) then
                if(.not. aux_image_init(N_d,M,down%N_d))call seperr('badness at bwfld_init(N_d)')
             end if
          else
             if(.not. image_init(N_d,down%N_d)) &
             call erexit("trouble initializing image N_d")
             output=.true.
          end if
          if(run%verb) write(run%io,*)'        ************************************************'
          if(run%verb) write(run%io,*) "        *      Writing normalization (DD*)    N_d=",down%normd
       end if
    end if

    if (down%gathers) then
       if (down%WIOu )then
          if(down%image_freqs)call seperr("piss off, you cant have everything")
          if(.not. bwfld_init(W,down%W))call seperr('badness at bwfld_init(W)')
          output=.true.
       else
          if(down%datum.and.(.not.output)) call seperr("Source modeling needs output. I suggest W=")
       end if
    endif
    if(down%sourcefn) then
       if ( down%wIOd) then
          if(down%image_freqs.or.down%datum)call seperr("piss off, you cant have everything")
          if(.not. bwfld_init(W_d,down%W_d))call seperr('badness at bwfld_init(W_d)')
          output=.true.
       end if
    end if

    if(.not.output) then
       call seperr("You must have some output idiot.")
    end if

    if(.not. slow_init(S,down%VS,down%VR))&
    call erexit("trouble initializing slowness/velocity")


  end subroutine file_io
!----------------------------------------------------------------
  subroutine read_image(izb)
    integer,intent(in)::izb
    integer:: dd,zz,idepth,n,z2

    dd=(izb)*down%nzs
    if(dd.lt.down%firstz) return
    if(.not.inim) then
       zz=1
       z2=down%firstz
       idepth=0
    else
       zz=down%firstz
       z2=1
       idepth=(izb-1)*down%zblock -zz+1
    end if
    n=min(down%nzs,down%nzs-z2+1+ (down%awtr%b-1)*down%nzs)

    if (down%image) then
       call image_location(R)
       call image_read(R,idepth,n)
    end if

    if(down%mult) then !! IS-SRMP output=UU
       if(down%image)then
          call update_index(M,R)
       else
          call image_location(M)
          call update_index(R,M)
       end if
       call image_read(M,idepth,n)
    end if

    if(down%normd) then !! normalization output=DD*
       if(down%image.or.down%mult)then
          call update_index(N_d,R)
       else
          call image_location(N_d)
          call update_index(R,N_d)
       end if
       call image_read(N_d,idepth,n)
    end if

  end subroutine read_image
!----------------------------------------------------------------
  subroutine write_image(izb)
    integer,intent(in)::izb
    integer:: dd,zz,idepth,n,z2

    dd=(izb)*down%nzs
    if(dd.lt.down%firstz) return

    if(.not.inim) then
       zz=1
       z2=down%firstz
       idepth=0
    else
       zz=down%firstz
       z2=1
       idepth=(izb-1)*down%zblock -zz+1
    end if
    n=min(down%nzs,down%nzs-z2+1+ (down%awtr%b-1)*down%nzs)

    if (.not. down%image_freqs) then
       if(down%image)call image_write(R ,idepth,n)!! migration output=UD*
       if(down%mult )call image_write(M ,idepth,n)!! IS-SRMP output=UU
       if(down%normd)call image_write(N_d,idepth,n)!! normalization output=DD*
    end if
    inim=.true.
  end subroutine write_image
!----------------------------------------------------------------
  subroutine start_timer(time)
    real ::time
    call date_and_time(values=time_array)
    time = time_array (5) * 3600 + time_array (6) * 60 &
         + time_array (7) + 0.001 * time_array (8)
  end subroutine start_timer
!----------------------------------------------------------------
  subroutine end_timer(ftime,stime)
    real::ftime,stime
    call date_and_time(values=time_array)
    ftime = ftime + (-stime+( time_array (5) * 3600 + time_array (6) * 60 &
                  + time_array (7) + 0.001 * time_array (8)))
  end subroutine end_timer
!----------------------------------------------------------------
! subroutine mig_outnull()
! !ZERO OUTPUT IF REQUESTED
!#ifdef DEBUG
! character(len=128), parameter :: unit='mig_outnull'
! if(run%debug) call in(unit)
!#endif
!
! if(down%adj) then !! R<-D
! call image_assign(R,(0.,0.))
! else !! R->D
! call data_assign(D ,(0.,0.))
! end if
!
!#ifdef DEBUG
! if(run%debug) call out(unit)
!#endif
! end subroutine mig_outnull
!----------------------------------------------------------------
end module SPmig_mod
