program CabsDivide
  !! half the complex axis-- all the time you got comin!
  use sep
  use flags
!  use timers  
!  use cos_taper_mod
  use smooth_mod

  implicit none

  integer :: i1,i2,i3,i4,i5,n2,n3,n4,n5,nt,nf,nmt,er,ompn,ompi,job,rect2
  real::df,maxf,minf,min1,max1,of,ot,dt,d1,maxt,o2,o3,o4,o5,d2,d3,d4,d5
  complex,  pointer:: cdata(:,:)
  real,     pointer:: denom(:,:)
  real::startFFT,finishFFT,startW,finishW,scale,small
  logical::smoothme
  integer, external :: omp_get_thread_num,omp_get_num_threads,sreed_window,srite_window

  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("rect",rect2,3);call from_param('small',small,0.001)
  call from_param('smooth',smoothme,.false.)
  call from_aux("traces","n1",nf);   call from_aux("traces","n2",n2);  call from_aux("traces","n3",n3,1)
  call from_aux("traces","d1",df);   call from_aux("traces","d2",d2);  call from_aux("traces","d3",d3,1.)
  call from_aux("traces","o1",of,0.);call from_aux("traces","o2",o2);  call from_aux("traces","o3",o3,0.)
  call from_aux("traces","n4",n4,1); call from_aux("traces","n5",n5,1)
  call from_aux("traces","d4",d4,1.);call from_aux("traces","d5",d5,1.); 
  call from_aux("traces","o4",o4,0.);call from_aux("traces","o5",o5,0.); 
  call auxclose("traces")

  !write(0,*)'n5,n4,n3,n2,n1',n5,n4,n3,n2,nf
  job=(n5*n4*n3*n2)/5
  job=job/n2
  job=n2*job
  ompn=1; ompi=1

#ifdef SEP_OMP  
!$OMP PARALLEL 
     ompn = omp_get_num_threads()  
!$OMP END PARALLEL 
#endif


  write(0,'(2x,a,i5,a,i8,a)') 'FFT: ',n2,' traces x ',n3*n4*n5,' panels.'
  write(0,'(4x,a,i2,a,i5 )') 'OMP threads',ompn,'  Panels:',job*5
  write(0,*) 'Smooth:',smoothme,'Smooth rect (2axis)',rect2

  allocate(cdata(nf,n2))
  allocate(denom(nf,n2))

  call sep_begin_prog()
  call to_history("n1",nf,"decon");  call to_history("n2",n2,"decon");  call to_history("n3",n3,"decon")
  call to_history("d1",df,"decon");  call to_history("d2",d2,"decon");  call to_history("d3",d3,"decon")
  call to_history("o1",of,"decon");  call to_history("o2",o2,"decon");  call to_history("o3",o3,"decon")
  call to_history("n4",n4,"decon");  call to_history("n5",n5,"decon"); 
  call to_history("d4",d4,"decon");  call to_history("d5",d5,"decon"); 
  call to_history("o4",o4,"decon");  call to_history("o5",o5,"decon");
  call to_history("esize",8,"decon")

  if(smoothme) then
     call to_history('smoothrect',rect2,"decon")
  else
     call to_history('smooth',smoothme,"decon")
  end if

  sy:do i5=1,n5
     sx:do i4=1,n4
        gy:do i3=1,n3
           call sep_read(cdata,"traces")
!!$OMP PARALLEL PRIVATE(ompi)
!#ifdef SEP_OMP
!           ompi = omp_get_thread_num()+1    
!#endif
!!$OMP DO PRIVATE( i2,i1)
!$OMP PARALLEL DO PRIVATE( i2,i1)
           do i2=1,n2
              do i1=1,nf
                 denom(i1,i2)=cabs(cdata(i1,i2))
                 if(denom(i1,i2).lt.small) denom(i1,i2)=1
                 denom(i1,i2)=1./denom(i1,i2)
              end do
           end do
!$OMP END PARALLEL DO 
!!$OMP END DO 
!!$OMP END PARALLEL 

           if (smoothme) call smooth(denom,1,rect2,0,1,0,0,0)

!$OMP PARALLEL DO PRIVATE( i1,i2)           
           do i2=1,n2
              do i1=1,nf
                 cdata(i1,i2)=cdata(i1,i2)*denom(i1,i2)
              end do
           end do
!$OMP END PARALLEL DO 
           call sep_write(cdata,"decon")
           er=(i5-1)*n4+(i4-1)*n3+(i3)*n2         
           if(mod(er,job)==0) then
              write(0,'(2x,a,i5,a,i5,a,i3,a)')"   Finished",i3,' of',n3*n4*n5,' panels:',100*er/(n5*n4*n3*n2),'%'
           end if
        end do gy
     end do sx
  end do sy

  call auxclose("decon")
  deallocate(cdata,denom)
  call sep_end_prog()
end program CabsDivide
