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 , .false. )
                call from_param("apply"  , apply   , .false. )
                if(compute .and.         apply) then
                        call time_shift_params_init();  call time_shift_apply_params_init(); call time_shift_array_alloc();  call time_shift_apply_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_apply_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%ay%n,vgrid%ax%n)       :: mod1, mod2
        real, dimension(vgrid%az%n,vgrid%ay%n,vgrid%ax%n)       :: ts, cc
                if(compute) call  compute_time_shift(adj,add,mod1,mod2,ts,cc)
                if(apply  ) call  apply_time_shift(adj,add,mod1,mod2,ts,cc)
                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%ay%n,vgrid%ax%n)       :: mod1, mod2
        real, dimension(vgrid%az%n,vgrid%ay%n,vgrid%ax%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
                iref=iref+1; if(iref.gt.nmix)iref=1
                do icol=1,ncol
                        k=(iref-1)*ncol +icol
                        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 
                                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
                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
        end subroutine
!++++++++++++++++++++
!++++++++++++++++++++
        !apply time shifts
        subroutine apply_time_shift(adj,add,mod1,mod2,ts)
        logical                                                 :: adj, add
        real, dimension(vgrid%az%n,vgrid%ay%n,vgrid%ax%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

end module
