Program BigCorr
  !##################################################
  !#                                              
  !#  cross-correlate every trace with every other
  !#     <in.H is complex, and n3 is freq, n1 is traces
  !#     >out.H bigger
  !##################################################
  use sep
  use flags
  use timers
  use cos_taper_mod 
  use smooth_mod
  implicit none
  integer :: j1,i1,j2,i2,i3,j3,a,b,c,d,i,j,k,k1,k2,k3,ompi,ompn,nof,pp
  integer :: n1,n2,n3,n4,n5,esize,n1out,blck,nb,er,nf,nt,trace,rect2,npad,nn
  real    :: o1,o2,o3,d1,d2,d3,df,maxf,minf,mem,scale,min1,max1,dt,small
  real   , pointer,dimension(:,:)  :: data,denom,tdat
  complex, pointer,dimension(:,:)  :: cdata
  complex, pointer,dimension(:,:,:):: res
  complex, pointer,dimension(:,:)  :: shottrace
  logical:: add,signum,scaled,tapered,decon
  integer*8,pointer:: planF(:)
  integer, external :: omp_get_thread_num,omp_get_num_threads,sreed_window,srite_window
  real:: startR,finishR,startW,finishW,startFFT,finishFFT, startMult, finishMult,start,finish
  integer,dimension(2)::n_cube,n_wind,f_wind,j_wind

  startR=0.;finishR=0.;startW=0.;finishW=0.;start=0.;finish=0.
  startFFT=0.;finishFFT=0.; startMult=0.; finishMult=0.;
  n3=1
  call start_timer(start)
  !----------------------------------------------------------------
  !!     header/parameter section
  call sep_init(SOURCE)
  call getch_add_string("head=/dev/null")
  call getch_add_string("noheader=y")
  call set_no_putch()
  call auxin("traces")
  call from_param('add',add,.false.); call from_param("signum",signum,.false.)
  call from_aux("traces",'n1',n1);call from_aux("traces",'n2',n2);call from_aux("traces",'n3',n3)
  call from_aux("traces",'o1',o1);call from_aux("traces",'o2',o2);!call from_aux("traces",'o3',o3)
  call from_aux("traces",'d1',d1);call from_aux("traces",'d2',d2);!call from_aux("traces",'d3',d3)
  call from_aux("traces",'esize',esize); if(esize==8) call seperr("need time input")
  call from_param("pad",npad,n1)
  if(npad.lt.n1) then
     write(0,*) ' npad is less than n1. I added',npad,'samples to n1=',n1
     npad=n1+npad
  end if
  if(n1.lt.npad) then
     write(0,*) 'padding 1 axis to ',npad
     allocate(tdat(n1,n2))
     nn=n1
     n1=npad
  else
     nn=n1
  end if
  esize=8
  call from_param("decon",decon,.false.);call from_param("rect",rect2,3)
  df=1./(d1*n1)
  call from_param('maxf',maxf,.5/d1); call from_param('minf',minf,0.)
  scale=0.025*(maxf-minf)
  call from_param('max1',max1,maxf-scale); call from_param('min1',min1,minf+scale)
  call from_param('blocks',blck,1); nb=n2/blck
  call from_param('small',small,0.0001)
  ompn=1
#ifdef SEP_OMP  
!$OMP PARALLEL 
  ompn = omp_get_num_threads()  
!$OMP END PARALLEL 
#endif
  write(0,*) 'OMP threads',ompn

  call init_cos_taper(n1/2+1,0.,df,minf,min1,maxf,max1,ompn)
  nof=i0_n+1
  n1out=i0_x-1
  nf=n1/2+1
  o1=minf

!write(0,*) 'nof',cos_taper_filter(nof-1,1),cos_taper_filter(nof,1),cos_taper_filter(nof+1,1)
!write(0,*) 'n1o',cos_taper_filter(n1out-1,1),cos_taper_filter(n1out,1),cos_taper_filter(n1out+1,1)
!write(0,*) ',nf,nof,n1out',nf,nof,n1out
!write(0,*) 'now i0_n,i1_n,i0_x,i1_x',i0_n,i1_n,i0_x,i1_x
!write(0,*) 'n1out',n1out,n1out-nof
  call auxinout("gathers")
  if(add) then
     call from_aux("gathers","esize",er)
     if(er.ne.8) call seperr("can not add frequency and time")
  else 
     call to_history('decon',decon,"gathers");call to_history('smoothrect',rect2,"gathers")
     call to_history('esize',esize,"gathers");call to_history('d1',df,"gathers")
     call to_history('n1',n1out-nof+1,"gathers");call to_history('o1',o1,"gathers")
     call to_history('n2',n2,"gathers");call to_history('o2',0.,"gathers");call to_history('d2',1.,"gathers")
     call to_history('n3',n2,"gathers");call to_history('o3',0.,"gathers");call to_history('d3',1.,"gathers")
     call to_history('n4',1,"gathers");call to_history('o4',0.,"gathers");call to_history('d4',1.,"gathers")
     call to_history('n5',1,"gathers");call to_history('o5',0.,"gathers");call to_history('d5',1.,"gathers")
     call to_history('maxf',minf,"gathers");  call to_history('maxf',maxf,"gathers")
     call to_history('maxf',min1,"gathers");  call to_history('maxf',max1,"gathers")
     call to_history('nt',n1,"gathers")
     call sep_close()
  end if
  write(0,*) '****************  blocks:',nb,blck
  write(0,*) 'Sign bit processing:',signum,'   Add :',add,'   Deconvolve :',decon

  !----------------------------------------------------------------
  !!     allocations
  allocate(data(n1,n2))
  call sep_begin_prog()
!  mem=mem+4*size(data)
!  write(0,*) 'allocate data',mem/1024/1024
  allocate(cdata(nf,n2))
!  mem=mem+8*size(cdata)
!  write(0,*) 'allocate cdata',mem/1024/1024
  allocate (planF(n2));planF=0
!  mem=mem+8*size(planF)
!  write(0,*) 'allocate planF',mem/1024/1024

  allocate(res(n1out-nof+1,n2,blck))
!  mem=mem+8*size(res)
!  write(0,*) 'allocate res',mem/1024/1024
  allocate(shottrace(nf,ompn))
!  mem=mem+8*size(shottrace)
!  write(0,*) 'allocate shottrace',mem/1024/1024
  shottrace=(0.,0.)
  allocate(denom(n1out,n2))
  denom=0.
  write(0,*) 'size(res,3),size(res,2),size(res,1)',size(res,3),size(res,2),size(res,1)
  !----------------------------------------------------------------
  !!    FFT Section
  call start_timer(startFFT)
  do i=1,n2
     call sfftw_plan_dft_r2c_1d(planF(i),n1,&
     data(1,i),cdata(1,i),FFTW_PATIENT)        
  end do
  do k2=1,size(cdata,2)
     do k1=1,size(cdata,1)
        cdata(k1,k2)=cmplx(0.,0.)
     end do
  end do
  do k2=1,size(data,2)
     do k1=1,size(data,1)
        data(k1,k2)=0.
     end do
  end do
  call end_timer(finishFFT,startFFT)

  !----------------------------------------------------------------  
  !! data shenanigans.
  !!   a) do the fft to freq.
  !!   b) scale the fft (2 way scaling)
  !!   c) make smoothed (x-axis) denom. for decon
  !!   d) apply cos_taper, scaling, and denom to cdata
do pp=1,n3
write(0,*)' Time panel',pp,'of',n3,'. Add :',add
  if(nn.ne.n1) then
     call sep_read(tdat,"traces")
     do i=1,n2
        do j=1,nn
           data(j,i)=tdat(j,i)
        end do
        do k=nn+1,n1
           data(j,i)=0.
        end do
     end do
  else
     call sep_read(data,"traces")
  end if
write(0,*)' sum(data)',sum(data)
  if (signum) then
     do k2=1,size(data,2)
        do k1=1,size(data,1)
           data(k1,k2)=data(k1,k2)/abs(data(k1,k2))
        end do
     end do
  end if

  scale=1./sqrt(real(n1))
  call start_timer(startFFT)
!$OMP PARALLEL DO PRIVATE(i,j)
  do i=1,n2
     call sfftw_execute(planF(i))
  end do
!$OMP END PARALLEL DO 
  call end_timer(finishFFT,startFFT)

  if(decon) then
!$OMP PARALLEL PRIVATE(ompi)
#ifdef SEP_OMP
     ompi = omp_get_thread_num()+1    
#endif
!$OMP DO PRIVATE( i2,i1)
     do i2=1,n2
        do i1=nof,n1out
           denom(i1,i2)=cabs(cdata(i1,i2))*scale
           if(denom(i1,i2).lt.small) denom(i1,i2)=1
           denom(i1,i2)=cos_taper_filter(i1,ompi)/denom(i1,i2)
        end do
     end do
!$OMP END DO 
!$OMP END PARALLEL 
     if(rect2.gt.1)call smooth(denom,1,rect2,0,1,0,0,0)
  else

!$OMP PARALLEL PRIVATE(ompi)
#ifdef SEP_OMP
     ompi = omp_get_thread_num()+1    
#endif
!$OMP DO PRIVATE( i2,i1)
     do i2=1,n2
        do i1=nof,n1out
           denom(i1,i2)=scale*cos_taper_filter(i1,ompi)
        end do
     end do
!$OMP END DO 
!$OMP END PARALLEL  
  end if
!$OMP PARALLEL DO PRIVATE( i2,i1)
  do i2=1,n2
     do i1=1,nof
        cdata(i1,i2)=cmplx(0.,0.)
     end do
     do i1=nof+1,n1out
        cdata(i1,i2)=cdata(i1,i2)*denom(i1,i2)
     end do
     do i1=n1out+1,size(cdata,1)
        cdata(i1,i2)=cmplx(0.,0.)
     end do
  end do
!$OMP END PARALLEL DO 
  !----------------------------------------------------------------
  !! main multiply, write loops
  outblock:do i3=1,nb
     if(add) then ! read previous result               
        call start_timer(startR)   !total cube , array size      ,f         , j ,esize,pointer location
        n1=size(res,1)
        er=sreed_window('gathers',3,(/n1,n2,n2/),(/n1,n2,blck/),(/0,0,(i3-1)*blck/),(/1,1,1/),8,res(1,1,1))
        if(er.lt.0) then
           write(0,*) 'sreed_window return=',er
           call seperr("Bad sreed_window")
        end if
        call end_timer(finishR,startR)
     else ! zero result array
        !write(0,*) 'zero result array'
        do k3=1,size(res,3)
           do k2=1,size(res,2)
              do k1=1,size(res,1)
                 res(k1,k2,k3)=cmplx(0.,0.)
              end do
           end do
        end do
     end if
     call start_timer(startMult)

     !----------------------------------------------------------------
     !! inner block (multiple panels before writing) 
     blocks:do i2=1,blck
              trace=i2+(i3-1)*i2
              ompi=1
!$OMP PARALLEL PRIVATE(ompi)
#ifdef SEP_OMP
              ompi = omp_get_thread_num()+1    
#endif
              do i=1,n1out !!select shot trace
                 shottrace(i,ompi)=conjg(cdata(i,trace))              
              end do
              !!  correlate with other traces
!$OMP DO PRIVATE(j2,j1,j3)
              x:do j2=1,n2
                 k1=0
                 t:do j1=nof,n1out
                    k1=k1+1
                    res(k1,j2,i2)=res(k1,j2,i2)+ shottrace(j1,ompi)*cdata(j1,j2)
                 end do t
              end do x
!$OMP END DO
!$OMP END PARALLEL

     end do blocks
     call end_timer(finishMult,startMult)
     !! writeout gathers
     call start_timer(startW)!total cube , array size      ,panel i3  , j ,esize,pointer location
     n1=size(res,1)
     er=srite_window('gathers',3,(/n1,n2,n2/),(/n1,n2,blck/),(/0,0,(i3-1)*blck/),(/1,1,1/),8,res(1,1,1))
     if(er.lt.0) then
        write(0,*) 'srite_window return=',er
        call seperr("Bad sreed_window")
     end if
     call end_timer(finishW,startW)

     er=((i3-1)*blck+i2-1)
     if(mod(er,n2/5)==0) write(0,'(2x,a,i5,a,i5,a,i3,a)') "   Finished",er,' of',n2,' panels:',100*er/n2,'%'
  end do outblock
  add=.true.
end do
  !----------------------------------------------------------------
  !!    Finish deallocate leave
  !        write(0,*) 'deallocate'
  deallocate(cdata,shottrace,res)
  call auxclose("gathers")
  call auxclose("traces")
  
  write(0,'(2x,a,f10.3,a)')'   Time reading       :', &
    finishR/60.,' (m)'
  write(0,'(2x,a,f10.3,a)')'   Time FFTing        :', &
    finishFFT,' (s)'
  write(0,'(2x,a,f10.3,a)')'   Time multing       :', &
    finishMult,' (s)'
  write(0,'(2x,a,f10.3,a)')'   Time writing       :', &
    finishW/60.,' (m)'
  call end_timer(finish,start)
  write(0,'(2x,a,f10.3,a)')'*  Total time       :', &
    finish/60.,' (m)'
  write(0,*) 'Program complete'
  call sep_end_prog()
  
end Program BigCorr


