module joint_encode_data_mod
  use joint_encode_data_history
  use joint_encode_data_io
  use joint_encode_data_array
  use joint_encode_data_param 
  use joint_encode_data_wavelet
  use joint_encode_data_grid
  use joint_encode_data_compute_mod
  use joint_encode_adjnull !new
  use joint_encode_init    !new
  use joint_encode_outputs !new 
  use joint_encode_add     !new
  use joint_encode_fft     !new
  use joint_encode_regs    !new
  use joint_encode_prc     !new 
  use joint_encode_hlx     !new
  use joint_encode_nhlx    !new
  use cgstep_module
  use solver_sep_mod
  use chain_mod 

  
  implicit none

contains
!=======================================================================================
  integer function data_init() result(stat)
    call green_input_history()
    call green_params()                                     !initialize data modeling parameters 
    !call data_params()                                      !initialize data modeling parameters 
    call green_grid()                                       !get greens grid 
    call velocity_grid()                                    !get velocity/refl grid      
    call data_wavelet_s()                                    !get wavelet
    call data_params()                                      !initialize data modeling parameters 
    call data_array_alloc()                                 !allocate memory 
    call data_velocity()                                    !get velocity/refl models
    if (wfld_d_out) then;                                   !save down wavefield in time
        call wfld_output_history(wfld_d_); 
        call wfld_io_init(wfld_d_); 
    endif;  
    call data_grid_born();
    call data_output_history(); 
    call data_io_init();
    call inv_param()

    call joint_encode_params_init()
    call joint_encode_aux_init()
    call joint_array_alloc()
    call joint_encode_all()

    if (adj .or. inv) then
            if (inv) write(*,*)'Initializing aux files for inversion...........'
            if (adj) write(*,*)'Initializing aux files for migratoin...........'
            call joint_encode_inout_init()
    else
            write(*,*)'Initializing aux files for modeling........'
            call joint_encode_inout_model_init()
    endif
    stat = 0
  end function
!=======================================================================================
  integer function data_clos() result(stat); call data_array_dealloc(); call joint_encode_array_clos(); stat = 0; end function
!!=======================================================================================
  integer function data_oper() result(stat)
  use chain_mod 

    stat = encode_compute_init1(nodes)
    call encode_scode_init()
    if(inv) then
        if(mode==1)then
            call solver_tiny(m=modl0, d=data0, Fop=data_lop, stepper=cgstepper, niter=niter, verb=T,resd=norm, timer=T)
        elseif(mode==2)then
            call solver_regs(m=modl0, d=data0, Fop=data_lop, Aop=grad_x_temp_oper, stepper=cgstepper, nAop=nmod*nsurv, niter=niter, eps=eps, verb=T, timer=T)
        elseif(mode==3)then
            tresh=maxval(abs(data0))/100.
            call IRLS_init(tresh)
            call IRLS(x=modl0, dat=data0, oper=data_lop, solv=cgstepper, niter=niter, nfreq=nfreq, wght=weight_IRLS, err=norm)
        elseif(mode==4)then
            call solver_regs2(m=modl0, d=data0, Fop=data_lop, Aop1=grad_oper, Aop2=grad_temp_oper, stepper=cgstepper, nAop=nmod*nsurv, niter=niter, epsi=epsi0, lmbd=lmbd0, verb=T, timer=T)
        elseif(mode==5)then
            !call polydiv_init(nmod,aa)
            call solver_prc( m=modl0, d=data0, Fop=data_lop, Sop=int_oper, stepper=cgstepper, nSop=nmod*nsurv, niter=niter, eps=eps, verb=T)
        elseif(mode==6)then
            call solver_prc2( m=modl0, d=data0, Fop=data_lop, Sop=int_oper, stepper=cgstepper, nSop=nmod*nsurv, niter=niter, eps=epsi0, verb=T, timer=T)
        elseif(mode==7)then
            call solver_prc2( m=modl0, d=data0, Fop=data_lop, Sop=row_int, stepper=cgstepper, nSop=nmod*nsurv, niter=niter, eps=(epsi0*0.+1.), verb=T, timer=T)
        elseif(mode==8)then
            hlx_conv=2
            call filter_init()
            call solver_prc2( m=modl0, d=data0, Fop=data_lop, Sop=hlx_it,  stepper=cgstepper, nSop=nmod*nsurv, niter=niter, eps=epsi0, verb=T, timer=T)
        elseif(mode==9)then
            call  ns_init(dips0(1:nmod))
            call solver_prc2( m=modl0, d=data0, Fop=data_lop, Sop=nhlx_it, stepper=cgstepper, nSop=nmod*nsurv, niter=niter, eps=epsi0, verb=T, timer=T)
        elseif(mode==10)then
            call  ns_init(dips0(1:nmod))
            call solver_prc2( m=modl0, d=data0, Fop=data_lop, Sop=nhlx_row1, stepper=cgstepper, nSop=nmod*nsurv, niter=niter, eps=(epsi0*0.+1.), verb=T, timer=T)
        elseif(mode==11)then
            call  ns_init(dips0(1:nmod))
            call solver_prc2( m=modl0, d=data0, Fop=data_lop, Sop=nhlx_chain1, stepper=cgstepper, nSop=nmod*nsurv, niter=niter, eps=(epsi0*0.+1.), verb=T, timer=T)
        endif
    else 
         stat=data_lop(adj,F,modl0,data0);
    end if 
    call joint_encode_output()

  end function
!!=======================================================================================
  integer function data_lop(adj, add, modl, data) result(stat)
    logical,intent(in)    :: adj,  add 
    real, dimension(:)    :: modl, data 
    integer                 :: isurv
        stat = 1
        call joint_adjnull(adj,add,modl,data)
        if(verb>5)then
            write(0,*);write(0,*)'+++++'
            write(0,'(A10,10F25.5)')'min',minval(data),minval(data1),minval(data2),minval(modl),minval(modl1),minval(modl2)
            write(0,'(A10,10F25.5)')'max',maxval(data),maxval(data1),maxval(data2),maxval(modl),maxval(modl1),maxval(modl2)
            write(0,*)'+++++';write(0,*)
        endif
        do isurv=1,nsurv;
            if(isurv==1)then; if(verb>4)write(0,*)'Starting survey 1: adj|dat|mod',adj,sum(data1),sum(modl1);stat = encode_compute(adj,modl1,data1,1); if(verb>5)write(0,*)'sum dat1:',sum(data1),'sum mod1:',sum(modl1);endif
            if(isurv==2)then; if(verb>4)write(0,*)'Starting survey 2: adj|dat|mod',adj,sum(data2),sum(modl2);stat = encode_compute(adj,modl2,data2,2); if(verb>5)write(0,*)'sum dat2:',sum(data2),'sum mod2:',sum(modl2);endif
            if(isurv==3)then; if(verb>4)write(0,*)'Starting survey 3: adj|dat|mod',adj,sum(data3),sum(modl3);stat = encode_compute(adj,modl3,data3,3); if(verb>5)write(0,*)'sum dat3:',sum(data3),'sum mod3:',sum(modl3);endif
            if(isurv==4)then; if(verb>4)write(0,*)'Starting survey 4: adj|dat|mod',adj,sum(data4),sum(modl4);stat = encode_compute(adj,modl4,data4,4); if(verb>5)write(0,*)'sum dat4:',sum(data4),'sum mod4:',sum(modl4);endif
            if(isurv==5)then; if(verb>4)write(0,*)'Starting survey 5: adj|dat|mod',adj,sum(data5),sum(modl5);stat = encode_compute(adj,modl5,data5,5); if(verb>5)write(0,*)'sum dat5:',sum(data5),'sum mod5:',sum(modl5);endif
            if(isurv==6)then; if(verb>4)write(0,*)'Starting survey 6: adj|dat|mod',adj,sum(data6),sum(modl6);stat = encode_compute(adj,modl6,data6,6); if(verb>5)write(0,*)'sum dat6:',sum(data6),'sum mod6:',sum(modl6);endif
        enddo
        call joint_add(adj,add,modl,data)

        stat=0
  end function
!!=======================================================================================
  integer function encode_compute(adj,  modl, data, isurv) result(stat)
    logical,intent(in) :: adj 
    real, dimension(dgrid%agx%n, dgrid%agy%n,dgrid%asx%n,dgrid%asy%n,dgrid%aw%n):: data 
    real, dimension(vgrid%ax%n, vgrid%ay%n, vgrid%az%n):: modl
    real, dimension(vgrid%ax%n, vgrid%ay%n, vgrid%az%n,ggrid%aw%n ):: model2
    integer         :: iiw, iiiw, nb_nw, max_nw, isurv, i
            stat = encode_beg(adj,data,modl)
            do i=1,ggrid%aw%n ; model2(:,:,:,i)=modl; enddo
            stat=encode_compute_init_s(isurv)
            nb_nw=nint(1.*ggrid%aw%n/nodes)  
            do iiw= 1,ggrid%aw%n,nodes
                max_nw=min(iiw+nodes-1,ggrid%aw%n)
                stat = encode_compute_init_blk(iiw,max_nw,isurv);
!$OMP parallel do private(inode, iw, iiiw, w, index_aw) schedule(dynamic,1) 
                do iw=iiw, max_nw                                ! freq.  blocks
                    inode=omp_get_thread_num()+1;   
                    iiiw=iw-iiw+1
                    w=ggrid%aw%o+(iw-1)*ggrid%aw%d; index_aw=1.5+(w-dgrid%aw%o)/dgrid%aw%d
                    if(sum(cabs(csou_full(:,:,index_aw)))<dgrid%asx%n*dgrid%asy%n*mcsou) cycle; 
                    call encode_compute_data(adj, model2(:,:,:,index_aw), fdata(:,:,:,:,index_aw),inode,iiiw,index_aw, isurv); 
                    !write(0,'(6I10,8F40.5)')inode,iw,iiw,iiiw,max_nw,index_aw, sum(fdata(:,:,:,:,index_aw))/maxval(cabs(fdata(:,:,:,:,index_aw))), sum(model2(:,:,:,index_aw))/maxval(model2(:,:,:,index_aw)) 
                    stat = encode_report(index_aw,w,sum(model2,4),cabs(fdata(:,:,:,:,index_aw)),adj,inode,isurv)
                end do
!$OMP end parallel do 
            end do 
            if(adj)modl=sum(model2,4)
        stat = encode_end(adj,data,modl)
  end function
!!=======================================================================================
  integer function  encode_beg(adj,data,modl) result(stat)
    logical,intent(in) :: adj
    real, dimension(dgrid%agx%n, dgrid%agy%n,dgrid%asx%n,dgrid%asy%n,dgrid%aw%n):: data
    real, dimension(vgrid%ax%n, vgrid%ay%n, vgrid%az%n):: modl

        stat=1
        if(adj.or.inv)then
            fdata=0.; fdata=data; 
            stat=encode_transf(1,fdata)
        else
            fdata=0.; 
        endif
        stat=0
  end function
!!=======================================================================================
  integer function  encode_end(adj,data,modl) result(stat)
    logical,intent(in) :: adj
    real, dimension(dgrid%agx%n, dgrid%agy%n,dgrid%asx%n,dgrid%asy%n,dgrid%aw%n):: data
    real, dimension(vgrid%ax%n, vgrid%ay%n, vgrid%az%n):: modl
        stat=1
        if(adj .and. .not. inv)then
            call data_io_img(modl)
        elseif (.not. (adj .or. inv))then
            call data_io_a('rite', fdata);
            stat=encode_transf(-1,fdata)
            data=fdata
        elseif(inv)then
            stat=encode_transf(-1,fdata)
            data=fdata
        endif 
        stat=0
    end function
!!=======================================================================================
  integer function  encode_init_data(data) result(stat)
    real, dimension(dgrid%agx%n, dgrid%agy%n,dgrid%asx%n,dgrid%asy%n,dgrid%aw%n):: data
        stat=1
        !stat=encode_transf(-1,fdata); data=fdata
        stat=0
  end function
!!=======================================================================================
  integer function lop(adj, add, modl, data) result(stat) !used for dottest
    logical,intent(in)    :: adj,  add
    real, dimension(:)    :: modl, data
    complex, dimension(size(data)):: ddata
    integer                 :: isurv
        stat = 1
        call adjnull(adj,F,modl,data)
        do isurv=1,nsurv
            if(isurv==1) stat = encode_compute(adj,modl1,data1,1)
            if(isurv==2) stat = encode_compute(adj,modl2,data2,2)
            if(isurv==3) stat = encode_compute(adj,modl3,data3,3)
            if(isurv==4) stat = encode_compute(adj,modl4,data4,4)
            if(isurv==5) stat = encode_compute(adj,modl5,data5,5)
            if(isurv==6) stat = encode_compute(adj,modl6,data6,6)
        enddo
        call joint_add(adj,add,modl,data)
        stat=0
  end function
!!=======================================================================================!!=======================================================================================
  integer function  encode_report(index_aw,w,modl,data,adj,inode,isurv) result(stat)
    real, dimension(vgrid%ax%n, vgrid%ay%n, vgrid%az%n):: modl
    real, dimension(dgrid%agx%n, dgrid%agy%n,dgrid%asx%n,dgrid%asy%n,dgrid%aw%n):: data
    integer                 :: index_aw, inode, isurv
    real                    ::  w
    logical,intent(in)      :: adj
    stat=1
    if(verb>2)then
        if(inv)then
            if(verb>3)then
                if(adj)then
                    write(0,'(A10,I5,A3,F7.3,A5,2F10.4, A6,I5, A15, E20.5, A7, I3)')'done iw: ', index_aw, '|', w,'=> (',sum(csou0(:,:,index_aw,isurv))/size(csou0(:,:,index_aw,isurv)),') of ',ggrid%aw%n, ' |  sum mig:', sum(modl), ' @ node', inode; 
                else
                    write(0,'(A10,I5,A3,F7.3,A5,2F10.4, A6,I5, A15, E20.5, A7, I3)')'done iw: ', index_aw, '|', w,'=> (',sum(csou0(:,:,index_aw,isurv))/size(csou0(:,:,index_aw,isurv)),') of ',ggrid%aw%n, ' |  sum dat:', sum(data), ' @ node', inode; 
                endif
            endif
        else
            if(adj)then
                write(0,'(A10,I5,A3,F7.3,A5,2F10.4, A6,I5, A15, E20.5, A7, I3)')'done iw: ', index_aw, '|', w,'=> (',sum(csou0(:,:,index_aw,isurv))/size(csou0(:,:,index_aw,isurv)),') of ',ggrid%aw%n, ' |  sum mig:', sum(modl), ' @ node', inode; 
            else
                write(0,'(A10,I5,A3,F7.3,A5,2F10.4, A6,I5, A15, E20.5, A7, I3)')'done iw: ', index_aw, '|', w,'=> (',sum(csou0(:,:,index_aw,isurv))/size(csou0(:,:,index_aw,isurv)),') of ',ggrid%aw%n, ' |  sum dat:', sum(data), ' @ node', inode; 
            endif
        endif
    endif
    stat=0
  end function
!!=======================================================================================
  integer function encode_io(modl, data) result(stat)
    real, dimension(dgrid%agx%n, dgrid%agy%n,dgrid%asx%n,dgrid%asy%n,dgrid%aw%n):: data
    real, dimension(vgrid%ax%n, vgrid%ay%n, vgrid%az%n):: modl
        stat = 1
        if (adj.or.inv)then
            call data_io_img(modl)  
        else
            call data_io_a('rite', fdata);
        endif
        stat = 0
  end function
!===========================================================================================================================================================================
end module
