Program LineCor
  !##################################################
  !#                                              
  !#  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 

  implicit none
  integer :: j2,j3,i1,i2,i3,i4,i5,a,b,c,d,er,shy,shx,nf,imp
  integer :: n1,n2,n3,n4,n5,esize,n1out,n4out,n5out,ompn, ompi
  real    :: o1,o2,o3,d1,d2,d3,maxt,mem
  real   , pointer,dimension(:,:,:)  :: data
  complex, pointer,dimension(:,:,:)  :: cdata
  complex, pointer,dimension(:,:)  :: res,tmpc
  real,    pointer,dimension(:,:)  :: tmpr
  complex, pointer,dimension(:,:)    :: shottrace
  logical:: mig,time,center,add
  integer*8,pointer:: planF(:,:),planI(:)
  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(5)::cube,panel,fwind,jwind
  real:: maxf,minf,scale,min1,max1,df
  logical:: signum

  startR=0.;finishR=0.;startW=0.;finishW=0.;start=0.;finish=0.
  startFFT=0.;finishFFT=0.; startMult=0.; finishMult=0.;
  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_param('filenum',er,9999)
  call from_param('maxf',maxf,.5/d1); call from_param('minf',minf,0.)
  scale=0.075*(maxf-minf)
  call from_param('max1',max1,maxf-scale); call from_param('min1',min1,minf+scale)
  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)
  call from_param('j2',j2,n2/2);  call from_param('j3',j3,n2/2);
  call auxclose("traces")
  ompn=1
  ompi=1

  df=1./(d1*n1)
  n1out=n1/2; n4out=n2/j2+1 ; n5out=n3/j3+1
  nf=n1/2+1
  write(0,*) 'n1,n2,n3,nf,n1out',n1,n2,n3,nf,n1out
  write(0,*) 'n1,n4,n5',n1out,n4out,n5out
  a=1; b=n3
  if(j3.ge.n3/2) then
     a=n3/2
     b=a
     n5out=1
  end if
  c=1; d=n2
  if(j2.ge.n2/2) then
     c=n2/2
     d=c
     n4out=1
  end if

  if(add) call auxinout("gathers")
  write(0,*) 'file number',er
  call to_history('filenum',er,"gathers")
  call to_history('esize',4,"gathers");call to_history('d1',d1,"gathers")
  call to_history('n1',n1out,"gathers");call to_history('o1',o1,"gathers")
  call to_history('n2',n2,"gathers");call to_history('o2',o2,"gathers");call to_history('d2',d2,"gathers")
  call to_history('n3',n3,"gathers");call to_history('o3',o3,"gathers");call to_history('d3',d3,"gathers")
  call to_history('n4',n4out,"gathers");call to_history('o4',o2,"gathers");call to_history('d4',d2*j2,"gathers")
  call to_history('n5',n5out,"gathers");call to_history('o5',o3,"gathers");call to_history('d5',d3*j3,"gathers")
  if(c.eq.d) then
  	call to_history('n4',1,"gathers");
	call to_history('o4',o2+(n2/2-1)*d2);
	call to_history('d4',1.,"gathers")
  end if
  if(a.eq.b) then 
     call to_history('n5',1,"gathers");
     call to_history('o5',o3+(n3/2-1)*d3);;
     call to_history('d5',1.,"gathers")
  end if
  call to_history('label3','sy');   call to_history('label4','sx'); 
  call to_history('label2','gy');   call to_history('label2','gx'); 

  call sep_begin_prog()
  mem=0.
#ifdef SEP_OMP  
!$OMP PARALLEL 
  ompn = omp_get_num_threads()  
!$OMP END PARALLEL 
#endif
  call init_cos_taper(nf,0.,df,minf,min1,maxf,max1,ompn)
  !----------------------------------------------------------------
  !!     Read traces section
  allocate(data(n1,n2,n3))
  mem=mem+4*size(data)
  write(0,*) 'allocate data',mem/1024/1024
  call sep_read(data,"traces")
  if(signum) then
     where(data.ne.0.) 
        data=data/abs(data)
     end where
  end if
  allocate(cdata(nf,n2,n3))!,cdata2(n1,n2))
  mem=mem+8*size(cdata)
  write(0,*) 'allocate cdata',mem/1024/1024
  allocate (planF(n2,n3),planI(n2))
  mem=mem+8*size(planF)+8*size(planI)
  write(0,*) 'allocate FFT plans',mem/1024/1024
   allocate(tmpr(n1,n2))
   allocate(tmpc(nf,n2))
 !----------------------------------------------------------------
  !!    FFT Section
  !  call start_timer(startFFT)
  write(0,*) ' planing FFTs'
   do i3=1,n3
      do i2=1,n2
         call sfftw_plan_dft_r2c_1d(planF(i2,i3),n1,&
         data(1,i2,i3),cdata(1,i2,i3),FFTW_PATIENT)        
      end do
   end do
write(0,*) 'fft data'
!$OMP PARALLEL DO PRIVATE( i3,i2)
   do i3=1,n3
      do i2=1,n2
         call sfftw_execute(planF(i2,i3))
      end do
   end do
!$OMP END PARALLEL DO 
   deallocate(data)
write(0,*) 'done fft'
   do i2=1,n2
      call sfftw_plan_dft_c2r_1d(planI(i2),n1,&
      tmpc(1,i2),tmpr(1,i2),FFTW_PATIENT)   
   end do


  allocate(shottrace(nf,ompn))
  mem=mem+8*size(shottrace)
  write(0,*) 'allocate shottrace',mem/1024/1024
  shottrace=cmplx(0.,0.)
  allocate(res(n1out,n2))
  mem=mem+8*size(res)
  write(0,*) 'allocate res',mem/1024/1024
  cube=(/n1out,n2,n3,n4out,n5out/)
  panel=(/n1out,n2,1,1,1/)
  jwind=(/1,1,1,1,1/)
  write(0,*) 'a,b',a,b
  write(0,*) 'c,d',c,d
  scale=1./real(n1)
  
  !----------------------------------------------------------------
  !! main multiply, write loops
  shx=-1
  sx:do i5=a,b,j3
     shx=shx+1; shy=-1
write(0,*) 'shoty',shx+1
     sy:do i4=c,d,j2
        shy=shy+1
write(0,*) '   shotx',shy+1
        do i1=1,nf !!select shot trace
write(0,*) '     i1',i1
           shottrace(i1,:)=conjg(cdata(i1,i4,i5))              
        end do
        gx:do i3=1,n3
write(0,*)'  panel',i3
           fwind=(/0,0,i3-1,shy,shx/)
           if(add) then ! read previous result               
              er=sreed_window('gathers',5,cube,panel,fwind,jwind,4,res(1,1))
              if(er.lt.0) then
                 write(0,*) 'sreed_window return=',er
                 call seperr("Bad sreed_window")
              end if
           else ! zero result array
write(0,*) 'zero'
!$OMP PARALLEL DO PRIVATE( i2,i1)
              do i2=2,size(res,2)
                 do i1=1,size(res,1)
                    res(i1,i2)=0.
                 end do;
              end do;
!$OMP END PARALLEL DO 
           end if !!add


write(0,*) 'correlate/iFFT'
           !!  correlate with other traces
!$OMP PARALLEL PRIVATE(ompi)
#ifdef SEP_OMP
           ompi = omp_get_thread_num()+1    
#endif
!$OMP DO PRIVATE( i1,i2)
           gy:do i2=1,n2
              t:do i1=1,nf
                 tmpc(i1,i2)= shottrace(i1,ompi)*cdata(i1,i2,i3)*cos_taper_filter(i1,ompi)*scale
              end do t
              call sfftw_execute(planI(i2))!!makes tmpr from tmc
              t2:do i1=1,n1out
                 res(i1,i2)=res(i1,i2)+tmpr(i1,i2)
              end do t2
           end do gy
!$OMP END DO 
!$OMP END PARALLEL 

write(0,*) 'writing'
           !! writeout gathers
           er=srite_window('gathers',5,cube,panel,fwind,jwind,4,res(1,1))
           if(er.lt.0) then
              write(0,*) 'srite_window return=',er
              call seperr("Bad sreed_window")
           end if
write(0,*) 'done '
        end do gx
        write(0,*) "Finished shot",shy+shx+1
     end do sy
  end do sx
  !----------------------------------------------------------------
  !!    Finish deallocate leave
  write(0,*) 'deallocate'
  deallocate(cdata,shottrace,res)
  call auxclose("gathers")
  !  write(0,'(8x,a,f10.3,a)')' *       Time FFTing       :', &
  !  finishFFT,' (s)'
  !        write(0,'(8x,a,f10.3,a)')' *       Time reading       :', &
  !         finishR/60.,' (m)'
  !        write(0,'(8x,a,f10.3,a)')' *       Time multing       :', &
  !         finishMult,' (s)'
  !        write(0,'(8x,a,f10.3,a)')' *       Time writing       :', &
  !         finishW/60.,' (m)'
        call end_timer(finish,start)
        write(0,'(8x,a,f10.3,a)')' *  Total time       :', &
           finish/60.,' (m)'
  !        write(0,*) 'Program complete'
  call sep_end_prog()
end Program LineCor
