program WEItst_main  
  !! 1 ref slowness   in the background wavefield (slo1 & mwc1)
  !! square root scattering operator with 1 mean reference slowness (bor1)
  use sep 
  use mpi_sep 
  use wei_util 
  use wei_par 
  use wei_pro 
  use wei_tag 
  use wei_mva

  use weiop_cig 
  use weiop_mwc 
  use weiop_ssf
  use weiop_slo
  use weiop_wem 
  use weiop_bor

  use oc_global_mod 
  use oc_utils_mod 
  use oc_cgstep_mod 
  use oc_laplacian_mod
  use oc_sos_mod

  use oc_tst_solver_mod 
  use oc_mpi_solverreg_mod
  use oc_mpi_solverpre_mod
  use oc_mpi_dottest_mod 

  implicit none
  integer :: stat
  integer :: mpierr 

  real               :: ocmaxmem, oceps
  integer            :: ocniter
  logical            :: ocverb,ocmovie,ocwght 
  integer            :: ocsolver
  integer            :: ocnlap

  character(len=128) :: res_, name 
  character(len=128) :: y_

  integer                          :: nf
  type(fileinfo)                   :: file
  integer,dimension(:),allocatable :: nnew
  !----------------------------------------------------------------
  call MPI_INIT(mpierr)
  call sep_init()
  !----------------------------------------------------------------
  !! init params for MPI and OMP
  call weipar_init()
  call weipar_report()

  if(wpar%impi==0) then
     call MPI_SEP_SEND_ARGS(wpar%nmpi,1,0)
  else
     call MPI_SEP_RECEIVE_ARGS()
  end if

  !! init tag names on all nodes
  call weitag_init()

  !! init process parameters on the MASTER
  if(wpar%impi==0) call weipro_init()
  !! broadcast parameters from MASTER to NODES
  if(wpar%nmpi >0) call weipro_bcast()
  !----------------------------------------------------------------
  if(wpar%impi==0) then
     call from_param("ocniter",ocniter,1)

     call from_param("ocmaxmem",ocmaxmem,     0.)
     !! convert Mbytes to bytes
     ocmaxmem = ocmaxmem * 1024*1024

     call from_param("oceps",   oceps,0.01)
     call from_param("ocverb",  ocverb,.false.)
     call from_param("ocmovie", ocmovie,.false.)
     call from_param("ocsolver",ocsolver,0)
     call from_param("ocnlap",  ocnlap,2)
     call from_param("ocwght",  ocwght,.false.)
  end if
  if(wpar%nmpi >0) call MPI_BCAST(ocniter, 1,MPI_INTEGER,0,MPI_COMM_WORLD,mpierr)
  !----------------------------------------------------------------
  call weimva_init(SLin=weiop_slo1_init	&
  ,                WCin=weiop_mwc1_init	&
  ,                SCin=weiop_bor1_init	&
  ,                FKin=weiop_wem3_init	&
  ,                FXin=weiop_ssf3_init	&
  ,                IGin=weiop_hcig_init	)
  !----------------------------------------------------------------
  if(wpar%impi==0) call weipro_report()
  !----------------------------------------------------------------

  if(wpar%impi==0) then

     !! create data residual file
     name = '_res'; res_ = oc_allocate(.false.,name)
     call auxputch("n1","i",ocniter,trim(res_))
     call auxputch("esize","i",4,trim(res_))

     call from_param("lapnf",nf,3)

     if(ocsolver==1 .or. ocsolver==2) then
        call oc_allocatefile(file,wtag%dS,nint(ocmaxmem))
        do while(2*nf+1 > file%n(1))
           nf=nf-1
        end do
        allocate(nnew(file%nd)); nnew = file%n
        call oc_deallocatefile(file)
        call oc_laplacian_init(wtag%dS,nf,15,0.01,nint(ocmaxmem),nlap=ocnlap,verb=ocverb)
     end if

     call oc_sos_init(wtag%dR)
  end if

  if(wpar%impi==0 .and. wpro%verb) call weimessage("start inversion")
  !----------------------------------------------------------------

  !! select a solver
  select case(ocsolver)
  case(0) !! simple
     if(wpar%impi==0 .and. wpro%verb) call weimessage("solving Lm=D")

     call oc_tst_solver_init(wpar%impi,ocniter,nint(ocmaxmem),ocverb,ocmovie,ocmovie,wght=ocwght)
     call oc_tst_solver( L=weimva &
     ,                   S=oc_cgstep &
     ,                   W=oc_sos &
     ,                   m_=wtag%dS &
     ,                   d_=wtag%dR &
     ,                   res_=res_ &
     ,                   op1=weiop_slo1	&
     ,                   op2=weiop_mwc1	&
     ,                   op3=weiop_bor1	&
     ,                   op4=weiop_wem3	&
     ,                   op5=weiop_ssf3	&
     ,                   op6=weiop_hcig ) 

  case(1) !! regularization
     if(wpar%impi==0 .and. wpro%verb) call weimessage("solving WLm=WD & eAm=0")

     call oc_mpi_solverreg_init(wpar%impi,ocniter,oceps,nint(ocmaxmem)&
     ,                          ocverb,ocmovie,ocmovie,wght=ocwght)
     call oc_mpi_solverreg( L=weimva &
     ,                      A=oc_laplacianN & 
     ,                      S=oc_cgstep &
     ,                      W=oc_sos &
     ,                      m_=wtag%dS,d_=wtag%dR,nreg=nnew,res_=res_ &
     ,                      op1=weiop_slo1 &
     ,                      op2=weiop_mwc1 &
     ,                      op3=weiop_bor1 &
     ,                      op4=weiop_wem3 &
     ,                      op5=weiop_ssf3 &
     ,                      op6=weiop_hcig )

  case(2) !! preconditioning
     if(wpar%impi==0 .and. wpro%verb) call weimessage("solving WLP p=WD &  eIp=0")

     call oc_mpi_solverpre_init(wpar%impi,ocniter,oceps,nint(ocmaxmem) &
     ,                          ocverb,ocmovie,ocmovie,wght=ocwght)
     call oc_mpi_solverpre( L=weimva &
     ,                      P=oc_ilaplacianN &
     ,                      S=oc_cgstep &
     ,                      W=oc_sos &
     ,                      m_=wtag%dS,d_=wtag%dR,npre=nnew,res_=res_ &
     ,                      op1=weiop_slo1 &
     ,                      op2=weiop_mwc1 &
     ,                      op3=weiop_bor1 &
     ,                      op4=weiop_wem3 &
     ,                      op5=weiop_ssf3 &
     ,                      op6=weiop_hcig )

  case default
     call seperr('unknown operation')
  end select

  if(wpar%impi==0 .and. wpro%verb) call weimessage("OK")
  !----------------------------------------------------------------  
  !! shutdown MPI
  call MPI_FINALIZE(mpierr) 
  !----------------------------------------------------------------
  call exit(0)
end program WEItst_main




