program main
  use sep
  use off2ang2d1_sd_mod
!  use off2ang2d2_sd_mod
  use solver_tiny_mod
  use solver_prc_mod
  use cgstep_mod
  use igrad1
  use ident_mod

  implicit none

  real, dimension(:), pointer :: modl(:), data(:), m(:), d(:),mask(:)
  integer                     :: nz, nhx, nax, nmx, imx, iax, ihx, iz, stat, style, mode
  integer		      :: niter,NAop
  real                        :: z0, hx0, mx0 
  real                        :: dz, dhx, dax, dmx
  real                        :: eps, axmin, axmax, hxmin, hxmax
  real, parameter             :: PI = 3.14159265359
  logical, parameter          :: T = .true. , F = .false.
  logical	              :: adj	
  call initpar()

  call from_param("style", style, 0)
  call from_param("adj", adj, T)
  call from_param("niter", niter, 50)
  call from_param("eps", eps, 1.)
  call from_param("mode",mode,1)

  if(adj) then
     call from_history("n1", nz ); call from_history("d1", dz ); call from_history("o1", z0 )
     call from_history("n2", nhx); call from_history("d2", dhx); call from_history("o2", hx0)  
     call from_history("n3", nmx); call from_history("d3", dmx); call from_history("o3", mx0)
     call from_param("amin", axmin,-60.)
     call from_param("amax", axmax, 60.)
     call from_param("dax", dax, 1.)
     nax = (axmax - axmin) / dax + 1
     call to_history("n1", nz) ; call to_history("d1", dz) ; call to_history("o1", z0)
     call to_history("n2", nax); call to_history("d2", dax); call to_history("o2", axmin) 
     call to_history("n3", nmx); call to_history("d3", dmx); call to_history("o3", mx0)

     allocate(mask(nz*nax), modl(nz*nax), data(nz*nhx), m(nz*nax), d(nz*nhx))

     mask=1
     call ident_init(nz, nax, mask) 
     nAop=size(modl)+3
     do imx=1, nmx
        data = 0;modl = 0.
        call sreed("in", data, 4*nz*nhx)
	if (style.eq.0) then
	  if (mode .eq. 0) then
             call off2ang2d1_sd_init(z0, hx0, dz, dhx, nz, nhx, nax, dax, axmin, 0)
             stat = off2ang2d1_sd_oper(adj, F, modl, data)
	  else if (mode .eq. 1) then
             call off2ang2d1_sd_init(z0, hx0, dz, dhx, nz, nhx, nax, dax, axmin, 1)
             call solver_tiny(m=modl, d=data, Fop=off2ang2d1_sd_oper, stepper=cgstep, niter=niter)
	  else
             call off2ang2d1_sd_init(z0, hx0, dz, dhx, nz, nhx, nax, dax, axmin, 1)
	     call solver_prc(m=modl, d=data,Fop=off2ang2d1_sd_oper,Sop=ident_lop,stepper=cgstep,nSop=nAop,niter=niter,eps=eps,&
	         verb=.true.)
	  end if
	else
           call off2ang2d1_sd_init(z0, hx0, dz, dhx, nz, nhx, nax, dax, axmin, 1)
           stat = off2ang2d1_sd_oper(adj, F, modl, data)
	end if	
        call srite("out", modl, 4*nz*nax)
     end do
  else
     call from_history("n1", nz) ; call from_history("d1", dz) ; call from_history("o1", z0)
     call from_history("n2", nax); call from_history("d2", dax); call from_history("o2", axmin) 
     call from_history("n3", nmx); call from_history("d3", dmx); call from_history("o3", mx0)
     call from_param("hmin", hxmin,-600.)
     call from_param("hmax", hxmax, 600.)
     call from_param("dhx", dhx, 25.)
     nhx = (hxmax - hxmin) / dhx + 1
     call to_history("n1", nz ); call to_history("d1", dz ); call to_history("o1", z0 )
     call to_history("n2", nhx); call to_history("d2", dhx); call to_history("o2", hxmin)  
     call to_history("n3", nmx); call to_history("d3", dmx); call to_history("o3", mx0)

     allocate(modl(nz*nax), data(nz*nhx), m(nz*nax), d(nz*nhx))

     do imx=1,nmx
        data = 0;modl = 0.
        call sreed("in", modl, 4*nz*nax)
!--------------------------------------
!	m=modl;modl=0.
!        do iax=1,nax
!	   do iz=1,nz
!	      modl((iz-1)*nax+iax)=m((iax-1)*nz+iz)
!	   end do
!	end do
!-----------------------------------------
	if (style.eq.0) then
           call off2ang2d1_sd_init(z0, hxmin, dz, dhx, nz, nhx, nax, dax, axmin, 0)
           stat = off2ang2d1_sd_oper(adj, F, modl, data)
	else
           call off2ang2d1_sd_init(z0, hxmin, dz, dhx, nz, nhx, nax, dax, axmin, 1)
           stat = off2ang2d1_sd_oper(adj, F, modl, data)
	end if

!-------------------------------------
!	d=data;data=0.
!       do ihx=1,nhx
!	   do iz=1,nz
!	      data((ihx-1)*nz+iz)=d((iz-1)*nhx+ihx)
!	   end do
!	end do
!---------------------------------------
        call srite("out", data, 4*nz*nhx)
     end do
  end if 

end program