module D_mod
  use util_mod
  use down_mod
  use external_mod
  use type_mod

  !---------------
  implicit none
  complex,dimension(:),allocatable::source
  complex,dimension(:,:),allocatable::tmp
  integer,dimension(:,:),allocatable::bincount
  type(txmap), public :: smap
  logical,private,save::planemodel
contains
  !--------------------------------------------------------------
  !  Function to initialize a data objects either from
  !    a file or from  parameters
  logical function data_init(X,tag)
    type(dat)        :: X
    character(len=128), intent(in)    :: tag

    data_init=F

    X%tag = tag
    if(isfile(X%tag)) then
       if(.not. data_history_init(X)) then
         write(run%io,*) "trouble initializing data_history"
         return
       end if
    else
       if(.not. data_param_init(X)) then
         write(run%io,*) "trouble initializing data_parameter"
         return
       end if
    end if
    data_init=T
  end function data_init
!----------------------------------------------------------------
  logical function workspace_params(X)
    type(dat):: X
    integer::n
    workspace_params=F
    call from_param("ox" ,down%amx%o,0.)
    call from_param("dx" ,down%amx%d,10.)
    call from_param("nx" ,n,200)
    down%amx%n=next_fast(n)

    call from_param("oy" ,down%amy%o,0.)
    call from_param("dy" ,down%amy%d,10.)
    call from_param("ny" ,n,1)
    down%amy%n=next_fast(n)

    if(down%green.and.(.not.down%headD_d)) then
       if(down%asx%n.eq.1) then
          call from_param("osx",down%asx%o,(down%amx%o+(down%amx%n-1)*down%amx%d)/2)
          call from_param("dsx",down%asx%d,10.)
          call from_param("nsx",down%asx%n,1)
       end if
       if(down%asy%n.eq.1) then
          call from_param("osy",down%asy%o,(down%amy%o+(down%amy%n-1)*down%amy%d)/2)
          call from_param("dsy",down%asy%d,10.)
          call from_param("nsy",down%asy%n,1)
       end if
       X%asx=down%asx
       X%asy=down%asy
       down%evnt%n=down%asx%n*down%asy%n
    else
       down%asx=X%asx
       down%asy=X%asy
    end if
    if(down%offset.and.(.not.down%roll)) then
       down%amx%o= real(smap%s(down%evnt%b))+down%amx%o
       down%amy%o=aimag(smap%s(down%evnt%b))+down%amy%o
    end if

    !! if rollalong, X%amx%o needs to be offset.  Not used if stationary.
    X%amx=down%amx
    X%amy=down%amy

    X%aw=down%aw
    X%tag=down%D
    X%from ='Fake-o-rama'
    X%esize=8
    X%indx=1;                      X%indy=1
    X%finx=X%indx + X%amx%n -1;    X%finy=X%indy + X%amy%n -1
    !! . . geometry check
    if((X%amx%d .lt. 0) .or. (X%amy%d .lt. 0))&
    call seperr("Please refrain from negative d1 and/or d2")
    workspace_params=T
    down%ax=down%amx;    down%ay=down%amy
  end function workspace_params
  !--------------------------------------------------------------
  !
  !  Function to initialize source object from input or 
  !          construct a wavelet
  ! 
  logical function source_init(X,tag)
    type(dat) :: X
    character(len=128), intent(in)    :: tag

    source_init=F

    X%tag = tag
    if(isfile(X%tag)) then
       if(.not. source_history_init(X) ) then
         write(run%io,*) "trouble at source_history_init"
         return
       end if
    else
!       if(.not. source_param_init(X)) then
!         write(run%io,*) "trouble making source function"
!         return
!       end if
       call erexit("Source module not built")
    end if

    source_init=T

  end function source_init

  !----------------------------------------------------------------
  !
  !
  !  Create an auxilary data file
  !
  logical function aux_data_init(X,tag)
    type(dat) :: X
    character(len=128), intent(in)    :: tag

    aux_data_init=F

    X%tag = tag
    if(.not. data_scratch_init(X)) then
       write(run%io,*) "trouble initializing scratch data file"
       return
    end if

    aux_data_init=T

  end function aux_data_init
!----------------------------------------------------------------
!   To use unsorted traces with headerfile locations (not f(z) yet)
!    file HD has indices to write traces into
  logical function data_sort_init(X,tag,Y,htag)!!Y=headers,X=gathers
    type(dat)::X,Y
    character(len=128),intent(in):: tag,htag
    integer::ix,iy,itr,stat
    real::a,mini,maxi,b
#ifdef DEBUG
    character(len=128), parameter :: unit='data_sort_init'
    if(run%debug) call in(unit)
#endif
    data_sort_init=F

    Y%tag=htag
    X%tag=tag
    call from_aux(tag,'esize',X%esize,4)
    if(X%esize.eq.4)call seperr("Real data not allowed puta")
    call hinaxis(tag,down%trc,2); X%amx%l='ntraces'
    call hinaxis(tag,down%evnt,3); X%amy%l='nevents'
    call hinaxis(tag,X%aw ,1); X%aw%l='w_'
    call hinaxis(htag,Y%amy,1); X%amy%l='xylocs'
    call hinaxis(htag,Y%amx,2); X%amx%l='ntraces'
    call hinaxis(htag,Y%asx,3);Y%amy%l='nevents'

    if(Y%amy%n/=2) call seperr("I want a 2D acquisition plane.")
    if(Y%amx%n/=down%trc%n) call seperr("Number traces D/=number headers HD.  No good.")
    if(Y%asx%n/=down%evnt%n.and.down%roll) call seperr("number shots D/= number shots HD.  Lame.")

    X%asx%n=down%evnt%n ;X%asy%n=1
	
    !! first shot
    if(size(smap%s).gt.down%evnt%n)& 
      down%evnt%b=int(down%evnt%o) 
    if(size(smap%s).eq.down%evnt%n)&
       down%evnt%b=1
    if(size(smap%s).lt.down%evnt%n)&
      call seperr("Must have as least as many source headers as shots in D.")

    stat=2*Y%amx%n*4
    if(.not.workspace_params(X)) &
     call seperr("Bad workspace_params in data_sort_init.")
    if(.not.data_setup(X,.false.)) &
     call seperr("Bad data_setup in data_sort_init. ")

    down%asx%o=minval(real(smap%s))
    down%asy%o=minval(aimag(smap%s))

    allocate(tmp(down%nws,down%trc%n))
    allocate(Y%bin(Y%amx%n,1,1))
    allocate(bincount(down%amx%n,down%amy%n))

    tmp=cmplx(0.,0.)
    Y%bin=cmplx(0.,0.)
    bincount=0
    call data_header_read(Y)
    stat=sseek(Y%tag,0,0)

    !! check to see if headers fit into workspace 
    if(down%amx%n.gt.1) then
       mini=minval(real(Y%bin))
       if(mini.lt.1.) &
       write(run%io,*)" ** Min trace location (x) less than ox parameter."
       mini=minval(real(smap%s))
       if(down%amx%o.gt.mini) &
       write(run%io,*) " ** Min shot (x) not on wfld grid (x)."

       maxi=maxval(real(Y%bin))     
       if(maxi.gt.real(down%amx%n))&
       write(run%io,*)" ** Max trace location (x) less than o,n,d (x) parameters."
       maxi=maxval(real(smap%s))
       if(maxi.gt.down%amx%o+(down%amx%n-1)*down%amx%d)&
       write(run%io,*)" ** Max shot (x) not on wfld grid (x)."
    end if
    if(down%amy%n.gt.1) then
       mini=minval(aimag(Y%bin))
       if(mini.lt.1.) &
       write(run%io,*)" ** Min trace location (y) less than oy parameter."
       mini=minval(aimag(smap%s))
       if(down%amy%o.gt.mini) &
       write(run%io,*) " ** Min shot (y) not on wfld grid (y)."

       maxi=minval(aimag(Y%bin))
       if(maxi.gt.real(down%amy%n))&
         write(run%io,*)" ** Max trace location (y) greater than o,n,d (y) parameters."
       maxi=maxval(aimag(smap%s)) 
       if(maxi.gt.down%amy%o+(down%amy%n-1)*down%amy%d)&
         write(run%io,*)" ** Max shot (y) not on wfld grid (y)."
    end if

    X%from = "parameter"
    X%amx = down%amx; X%amy = down%amy; 
    X%asx = down%asx; X%asy = down%asy; 
    X%aw  = down%aw;  
    X%amx%o=down%amx%o!-down%asx%o
    X%amy%o=down%amy%o!-down%asy%o

    a=8*(size(tmp)+size(Y%bin))+4*size(bincount)
    down%mem=down%mem+a
#ifdef MEMCHECK
    if(run%debug) write(run%io,*) 'D.F90 HD%bin, tmp,headers,bincount: allocate',a/1024./1024.,'Mbytes.  Total=',down%mem/1024./1024.
#endif

    Y%from="history"
    X%from='Internally sorted'

    write(run%io,*)'        *     ',trim(X%from),'     Headers in file : ',trim(Y%tag)
    write(run%io,'(a,i5,a,i5)') "         *     Traces: ",down%trc%n, '                 Shots: ',down%evnt%n
    write(run%io,*)                         "        *                [n]    [o]      [d] "
    write(run%io,'(a,i4,a,f7.3,f9.3,a,f5.2,a)') "         *           W: ",X%aw%n," ",X%aw%o,X%aw%d,&
    " =",real(down%aw%n)/real(down%nws)," blocks"
    if(allocated(tmp)) &
    write(run%io,'(a,f6.2,a)') '         *      Temporary read array: ',size(tmp)/1024./1024.," Mbytes"
    write(run%io,*) '        *      Headers: ', size(Y%bin)/1024. ," Kbytes"

    stat=sseek(Y%tag,0,0)
    Y%bin=(0.,0.)
    data_sort_init=T

#ifdef DEBUG
    if(run%debug) call out(unit)
#endif
  end function data_sort_init
!----------------------------------------------------------------
!  uses header file HD_d to locate sources
  logical function source_sort_init(Y,tag)
    type(dat)::Y
    character(len=128),intent(in):: tag
    integer::ix,iy,i,stat
    real::a,lx,ly,minx,maxx,miny,maxy
#ifdef DEBUG
    character(len=128), parameter :: unit='source_sort_init'
    if(run%debug) call in(unit)
#endif
    source_sort_init=F
    Y%tag=tag
    call auxinout(Y%tag)

    call hinaxis(tag,Y%amy,1); 
    call hinaxis(tag,Y%amx,2); 

    if(Y%amy%n/=2) call seperr("I want a 2D acquisition plane.")
    i=2*Y%amx%n*4

    allocate(Y%bin(Y%amx%n,1,1))
    Y%bin=cmplx(0.,0.)

    if(run%lock) call sep_open_lock(Y%tag)
    stat=sreed(Y%tag,Y%bin,i)
    if(run%lock) call sep_close_lock(Y%tag)
    if(stat/=i)  call seperr("trouble reading in source map (HD_d)")

    if(down%xyr) then
       do i=1,Y%amx%n
          lx=aimag(Y%bin(i,1,1))
          ly=real(Y%bin(i,1,1))
          Y%bin(i,1,1)=cmplx(lx,ly)
       end do
    end if

    if(Y%amx%n.ne.down%evnt%n.and.down%gathers) then
       write(run%io,*)'        * Image space may be Huge! Source headers '
       write(run%io,*)'        *  (HD_d) 2axis /= number shots in data (D).'
    end if

    write(run%io,*)'        *   Source locations provided in file: ',trim(Y%tag)
    write(run%io,'(a,i4)')'         *        number of shots:',Y%amx%n
    call smap_init_HD_d(Y)
    if(.not.down%gathers) down%evnt%n=Y%amx%n

    a=8*(size(Y%bin))
    down%mem=down%mem+a
    source_sort_init=T
#ifdef MEMCHECK
    if(run%debug) write(run%io,*) 'D.F90 HD_d%bin: allocate',a/1024./1024.,'Mbytes.  Total=',down%mem
#endif
#ifdef DEBUG
    if(run%debug) call out(unit)
#endif
  end function source_sort_init
!----------------------------------------------------------------
!
!  Initialize a tag from history
!   check to make sure that is dimensions dont disagree
!    with the general params
  logical function data_history_init(X)
    type(dat) :: X
    integer:: bogus
#ifdef DEBUG
    character(len=128), parameter :: unit='data_history_init'
    if(run%debug) call in(unit)
#endif

    call auxinout(X%tag)
    data_history_init=F

    X%from = "history"
    call from_aux(X%tag,'esize',X%esize,4)
    if(X%esize==4) call seperr("data has wrong esize, input frequency")

    !! . . Initialize Axes calls from history file
    call hinaxis(X%tag,X%amx,1); X%amx%l='mx'
    call hinaxis(X%tag,X%amy,2); X%amy%l='my'
    call hinaxis(X%tag,X%aw ,3); X%aw%l='w_'
    call hinaxis(X%tag,X%asx,4); X%asx%l='sx'
    call hinaxis(X%tag,X%asy,5); X%asy%l ='sy'

    if(.not.down%headD_d) then !! source positions from D file
       if(.not.down%offset) then!! turn it into offset
          X%amx%o=X%amx%o-X%asx%o
          X%amy%o=X%amy%o-X%asy%o
       end if
    end if

    !! workspace defined by input parameters if using headers
    if(.not.down%headD) then  
       down%amx=X%amx;down%amy=X%amy
       down%amx%o=X%amx%o+X%asx%o
       down%amy%o=X%amy%o+X%asy%o
    end if

    !! . . geometry check
    if((X%amx%d .lt. 0.).or.(X%amy%d .lt. 0.))&
    call seperr("Please refrain from negative d1 and/or d2")

    !! . . Default indices to read into prop. grid
    X%indx=1;                      X%indy=1
    X%finx=X%indx + X%amx%n -1;    X%finy=X%indy + X%amy%n -1

    if(.not. data_setup(X,.false.)) then
       write(run%io,*) "Trouble setting up data in data_history_init"
       return
    end if

    if(.not.down%headD_d)call smap_init(X)

    if(run%verb.gt.2) call data_report(X)

    data_history_init=T
#ifdef DEBUG
    if(run%debug) call data_report(X)
    if(run%debug) call out(unit)
#endif
  end function data_history_init
  !----------------------------------------------------------------
  !
  ! almost the same as data_history_init, leaving out assigns to 
  !     down structure
  logical function source_history_init(X)
    type(dat) :: X
    integer :: er
    real::a
#ifdef DEBUG
    character(len=128), parameter :: unit='source_history_init'
    if(run%debug) call in(unit)
#endif
    source_history_init=F
    call auxinout(X%tag)
    X%from = "history"
    call from_aux(X%tag,'esize',X%esize,4)
    if(X%esize==4) call seperr("source has wrong esize, input frequency")
!! . . Initialize Axes calls from history file

    X%amx = down%fatx!    call hinaxis(X%tag,X%amx,1)
    X%amy = down%faty!    call hinaxis(X%tag,X%amy,2)
    call hinaxis(X%tag,X%aw ,3)
    call hinaxis(X%tag,X%asx,4)
    call hinaxis(X%tag,X%asy,5)
    if(.not.down%gathers) then
       down%asx=X%asx; down%asy=X%asy
    end if

!! Frequency comparisson to data file
    if (X%aw%n /= down%aw%n) call erexit("Source file wrong length")   
    if (X%aw%d /= down%aw%d) call erexit("Source file freq. increment wrong")
    if (X%aw%o /= down%aw%o) call erexit("Source file freq. origin wrong")

!! . . Default indices to read into prop. grid
    X%indx=1;          X%indy=1
    X%finx=X%amx%n;    X%finy=X%amy%n

    if(.not. data_setup(X,.false.)) then
       write(run%io,*) "trouble setting up source files"
       return
    end if
    
    !!re-allocate array after distribution to size of
    !!propagation grid
    a=8*size(X%bin)
write(0,*)' source_history_init',a,size(X%bin)
    down%mem=down%mem-a
    deallocate(X%bin)

    allocate(X%bin(down%amx%n,down%amy%n,down%nws),stat=er)
    if(er/=0) call seperr('trouble at source_history_init reallocate D_d.bin')
    X%bin=(0.0,0.0)
write(0,*)' source_history_init',a,size(X%bin)
    !! little read in array
    if(down%fatx%n.eq.1.and.down%faty%n.eq.1) then
       allocate(source(down%nws))
    end if
    a=8*(size(X%bin)+size(source))
    down%mem=down%mem+a
#ifdef MEMCHECK
    if(run%debug) write(run%io,*) 'D.F90 D_d%bin: allocate',a,'Mbytes.  Total=',down%mem/1024./1024.
#endif

    if(run%verb.gt.2) call data_report(X)
    source_history_init=T
#ifdef DEBUG
!    if(run%debug) call data_report(X)
    if(run%debug) call out(unit)
#endif
  end function source_history_init
!----------------------------------------------------------------
  subroutine image_size_calc()
    real::n,m,g,minsx,minsy,maxsx,maxsy
#ifdef DEBUG
    character(len=128), parameter :: unit='image_size_calc'
    if(run%debug) call in(unit)
#endif

    if(.not.down%add) then !!if add, coordinates are set in init. of output file
    !! . . Define propagation wavefield on model surface
       minsx=down%asx%o
       minsy=down%asy%o
       minsx=minval(real(smap%s))
       minsy=minval(aimag(smap%s))
       maxsx=maxval(real(smap%s))
       maxsy=maxval(aimag(smap%s))
       
       if(down%roll) then !! rolling acqu. total image space bigger than each shot 
          if(down%amx%n > 1)then
             !! total image space
             down%aix%d=down%amx%d
!             down%aix%o=min(down%amx%o,minsx)
             down%aix%o=minsx+down%amx%o
!             n=max(maxsx,maxsx-(minsx-down%amx%o)+(down%amx%n-1)*down%amx%d)!!last possible trace
             n=max(maxsx,maxsx+down%amx%o+(down%amx%n-1)*down%amx%d)!!last possible trace
             down%aix%n=1.5+(n-down%aix%o)/down%aix%d
          else
             down%aix%n=1;down%aix%d=1.;down%aix%o=0.
          end if
          if(down%amy%n > 1)then
             down%aiy%d=down%amy%d
             down%aiy%o=minsy+down%amy%o
             n=max(maxsy,maxsy+down%amy%o+(down%amy%n-1)*down%amy%d)!!last possible trace
!             n=max(maxsy,maxsy-(minsy-down%amy%o)+(down%amy%n-1)*down%amy%d)!!last possible trace
             down%aiy%n=1.5+(n-down%aiy%o)/down%aiy%d
          else
             down%aiy%n=1;down%aiy%d=1.;down%aiy%o=0.
          end if
          !! define amx%o as a real earth location
          down%amx%o=real(smap%s(down%evnt%b))+down%amx%o
          down%amy%o=aimag(smap%s(down%evnt%b))+down%amy%o
       else !!  stationary acqu. total image space. workspace=imagespace
          if(down%amx%n > 1)then
             down%aix=down%amx
          else
             down%aix%n=1;down%aix%d=1.;down%aix%o=0.
          end if
          if(down%amy%n > 1)then
             down%aiy=down%amy
          else
             down%aiy%n=1;down%aiy%d=1.;down%aiy%o=0.
          end if
       end if
       n=down%aix%o
       m=down%aix%o+(down%aix%n-1)*down%aix%d
       down%aix%o=0.
       call from_param("minx",down%aix%o,n)
       call from_param("maxx",g,m)
       if(g.lt.down%aix%o) call seperr("maxx is less than minx")
       down%aix%n=1.5+(g-down%aix%o)/down%aix%d

       n=down%aiy%o
       m=down%aiy%o+(down%aiy%n-1)*down%aiy%d
       down%aiy%o=0.
       call from_param("miny",down%aiy%o,n)
       call from_param("maxy",g,m)
       if(g.lt.down%aiy%o) call seperr("maxy is less than miny")
       down%aiy%n=1.5+(g-down%aiy%o)/down%aiy%d
    end if

    !! sometimes you wont use the whole wavefield- so introduce ax,ay
    down%ax=down%amx;down%ay=down%amy
#ifdef DEBUG
    if(run%debug) call out(unit)
#endif
  end subroutine image_size_calc
!----------------------------------------------------------------
  subroutine source_pos(X,isx)
    type(dat):: X
    integer :: isx
    real::ly,lx
    if(down%passive.or.down%datum)then
       X%indx=1;X%finx=down%amx%n
       X%indy=1;X%finy=down%amy%n
    else
       lx=real(smap%s(isx))
       ly=aimag(smap%s(isx))
       if(down%amx%n.gt.1)then
          if(X%amx%n .lt. down%amx%n/2 +1 )then
             X%indx= (1.5+  (lx- down%amx%o)/down%amx%d) -(X%amx%n-1)/2
          end if
          X%finx=X%indx+X%amx%n-1
          if(X%indx.lt.1.or.X%finx.gt.down%amx%n) call seperr("Source location not in wavefield (x). ")
       end if

       if(down%amy%n.gt.1)then
          if(X%amy%n .lt. down%amy%n/2 +1 )then
             X%indy= (1.5+ (ly - down%amy%o)/down%amy%d) -(X%amy%n-1)/2
          end if
          X%finy=X%indy+X%amy%n-1
          if(X%indy.lt.1.or.X%finy.gt.down%amy%n) call seperr("Source location not in wavefield (y). ")
       end if
    end if
    call source_report(X)
  end subroutine source_pos
!----------------------------------------------------------------
!   relocate wavefield within the model
!    remember: amx%o is a relative distance from the shot!
  subroutine data_position(X,isx)
    type(dat) :: X
    integer , intent(in) :: isx
    if(down%roll) then          
       down%amx%o= X%amx%o + real(smap%s(isx))
       down%amy%o= X%amy%o + aimag(smap%s(isx))
    end if
    call data_position_report()
  end subroutine data_position
  !----------------------------------------------------------------
  subroutine data_position_report()
    if (run%verb.gt.3) then
       if(down%amx%n.gt.1)&
       write(run%io,'(a,f11.3)')"         *  Working grid positioned at x=",down%amx%o
       if (down%amy%n.gt.1)&
       write(run%io,'(a,f11.3)')"         *                             y=",down%amy%o
    end if
  end subroutine data_position_report
  !----------------------------------------------------------------
  !
  !
  !
  logical function create_data_file(X,tag)
    type(dat):: X
    character(len=*)         :: tag
       create_data_file=F 
       call auxscr(tag)
       call to_history("esize",8,tag)
       call houaxes(tag,X%amx,X%amy,X%aw,X%asx,X%asy)
       create_data_file=T 
     end function create_data_file
  !----------------------------------------------------------------
  logical function data_param_init(X)
    type(dat) :: X
#ifdef DEBUG
    character(len=128), parameter :: unit='data_param_init' 
    if(run%debug) call in(unit)
#endif
    data_param_init=F

    X%from = "parameter"
    X%amx = down%amx; X%amy = down%amy; 
    X%asx = down%asx; X%asy = down%asy; 
    X%aw  = down%aw;  
    if(.not. create_data_file(X,X%tag)) then
      write(run%io,*) "trouble creating data file"
      return
    end if
    if(.not. data_setup(X,.true.)) then
      write(run%io,*) "trouble setting up data"
      return
    end if

    data_param_init=T
#ifdef DEBUG
    if(run%debug) call data_report(X)
    if(run%debug) call out(unit)
#endif
  end function data_param_init

!----------------------------------------------------------------
  subroutine data_assign(X,val)
   type(dat) :: X
   complex :: val
   logical ::a

    a=file_val_assign(X%tag,val)
  end subroutine data_assign
  !----------------------------------------------------------------
  !
  ! data_scratch_init
  !
  !  Create a data scratch file
  !
  !
  logical function data_scratch_init(X)
    type(dat) :: X
#ifdef DEBUG
    character(len=128), parameter :: unit='data_scratch_init' 
    if(run%debug) call in(unit)
#endif

    data_scratch_init=F
    if(down%passive) then
       X%from = "Data"
    else
       X%from= "parameter"
    end if
    X%amx = down%amx; X%amy = down%amy; 
    X%asx = down%asx; X%asy = down%asy
    X%aw  = down%aw; 
    if (down%passive) then
       X%asx%d=0.  ; X%asy%d = 0.
       X%asx%n=down%asx%n; X%asy%n=1
    endif

    X%indx=1;                      X%indy=1
    X%finx=X%amx%n;    X%finy=X%amy%n

    if(.not. create_data_file(X,X%tag)) then
      write(run%io,*) "trouble  creating data file from scratch: ",trim(X%tag)
      return
    end if

    if(.not. data_setup(X,.true.)) then
      write(run%io,*) "trouble setting up data from scratch: ",trim(X%tag)
      return
    end if

    data_scratch_init=T
    if(run%verb.gt.2)call data_report(X)
#ifdef DEBUG
    if(run%debug) call out(unit)
#endif
  end function data_scratch_init

  !----------------------------------------------------------------
  !  The core  routine to setup data
  logical function data_setup(X,zero)
    type(dat):: X
    logical :: zero
    integer :: er
    real::a
#ifdef DEBUG
    character(len=128), parameter :: unit='data_setup' 
    if(run%debug) call in(unit)
#endif
    data_setup=F
    X%bites = 8* product( (/X%amx%n,X%amy%n,down%nws/) )

    !allocate array
    allocate(X%bin(X%amx%n,X%amy%n,down%nws),stat=er)
    if(er/=0) call seperr('cant allocate D.bin')
    X%bin=(0.0,0.0)

    a=8*size(X%bin)
    down%mem=down%mem+a
#ifdef MEMCHECK
    if(run%debug) write(run%io,*) 'D.F90 binary allocate ',trim(X%tag),' ',a/1024./1024.,'Mbytes.  Total=',down%mem/1024./1024.
#endif

    data_setup=T
#ifdef DEBUG
    if(run%debug) call out(unit)
#endif
  end function data_setup
!----------------------------------------------------------------
! Information on the data
!
  subroutine data_report(X)
    type(dat) :: X
    real::bb
    bb=real(down%aw%n)/real(down%nws)
    write(run%io,*)'        ************************************************'
    write(run%io,*)'        *   ',trim(X%tag),' from ',trim(X%from)
    write(run%io,*)                         "        *                [n]     [o]      [d] "       
    write(run%io,'(a,i5,f9.2,f9.2)       ') "         *        Offx: ",X%amx%n, X%amx%o, X%amx%d
    write(run%io,'(a,i5,f9.2,f9.2)       ') "         *        Offy: ",X%amy%n, X%amy%o, X%amy%d
    write(run%io,'(a,i5,f9.2,f9.2)       ') "         *        Shtx: ",X%asx%n, X%asx%o, X%asx%d
    write(run%io,'(a,i5,f9.2,f9.2)       ') "         *        Shty: ",X%asy%n, X%asy%o, X%asy%d
    write(run%io,'(a,i5,f9.5,1x,f9.7,a,f5.2,a)') "         *   Frequency: ",X%aw%n,X%aw%o,X%aw%d,&
    " =",bb," blocks"
    write(run%io,'(a,f6.2,a)        ') "         *       block: ",8.*real(size(X%bin))/1024./1024.," Mbytes"
  end subroutine data_report
!--------------------------------------------------------------
! Read the data 
  logical function data_read(X,k)
    type(dat) :: X
    integer, intent(in) :: k
    integer :: stat,i,j,l
#ifdef DEBUG
    character(len=128), parameter :: unit='data_read'
    if(run%debug) call in(unit)
#endif

    data_read=F
    X%bin=(0.0,0.0)
    X%bites=down%amx%n*down%amy%n*down%nws*8
    if(run%lock) call sep_open_lock(X%tag)
    stat= sreed(X%tag, X%bin, X%bites)
    if(run%lock) call sep_close_lock(X%tag)
    if(stat.ne.X%bites)return

    if(down%scale.ne.1.0) then
!$OMP PARALLEL DO PRIVATE(i,j,l)
       do i=1,down%nws
       do j=1,down%amy%n
       do l=1,down%amx%n
          X%bin(l,j,i)=down%scale*X%bin(l,j,i)
       end do
       end do
       end do
!$OMP END PARALLEL DO 
    end if
    data_read=T
#ifdef DEBUG
    if(run%debug) call out(unit)
#endif
  end function data_read
!----------------------------------------------------------------
  subroutine data_header_read(Y)
    type(dat)::Y
    integer:: ix,iy,stat,itr

       Y%bin=cmplx(0.,0.)
       stat=down%trc%n*8
       if(run%lock) call sep_open_lock(Y%tag)
       if(stat/= sreed(Y%tag,Y%bin,down%trc%n*8)) &
       call seperr("Poor reading of headers in data_header_read.")
       if(run%lock) call sep_close_lock(Y%tag)
       if(down%offset) Y%bin=Y%bin+smap%s(1)  !! add source pos. if input in offset coords

       if(.not.down%xyr) then
          do itr=1,down%trc%n
             ix=1.5+(real (Y%bin(itr,1,1))-down%amx%o)/down%amx%d
             iy=1.5+(aimag(Y%bin(itr,1,1))-down%amy%o)/down%amy%d
             Y%bin(itr,1,1)=cmplx(real(ix),real(iy))
          end do
       else
          do itr=1,down%trc%n
             iy=1.5+(real (Y%bin(itr,1,1))-down%amx%o)/down%amx%d
             ix=1.5+(aimag(Y%bin(itr,1,1))-down%amy%o)/down%amy%d
             Y%bin(itr,1,1)=cmplx(real(ix),real(iy))
          end do
       end if
    if(Y%asx%n.eq.1)then!!receivers changing position
       stat=sseek(Y%tag,0,0)
    end if
  end subroutine data_header_read  
!----------------------------------------------------------------
! read and sort traces
  logical function data_sort_read(X,X_d,k,Y)
    type(dat) :: X,X_d !!D,D_d
    type(dat) :: Y !!receiver headers
    integer, intent(in) :: k
    integer :: stat,ix,iy,iw,itr,i1,i2,i3,n_dim_cube,badtraces
    integer,dimension(3)::n_cube,n_wind,f_wind,j_wind
    real::a
#ifdef DEBUG
    character(len=128), parameter :: unit='data_sort_read'
    if(run%debug) call in(unit)
#endif

    data_sort_read=F
    X%bin=(0.0,0.0)
    bincount=0
    badtraces=0
    n_dim_cube=3
    n_cube=(/down%aw%n    ,down%trc%n,X%asx%n     /)
    j_wind=(/1            ,1         ,1           /)
    n_wind=(/down%nws     ,down%trc%n,1           /)
    f_wind=(/k*down%wblock,0         ,down%shotnum/)
    if(run%lock) call sep_open_lock(X%tag)
    stat=sreed_window(X%tag,n_dim_cube,n_cube,n_wind,f_wind,j_wind,8,tmp(1,1))
    if (stat.lt.0) then
       write(run%io,*) "Badness read in data_sort_read for file ",trim(X%tag)
    end if
    if(run%lock) call sep_close_lock(X%tag)
    if(.not.down%passive) then
       do itr=1,down%trc%n
          ix=int(real (Y%bin(itr,1,1)))
          if(ix.gt.down%amx%n.or.ix.lt.1)then
             badtraces=badtraces+1
             write(run%io,'(10x,a,i5,a,i5,a,i4)') & 
             'Trace ',itr,' out of range.  X Index:',ix,' Bad traces=',badtraces
          else
             iy=int(aimag(Y%bin(itr,1,1)))
             if(iy.gt.down%amy%n.or.iy.lt.1)then
                badtraces=badtraces+1
                write(run%io,'(10x,a,i5,a,i5,a,i4)') &
                'Trace ',itr,' out of range. Y Index:',iy,' Bad traces=',badtraces
             else
                bincount(ix,iy)=bincount(ix,iy)+1
                do iw=1,down%nws
                   X%bin(ix,iy,iw)=tmp(iw,itr)
                end do
             end if
          end if
       end do
    else
       do itr=1,down%trc%n
          ix=int(real (Y%bin(itr,1,1)))
          if(ix.gt.down%amx%n.or.ix.lt.1)then
             badtraces=badtraces+1
             if(run%verb.gt.0) write(run%io,'(10x,a,i5,a,i5,a,i4)') &
             'Trace ',itr,' out of range. Index:',ix,' Bad traces=',badtraces
          else
             iy=int(aimag(Y%bin(itr,1,1)))
             if(iy.gt.down%amy%n.or.iy.lt.1)then
                badtraces=badtraces+1
                if(run%verb.gt.0) write(run%io,'(10x,a,i5,a,i5,a,i4)') &
                'Trace ',itr,' out of range. Index:',iy,' Bad traces=',badtraces
             else
                bincount(ix,iy)=bincount(ix,iy)+1
                do iw=1,down%nws
                   X%bin(ix,iy,iw)=      tmp(iw,itr)
                   X_d%bin(ix,iy,iw)=conjg(tmp(iw,itr))
                end do
             end if
          end if
       end do

       do i3=1,down%nws
          do i2=1,down%amy%n
             do i1=1,down%amx%n
                if(bincount(i1,i2).gt.0) then
                   X_d%bin(i1,i2,i3)=X_d%bin(i1,i2,i3)/bincount(i1,i2)
                end if
             end do
          end do
       end do
    end if

    do i3=1,down%nws
       do i2=1,down%amy%n
          do i1=1,down%amx%n
             if(bincount(i1,i2).gt.0) then
                X%bin(i1,i2,i3)=X%bin(i1,i2,i3)/bincount(i1,i2)
             end if
          end do
       end do
    end do
    if(badtraces.gt.0) then
       a=100*real(badtraces)/real(down%trc%n)
       write(run%io,'(5x,a,i5,a,i5,a,i2,a)')&
       'Are you OK with chucking',badtraces,' out of',down%trc%n,' (',nint(a),'% of your data)'
    end if

#ifdef DEBUG
    if(run%debug) call out(unit)
#endif
    data_sort_read=T
  end function data_sort_read
!--------------------------------------------------------------
! Read the source
  logical function source_read(X,k)
    type(dat) :: X
    integer, intent(in) :: k
    integer :: stat,i,j,l
    real:: tw

#ifdef DEBUG
    character(len=128), parameter :: unit='source_read'
    if(run%debug) call in(unit)
#endif

    source_read=F
!$OMP PARALLEL DO PRIVATE(i,j,k)
       do i=1,down%nws
          do j=1,down%amy%n
             do l=1,down%amx%n
                X%bin(l,j,i)=cmplx(0.0,0.0)
             end do
          end do
       end do
!$OMP END PARALLEL DO 

    if(down%fatx%n.eq.1.and.down%faty%n.eq.1) then
       if(run%lock) call sep_open_lock(X%tag)
       stat= sreed(X%tag,source,8*down%nws)
       if(run%lock) call sep_close_lock(X%tag)
       if(stat.ne.8*down%nws) return

       if(down%modelplane) then
          write(run%io,*)'        * Modeling source function as a refl. from a plane.'
          call planemod(source,X%bin)
       else
          if(down%tshift.ne.0.) then
             do i=1,down%nws
                tw=down%tshift*(down%bw%o+(i-1)*down%aw%d)
                X%bin(X%indx,X%indy,i)=source(i)*cmplx(cos(tw),sin(tw))
             end do
          else
             do i=1,down%nws
                X%bin(X%indx,X%indy,i)=conjg(source(i))
             end do
          end if
       end if
    else
       if(run%lock) call sep_open_lock(X%tag)
       stat=sreed(X%tag, X%bin, 8*size(X%bin))
       if(run%lock) call sep_close_lock(X%tag)
       if(stat.ne.8*size(X%bin)) return
!$OMP PARALLEL DO PRIVATE(i,j,k)
       do i=1,down%nws
          do j=1,down%amy%n
             do l=1,down%amx%n
                X%bin(l,j,i)=conjg(X%bin(l,j,i))
             end do
          end do
       end do
!$OMP END PARALLEL DO 
    end if
    source_read=T
#ifdef DEBUG
    write(run%io,*) 'source_read',source_read
    write(run%io,*)' X%indx:X%finx,X%indy:X%finy',X%indx,X%finx,X%indy,X%finy
    if(run%debug) call out(unit)
#endif
  end function source_read
!----------------------------------------------------------------
subroutine source_report(X)
  type(dat)::X
  if(down%passive.or.down%datum) return

  if( run%verb.gt.2) then
     if (X%amx%n==1) then !! single trace source
        write(run%io,'(a,f11.3)') '         *           Source located at x=',&
        down%amx%o + (X%indx-1)*down%amx%d +((X%amx%n-1)/2)*down%amx%d
     else  !! areal source
        write(run%io,'(a,f11.3)') '         *           Source located at x=',&
        down%amx%o + X%amx%o
     end if
     
     if  (X%amy%n==1) then
        write(run%io,'(a,f11.3)')    '         *                             y=',&
        down%amy%o + (X%indy-1)*down%amy%d +((X%amy%n-1)/2)*down%amy%d
     else
        write(run%io,'(a,f11.3)') '         *                             y=',&
        down%amy%o + X%amy%o
     end if
  end if
end subroutine source_report

!----------------------------------------------------------------
!  Copy data
  subroutine data_copy_conjg(X,Y)
    type(dat) :: X
    type(dat) :: Y
    integer :: iwb,i,j,k
!    logical :: res
#ifdef DEBUG
    character(len=128), parameter :: unit='data_copy_conjg'
    if(run%debug) call in(unit)
#endif
    if(size(Y%bin).ne.size(X%bin))call seperr("Wavefields D and D_d are different size.")
    if (.not.down%fwdscat) then
!$OMP PARALLEL DO PRIVATE(i,j,k) 
       do i=1,down%nws
       do j=1,down%amy%n
       do k=1,down%amx%n
          Y%bin(k,j,i) = conjg(X%bin(k,j,i))
       end do; end do; end do
!$OMP END PARALLEL DO 
    else
!$OMP PARALLEL DO PRIVATE(i,j,k) 
       do i=1,down%nws
       do j=1,down%amy%n
       do k=1,down%amx%n
          Y%bin(k,j,i) = X%bin(k,j,i)
       end do; end do; end do
!$OMP END PARALLEL DO 
    end if

#ifdef DEBUG
    if(run%debug) call out(unit)
#endif
  end subroutine data_copy_conjg
!----------------------------------------------------------------
!  Copy data
  subroutine data_conjg(X)
    type(dat) :: X
    integer :: i,j,k
#ifdef DEBUG
    character(len=128), parameter :: unit='data_conjg'
    if(run%debug) call in(unit)
#endif
    if (down%fwdscat) then
       return
    else
!$OMP PARALLEL DO PRIVATE(i,j,k) 
       do i=1,down%nws
       do j=1,down%amy%n
       do k=1,down%amx%n
          X%bin(k,j,i) = conjg(X%bin(k,j,i))
       end do; end do; end do
!$OMP END PARALLEL DO 
    end if
#ifdef DEBUG
    if(run%debug) call out(unit)
#endif
  end subroutine data_conjg

!----------------------------------------------------------------
  subroutine smap_init(X)
    type(dat) :: X
    integer :: er,i,iy,ix,nx,ny
    real::a,ly,lx
#ifdef DEBUG
    character(len=128), parameter :: unit='smap_init'
    if(run%debug) call in(unit)
#endif

    down%evnt%b=1
    nx=max(X%asx%n,down%nx_shots,down%asx%n)
    ny=max(X%asy%n,down%ny_shots)
    down%asx%o=min(X%asx%o,down%ox_shots)
    down%asy%o=min(X%asy%o,down%oy_shots)
    down%asx%d=X%asx%d;    down%asy%d=X%asy%d
    
    allocate(smap%s(nx*ny))
    i=1
    do iy=1,ny
       ly=down%asy%o + (iy-1)* down%asy%d
       do ix=1,nx
          lx=down%asx%o + (ix-1)* down%asx%d
          smap%s(i)=cmplx(lx,ly)
          if((abs(ly-X%asy%o).lt.X%asy%d/2.).and.(abs(lx-X%asx%o).lt.X%asx%d/2.)) down%evnt%b=i
          i=i+1
       end do
    end do
    
    a=8*size(smap%s)
    down%mem=down%mem+a
#ifdef MEMCHECK
    if(run%debug) write(run%io,*) 'D.F90 smap allocate ',a/1024./1024.,'Mbytes.  Total=',down%mem/1024./1024.
#endif

#ifdef DEBUG
    if(run%debug) call out(unit)
#endif
  end subroutine smap_init
!----------------------------------------------------------------
subroutine smap_init_HD_d(X)
    integer :: er
    type(dat) :: X
#ifdef DEBUG
    character(len=128), parameter :: unit='smap_init_HD_d'
    if(run%debug) call in(unit)
#endif
    allocate(smap%s(size(X%bin,1)))
    smap%s= X%bin(:,1,1)

    down%asx%o=minval(real(smap%s))
    down%asy%o=minval(aimag(smap%s))

    down%asx%n=down%evnt%n
    down%asy%n=1
#ifdef DEBUG
    if(run%debug) call out(unit)
#endif
  end subroutine smap_init_HD_d
!---------------------------------------------------------------- 
  subroutine planemod(func,wfld)
    complex,dimension(:,:,:):: wfld
    complex,dimension(:)::func
    real::sx,sy,x,y,lx,ly,ll,lt,z,z2,t0,s,s2,x2,a,az,d,ww,la
    integer:: ix,iy,iw
    x=real(smap%s(down%shotnum+1))
    y=aimag(smap%s(down%shotnum+1))

    call from_param("plane_z",z,0.)
    call from_param("plane_v",s,1465.)
    s=1./s
    call from_param("plane_az",az,0.)
    call from_param("plane_d",d,0.)
    call from_param("plane_x",sx,18050.)
    call from_param("plane_y",sy,-34350.)

    !! remove azimuth from trial source loc to receiver 
    la=atan2(sy-y,sx-x)

    a=la-az
    if(a.lt.0.)a=pi+a
    a=cos(a)
    d=sin(d)

    !! distances
    lx=(x-sx)**2
    ly=(y-sy)**2
    x2=lx+ly    
    lx=sqrt(x2)
    z2=z*z

    !! vel squared
    s2=s*s
    !! term 1
    lt=4*s2*z2+s2*x2
    !! term 2
    ll=4*z*lx*d*a*s2

    !! event at source trace (correlated gathers) always at t=0
    t0=sqrt(lt+ll)

!!$OMP PARALLEL DO PRIVATE(w,y,ly,x,lx,x2,a,la,lt,ll)
    do iw=1,down%nws
       ww= down%bw%o + (iw-1)*down%aw%d
       do iy=1,down%amy%n
          y=down%amy%o+(iy-1)*down%amy%d
          ly=(y-sy)**2
          do ix=1,down%amx%n
             x=down%amx%o+(ix-1)*down%amx%d
             la=atan2(sy-y,sx-x)
             a=la-az
             if(a.lt.0.)a=pi+a
             a=cos(a)
             lx=(x-sx)**2
             x2=lx+ly    
             lx=sqrt(x2)
             lt=4*s2*z2+s2*x2
             ll=4*z*lx*d*a*s2
             a=ww*(sqrt(lt+ll)-t0)
             wfld(ix,iy,iw)=(func(iw)*cmplx(cos(a),-sin(a)))
          end do
       end do
    end do
!!$OMP END PARALLEL DO 
  end subroutine planemod
!----------------------------------------------------------------
!  Find the next good fft length
!
integer function next_fast(nlow)
integer ::i,nlow
integer ,allocatable::igood_len(:)
logical :: found

allocate(igood_len(800))

igood_len=(/ 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,18,20,21,22,24,25,26,27,28,30,32,33,35,36,39,40,42,44,45,48,&
	& 49,50,52,54,55,56,60,63,64,65,66,70,72,75,77,78,80,81,84,88,90,91,96,98,99,100,104,105,108,110,112,117,120,&
	& 125,126,128,130,132,135,140,144,147,150,154,156,160,162,165,168,175,176,180,182,189,192,195,196,198,200, &
	& 208,210,216,220,224,225,231,234,240,243,245,250,252,256,260,264,270,273,275,280,288,294,297,300,308,312, &
	& 315,320,324,325,330,336,343,350,351,352,360,364,375,378,384,385,390,392,396,400,405,416,420,432,440,441, &
	& 448,450,455,462,468,480,486,490,495,500,504,512,520,525,528,539,540,546,550,560,567,576,585,588,594,600, &
	& 616,624,625,630,637,640,648,650,660,672,675,686,693,700,702,704,720,728,735,750,756,768,770,780,784,792, &
	& 800,810,819,825,832,840,864,875,880,882,891,896,900,910,924,936,945,960,972,975,980,990,1000,1008,1024,  &
	& 1029,1040,1050,1053,1056,1078,1080,1092,1100,1120,1125,1134,1152,1155,1170,1176,1188,1200,1215,1225,1232, &
	& 1248,1250,1260,1274,1280,1296,1300,1320,1323,1344,1350,1365,1372,1375,1386,1400,1404,1408,1440,1456,1470, &
	& 1485,1500,1512,1536,1540,1560,1568,1575,1584,1600,1617,1620,1625,1638,1650,1664,1680,1701,1715,1728,1750, &
	& 1755,1760,1764,1782,1792,1800,1820,1848,1872,1875,1890,1911,1920,1925,1944,1950,1960,1980,2000,2016,2025, &
	& 2048,2058,2079,2080,2100,2106,2112,2156,2160,2184,2200,2205,2240,2250,2268,2275,2304,2310,2340,2352,2376, &
	& 2400,2430,2450,2457,2464,2475,2496,2500,2520,2548,2560,2592,2600,2625,2640,2646,2673,2688,2695,2700,2730, &
	& 2744,2750,2772,2800,2808,2816,2835,2880,2912,2925,2940,2970,3000,3024,3072,3080,3087,3120,3125,3136,3150, &
	& 3159,3168,3185,3200,3234,3240,3250,3276,3300,3328,3360,3375,3402,3430,3456,3465,3500,3510,3520,3528,3564, &
	& 3584,3600,3640,3675,3696,3744,3750,3773,3780,3822,3840,3850,3888,3900,3920,3960,3969,4000,4032,4050,4095, &
	& 4096,4116,4125,4158,4160,4200,4212,4224,4312,4320,4368,4375,4400,4410,4455,4459,4480,4500,4536,4550,4608, &
	& 4620,4680,4704,4725,4752,4800,4851,4860,4875,4900,4914,4928,4950,4992,5000,5040,5096,5120,5145,5184,5200, &
	& 5250,5265,5280,5292,5346,5376,5390,5400,5460,5488,5500,5544,5600,5616,5625,5632,5670,5733,5760,5775,5824, &
	& 5850,5880,5940,6000,6048,6075,6125,6144,6160,6174,6237,6240,6250,6272,6300,6318,6336,6370,6400,6468,6480, &
	& 6500,6552,6600,6615,6656,6720,6750,6804,6825,6860,6875,6912,6930,7000,7020,7040,7056,7128,7168,7200,7280, &
	& 7350,7371,7392,7425,7488,7500,7546,7560,7644,7680,7700,7776,7800,7840,7875,7920,7938,8000,8064,8085,8100, &
	& 8125,8190,8192,8232,8250,8316,8320,8400,8424,8448,8505,8575,8624,8640,8736,8750,8775,8800,8820,8910,8918, &
	& 8960,9000,9072,9100,9216,9240,9261,9360,9375,9408,9450,9504,9555,9600,9625,9702,9720,9750,9800,9828,9856, &
	& 9900,9984,10000,10080,10125,10192,10240,10290,10368,10395,10400,10500,10530,10560,10584,10692,10752,10780,&
	&10800,10920,10976,11000,11025,11088,11200,11232,11250,11264,11319,11340,11375,11466,11520,11550,11648,11700,&
	&11760,11880,11907,12000,12096,12150,12250,12285,12288,12320,12348,12375,12474,12480,12500,12544,12600,12636,&
	&12672,12740,12800,12936,12960,13000,13104,13125,13200,13230,13312,13365,13377,13440,13475,13500,13608,13650,&
	&13720,13750,13824,13860,14000,14040,14080,14112,14175,14256,14336,14400,14553,14560,14625,14700,14742,14784,&
	&14850,14976,15000,15092,15120,15288,15360,15400,15435,15552,15600,15680,15750,15795,15840,15876,15925,16000,&
	&16128,16170,16200,16250,16380,16384,16464,16500,16632,16640,16800,16848,16875,16896,17010,17150,17199,17248,&
	&17280,17325,17472,17500,17550,17600,17640,17820,17836,17920,18000,18144,18200,18375,18432,18480,18522,18711,&
	&18720,18750,18816,18865,18900,19008,19110,19200,19250,19404,19440,19500,19600,19656,19712,19800,19845,19968,&
	&20000,20160,20250,20384,20475,20480,20580,20625,20736,20790,20800,21000,21060,21120,21168,21384,21504,21560,&
	&21600,21840,21875,21952,22000,22050,22113,22176,22275,22295,22400,22464,22500,22528,22638,22680,22750,22932,&
	&23040,23100,23296,23400,23520,23625,23760,23814,24000,24192,24255,24300,24375,24500,24570,24576,24640,24696,&
	&24750,24948,24960 /)

  i=1; found=.false.
  do while(i < size(igood_len) .and. .not. found)
    if(igood_len(i)< nlow) then
     i=i+1
    else
     found=.true.
    end if
  end do
  next_fast=igood_len(i)
  deallocate(igood_len)
end function next_fast
  !----------------------------------------------------------------  
end module D_mod
