program rnph

 use sep
 use inicia
 use solver_irls_mod
 use cgstep_mod
 use opera
 use fileio
 use causint1
 use igrad1
 use irls

 implicit none

 character*3			:: reg
 logical				:: T=.true.,F=.false.
 real, allocatable, dimension(:)	:: m,m0
 real, allocatable, dimension(:)	:: data
 real				:: pi,eps
 integer				:: nw,nx,nz,nreal
 integer				:: niter, nmem, nfreq
 integer				:: dados,nwind

 call sep_init()
 call sep_begin_prog()

 call from_param("nw",nw,20)
 call from_param("nx",nx,50)
 call from_param("nz",nz,20)
 call from_param("nreal",nreal,1)
 call from_param("nwind",nwind,7)
 call from_param("niter",niter,100)
 call from_param("nmem",nmem,20)
 call from_param("nfreq",nfreq,10)
 call from_param("dados",dados,1)
 call from_param("eps",eps,0.1)
 call from_param("reg",reg,"nul")

 pi=4*atan(1.)

 open(unit=2,file="lixo")

  allocate(m0(nw*nx*nz*nreal),m(nw*nx*nz*nreal))
  allocate(data(nx*nx*nz*nz));data=0.

  call genmodel(m0)
  if (exist_file("m0")) then
     call io("m0",m0,nz,nx,nw,nreal)
  end if

  call gendata_init(nw,nx,nz,nwind,nreal,dados)
  call gendata(data)
  if (exist_file("data0")) then
     call io("d0",data,nz,nx)
  end if

! Nonlinear solver
  call opera_init(nw,nx,nz,nreal)
  if (reg.eq."int") then
     call solver_irls(m=m, d=data,Fop=opera_lop,Aop=igrad1_lop,stepper=cgstep,&
                      nAop=nw*nx*nz*nreal,niter=niter, eps=eps, m0=m0, nmem=nmem,&
                      nfreq=nfreq, verb=T)
  else if (reg.eq."der") then
     call solver_irls(m=m, d=data,Fop=opera_lop,Aop=igrad1_lop,stepper=cgstep,&
                      nAop=nw*nx*nz*nreal,niter=niter, eps=eps, m0=m0, nmem=nmem,&
                      nfreq=nfreq, verb=T)
  else if (reg.eq."nul") then
     call solver_irls(m=m, d=data,Fop=opera_lop,stepper=cgstep, Wdiag=l1, &
                      niter=niter, nAop=0, eps=eps, m0=m0, nmem=nmem,&
                      nfreq=nfreq, verb=T)
  end if
  deallocate(m0)
  m=m/maxval(abs(m))*pi
  if (exist_file("m")) then
     call io("mf",m,nz,nx,nw,nreal)
  end if

 if (exist_file("data")) then
    call opera_lop2(F,F,m,data,m)
    call io("df",data,nz,nx)
 end if
 close(2) 
 call sep_end_prog()
end program 
!
!  call solver_irls(m=rand2,d=data,Fop=opera_lop,stepper=cgstep,niter=niter, &
!                    Wop=,Jop=,Wdiag=,m0=m0, nmem=, nfreq=, &
!                    err=, resd=, mmov=, rmov=, verb=.true.)
