module time_shift_mod
        use time_shift_params
        use time_shift_utils
        implicit none

        contains
!====================
      integer function time_shift_init() result(stat)
                call from_param("compute", compute , .true. )
                call from_param("apply"  , apply   , .true. )
                if(compute .and.         apply) then
                        call time_shift_params_init();  call time_shift_apply_params_init(); call time_shift_array_alloc()
                elseif(compute .and. (.not. apply)) then
                        call time_shift_params_init();  call time_shift_array_alloc()
                elseif(.not. compute .and. apply)then
                                                        call time_shift_apply_params_init(); call time_shift_array_alloc()
                else
                        call seperr('You dont know what you want to do??? compute or apply must be .true.')  
                endif        
                stat=0
        end function
!++++++++++++++++++++
        integer function time_shift_clos() result(stat)
                if(compute) call time_shift_array_dealloc()
                if(apply  ) call time_shift_apply_array_dealloc()
                stat=0
        end function
!++++++++++++++++++++
        integer function time_shift_oper() result(stat)
        logical                                                 :: adj, add
        real, dimension(vgrid%az%n,vgrid%ax%n,vgrid%ay%n)       :: mod1, mod2
        real, dimension(vgrid%az%n,vgrid%ax%n,vgrid%ay%n)       :: ts, cc
                !if(compute) call  compute_time_shift(adj,add,mod1,mod2,ts,cc,isurv)
                !if(apply  ) call    apply_time_shift(adj,add,mod1,mod2,ts   ,isurv)
                stat=0
        end function
!++++++++++++++++++++
        !compute time shifts
        subroutine compute_time_shift(adj,add,mod1,mod2,ts,cc)
        logical                                                 :: adj, add
        real, dimension(vgrid%az%n,vgrid%ax%n,vgrid%ay%n)       :: mod1, mod2
        real, dimension(vgrid%az%n,vgrid%ax%n,vgrid%ay%n)       :: ts, cc
        integer                                                 :: time_array_0(8), time_array_1(8)
        real                                                    :: start_time, finish_time 
        integer                                                 :: stat 
        stat = 1

        iref=0; nhalf=nmix/2
        !intialize the fisrt 1 to nmix-1 rows
        do irow=-nhalf+1,nhalf
                write(0,*)'irow',irow
                iref=iref+1; if(iref.gt.nmix)iref=1
                do icol=1,ncol
                write(0,*)'icol',icol
                        k=(iref-1)*ncol +icol
                write(0,*)'k',k
                        if((irow.lt.1.or.irow.gt.nrow).or.(icol.lt.1.or.icol.gt.ncol)) then
                                iexist(k)=0                                             !does not exist??
                        else
                                iexist(k)=1                                             !it exists...
                        end if
                        if(iexist(k).ne.0) then                                         !if trace(s) exists
                                trace1(:,k)=mod1(:,irow,icol);  call apply_agc(trace1(:,k), trace3, vgrid%az%n, nagc);                  !read and apply agc to base
                                trace2(:,k)=mod2(:,irow,icol);  call apply_agc(trace2(:,k), trace3, vgrid%az%n, nagc);                  !read and apply agc to moni
!write(0,*)'=========',sum(mod1),sum(mod2)
 !                       do iz=1,vgrid%az%n,10
  !                              write(0,*)irow,iz,trace2(iz,k),trace1(iz,k)
   !                     enddo
!write(0,*)'++++++++'
!stop
                                do iout=nskip,vgrid%az%n,nskip                                                                          !loop through the trace 
                                        iskip=iout/nskip                                                                                !skip nskip samples
                                        if(iout.gt.(mlag+igate+1+ilag_bias).and.(iout+mlag+igate+ilag_bias).lt.vgrid%az%n) then         !sample within the limits??
                                                do ilag=-mlag,mlag                                                                      !correlation lag...from -mlag to mlag 
                                                        index_tcor=iskip+(ilag+mlag+1-1)*IMAX + (k-1)*IMAX*JMAX                         !compute index in the BIG tccor trace
                                                        tccors(index_tcor) = 0.; den1=0.;den2=0; tmp_ccor=0.;                           !initialize
                                                        do i1=iout-igate, iout+igate                                                    !correlation loop from -igate to igate (around sample iout)  
                                                                i2=i1+ilag + ilag_bias                                                  !sample in trace 2
                                                                tmp1=trace1(i1,k); tmp2=trace2(i2,k)                                    !initilize...saves reading vector
                                                                wt=(igate - abs(1.*(i1-iout)) )/(1.*igate); if(.not. lwght) wt=1        !compute weight       
                                                                tmp_ccor = tmp_ccor + tmp1*tmp2*wt                                      !compute correlation
                                                                den1     = den1     + tmp1*tmp1*wt                                      !compute denom. 1
                                                                den2     = den2     + tmp2*tmp2*wt                                      !compute denom. 2
                                                        end do
                                                        den1=sqrt(den1)*sqrt(den2); if(den1.ne.0) tccors(index_tcor)=tmp_ccor/den1      !normalize correlation
                                                end do
                                        end if
                                end do
                        end if
                end do
        end do
        do irow=1,nrow
                if(verb>6)write(0,*)'Working on row',irow,' of ',nrow
                iref=iref+1; if(iref.gt.nmix)iref=1
                do icol=1,ncol
                        k=(iref-1)*ncol +icol
                        !read in row irow+nhalf for each col
                        if((irow+nhalf.lt.1.or.irow+nhalf.gt.nrow).or.(icol.lt.1.or.icol.gt.ncol)) then
                                iexist(k)=0                           !does not exist
                        else
                                iexist(k)=1                           !exists...
                        end if
                        if(iexist(k).ne.0) then
                                iexist(k)=1;
                                trace1(:,k)=mod1(:,irow,icol);  call apply_agc(trace1(:,k), trace3, vgrid%az%n, nagc); !read and apply agc to base
                                trace2(:,k)=mod2(:,irow,icol);  call apply_agc(trace2(:,k), trace3, vgrid%az%n, nagc); !read and apply agc to moni
                                do iout=nskip,vgrid%az%n,nskip
                                        iskip=iout/nskip
                                        if(iout.gt.(mlag+igate+1+ilag_bias).and.(iout+mlag+igate+ilag_bias).lt.vgrid%az%n) then
                                                do ilag=-mlag,mlag
                                                        index_tcor=iskip +(ilag+mlag+1-1)*IMAX + (k-1)*IMAX*JMAX
                                                        tccors(index_tcor)=0.; den1=0.; den2=0.; tmp_ccor=0.
                                                        do i1=iout-igate, iout+igate
                                                                i2=i1+ilag + ilag_bias
                                                                tmp1=trace1(i1,k); tmp2=trace2(i2,k)
                                                                wt=(igate - abs(1.*(i1-iout)) )/(1.*igate); if(.not. lwght) wt=1
                                                                tmp_ccor = tmp_ccor + tmp1*tmp2*wt
                                                                den1     = den1     + tmp1*tmp1*wt
                                                                den2     = den2     + tmp2*tmp2*wt
                                                        end do
                                                        den1=sqrt(den1)*sqrt(den2); if(den1.ne.0) tccors(index_tcor)=tmp_ccor/den1
                                                end do
                                        end if
                                end do
                        end if
                end do

                do icol=1,ncol                                                                                                  !for all columns in this row 
                        do iout=nskip,vgrid%az%n,nskip                                                                          !loop thru all samples 
                                ccmax=-1; iskip=iout/nskip; ntraces=0; tshift=-ilag_bias*vgrid%az%d;exit1=F                     !skip nkip samples for output
                                if(iout.ge.(mlag+igate+1).and.(iout+mlag+igate).le.vgrid%az%n) then                             !sample within limits??
                                        do ilag=-mlag,mlag                                                                      !loop thru the lags..-malg to mlag
                                                iilag=mlag+1+ilag; ccval(iilag)=0; nsum=0                                       !actual index within correlation series for that point
                                                do imix=1,nmix                                                                  !row loop for mixing
                                                        do icave=icol-nhalf,icol+nhalf                                          !col loop for mixing
                                                                if(icave.lt.1 .or. icave.gt.ncol) cycle                         
                                                                k=(imix-1)*ncol + icave
                                                                if(iexist((imix-1)*ncol + icave).ne.0) then                     !trace exists??
                                                                        index_tcor=iskip +(iilag-1)*IMAX + (k-1)*IMAX*JMAX      !index in BIG tccor trace        
                                                                        ccval(iilag)=ccval(iilag) + tccors(index_tcor)          !add up nmix*nmix traces 
                                                                        nsum=nsum+1; ntraces=1
                                                                end if
                                                        end do
                                                end do
                                                if(nsum.ne.0) ccval(iilag)=ccval(iilag)/nsum                                    !average the mixed ccors values
                                        end do
                                        if(ntraces.eq.0) then; exit1=T; exit; endif                     
                                        call mxcval(ccval,2*mlag+1,tmax,ccmax)
                                        if(tmax.ge.0)tshift=(tmax-mlag-1)*vgrid%az%d                                            !compute time-shift
                                endif
                                trace3(iout)=tshift+ilag_bias*vgrid%az%d                                                        !adust by lag bias
                                if(lccor) trace4(iout)=ccmax                                                                    !store max. ccor
                        end do
                        if(exit1==T) cycle
                        if(iskip.ne.1) then
                                do iout=1,vgrid%az%n
                                        i1=iout/nskip  ; i1=nskip*i1
                                        i2=iout/nskip+1; i2=nskip*i2
                                        if(i1.lt.1) then
                                                trace3(iout)=0.; if(lccor) trace4(iout)=-1
                                        else if(i2.gt.vgrid%az%n) then
                                                trace3(iout)=0.; if(lccor) trace4(iout)=-1
                                        else
                                                          trace3(iout)=trace3(i1) + ((trace3(i2)-trace3(i1))*(iout-i1))/(i2-i1) !interpolate 
                                                if(lccor) trace4(iout)=trace4(i1) + ((trace4(i2)-trace4(i1))*(iout-i1))/(i2-i1) 
                                        end if
                                end do
                        end if
                        ts(:,irow,icol)=trace3
                        cc(:,irow,icol)=trace4
                end do
        end do
        write(0,*)'Done'
        call sep_write(mod1,ts1_)
        call sep_write(mod2,ts2_)
        call sep_write(ts,cc1_)
        call sep_write(cc,cc2_)

        stop
        end subroutine
!++++++++++++++++++++
!++++++++++++++++++++
        !apply time shifts
        subroutine apply_time_shift(adj,add,mod1,mod2,ts)
        logical                                                 :: adj, add
        real, dimension(vgrid%az%n,vgrid%ax%n,vgrid%ay%n)       :: mod1, mod2, ts    
        integer                                                 :: time_array_0(8), time_array_1(8)
        real                                                    :: tsamp
        real                                                    :: start_time, finish_time

        do irow=1,nrow
                if(irow==1.or.modulo(irow,nrow/50)==0) write(0,*)'Working on row',irow,' of ',nrow     
                do icol=1,ncol
                        trace_tsft=   ts(:,irow,icol); 
                        trace_uvol=mod2(:,irow,icol); 
                        do iout=1, vgrid%az%n
                                tsamp=iout+trace_tsft(iout)/(1.*vgrid%az%d);                    !compute num. of time-shift sample
                                call cubic(trace_uvol, tsamp, trace_svol(iout), vgrid%az%n)     !interpolate
                        end do
                        mod2(:,irow,icol)=trace_svol                                           !write shifted volume
                end do
        end do
        end subroutine
!#####################################################  conv or deconv
        integer function shift_it(adj, add, xx, yy) result(stat)
                logical, intent(in)             ::  adj, add 
                real   , dimension (:), target  ::  xx , yy
                real   , dimension (nmod)       :: m, m0, d, d0, ts, cc
                integer                         :: isurv
                !call joint_adjnull2 (adj, add, xx, yy)
        
               

                m0=xx(1:nmod)
                d0=yy(1:nmod)


                do isurv=2, nsurv
                          m=xx ((isurv-1)*nmod+1:isurv*nmod)
                          d=yy ((isurv-1)*nmod+1:isurv*nmod)
                         ts=ts0((isurv-1)*nmod+1:isurv*nmod)
                         cc=cc0((isurv-1)*nmod+1:isurv*nmod)
                
                        if(adj)then
                                call time_shift_it(adj, add, m0, m, ts, cc, isurv); 
                        else
                                call time_shift_it(adj, add, d0, d, ts, cc, isurv); 
                        endif
                        call shift_transp(-1,m ,d ); xx((isurv-1)*nmod+1:isurv*nmod)=m; yy((isurv-1)*nmod+1:isurv*nmod)=d
                        call shift_transp(-1,m0,d0);
                enddo
                if(ltsft)call shift_transp(-1,ts);  ts0((isurv-1)*nmod+1:isurv*nmod)=ts; 
                if(lccor)call shift_transp(-1,cc);  cc0((isurv-1)*nmod+1:isurv*nmod)=cc; 
                stat=0
        end function
!++++++++++++++++++++
        subroutine time_shift_it(adj,add,mod1,mod2,ts,cc,isurv) 
        logical                                                 :: adj, add
        real, dimension(vgrid%az%n,vgrid%ax%n,vgrid%ay%n)       :: mod1, mod2
        real, dimension(vgrid%az%n,vgrid%ax%n,vgrid%ay%n)       :: ts, cc
        integer                                                 :: isurv
        optional                                                :: isurv
        write(0,*)'eerer'
                if(compute) call  compute_time_shift(adj,add,mod1,mod2,ts,cc)
        write(0,*)'eerer222'
                if(apply  ) call    apply_time_shift(adj,add,mod1,mod2,ts   )
                if(present(isurv).and.ltsft)call ts_io2(ts,isurv);
                if(present(isurv).and.lccor)call cc_io2(cc,isurv);
        end subroutine
!#####################################################  transp before filter and back
        subroutine shift_transp(fwd,xx,yy)
                integer,intent(in)                      :: fwd 
                real, dimension (:)                     :: xx, yy 
                real, dimension(:,:), allocatable       :: temp
                optional                                :: yy
                
                if (fwd==1)then
                        allocate (temp(vgrid%ax%n,vgrid%az%n))
                        temp=reshape(xx,(/vgrid%ax%n,vgrid%az%n /)); xx=reshape(transpose(temp),(/ nmod /))
                        if(present(yy))then; temp=reshape(yy,(/vgrid%ax%n,vgrid%az%n /)); yy=reshape(transpose(temp),(/ nmod /)); endif
                else
                        allocate (temp(vgrid%az%n,vgrid%ax%n))
                        temp=reshape(xx,(/vgrid%az%n,vgrid%ax%n /)); xx=reshape(transpose(temp),(/ nmod /))
                        if(present(yy))then; temp=reshape(yy,(/vgrid%az%n,vgrid%ax%n /)); yy=reshape(transpose(temp),(/ nmod /));endif
                endif
                deallocate (temp)
        end subroutine

end module
