!!
!! The routine calculates synthetic primaries from a dataset
!! previously submitted to any noise attenuation process. It
!! is based on the Aki-Richards AVO equation
!!
!!               r()=A + B[sin()]^2 + C[sin()*tan()]^2
!!
!!          
!!
!!

  program AVO

  use sep
  use mute_mod
  use cgstep_mod
  use adj_mod
  use solver_reg_mod
  use solver_prc_mod
  use solver_tiny_mod
  use avomodel
  use tcai1
  use igrad1

  implicit none
 
  character*20					:: date,tempo
  real		  		    		:: theta, time, pi
  integer	  			    	:: n, npe, npi, nsolv
  integer					:: niter
  integer                         		:: i1,i2,i3,j1,j2,nAop
  integer					:: n1,n2,n3
  real                            		:: o1,o2,o3,d1,d2,d3,eps
  real				    		:: emt(10),emx(10),imt(10),imx(10)
  real,dimension(:,:),allocatable 		:: A,B,C,mask
  real,dimension(:),allocatable   		:: m,mask1,data,dado
  real,dimension(:),pointer			:: aa

  call date_and_time(date,tempo)
  write(0,*) date,'  ',tempo

  call sep_init()
  call sep_begin_prog()
  !! get info from history file
  call from_history("n1",n1)
  call from_history("n2",n2)
  call from_history("n3",n3)
  call from_history("d1",d1)
  call from_history("d2",d2)
  call from_history("d3",d3)
  call from_history("o1",o1)
  call from_history("o2",o2)
  call from_history("o3",o3)

  !! get info from parameter
  call from_param("niter",niter,50)
  call from_param("eps",eps,1.)
  call from_param("solv",nsolv,3)

  !! write info out to history file
  call to_history("n1",n1)
  call to_history("n2",n2)
  call to_history("n3",n3)
  call to_history("d1",d1)
  call to_history("d2",d2)
  call to_history("d3",d3)
  call to_history("o1",o1)
  call to_history("o2",o2)
  call to_history("o3",o3)

!!  get mute parameters
  call from_param("npe",npe,0) 
  call from_param("emt1",emt(1))
  call from_param("emt2",emt(2))
  call from_param("emt3",emt(3))
  call from_param("emt4",emt(4))
  call from_param("emt5",emt(5))
  call from_param("emt6",emt(6))
  call from_param("emt7",emt(7))
  call from_param("emt8",emt(8))
  call from_param("emt9",emt(9))
  call from_param("emt10",emt(10))
  call from_param("emx1",emx(1))
  call from_param("emx2",emx(2))
  call from_param("emx3",emx(3))
  call from_param("emx4",emx(4))
  call from_param("emx5",emx(5))
  call from_param("emx6",emx(6))
  call from_param("emx7",emx(7))
  call from_param("emx8",emx(8))
  call from_param("emx9",emx(9))
  call from_param("emx10",emx(10))
  call from_param("npi",npi,0)
  call from_param("imt1",imt(1))
  call from_param("imt2",imt(2))
  call from_param("imt3",imt(3))
  call from_param("imt4",imt(4))
  call from_param("imt5",imt(5))
  call from_param("imt6",imt(6))
  call from_param("imt7",imt(7))
  call from_param("imt8",imt(8))
  call from_param("imt9",imt(9))
  call from_param("imt10",imt(10))
  call from_param("imx1",imx(1))
  call from_param("imx2",imx(2))
  call from_param("imx3",imx(3))
  call from_param("imx4",imx(4))
  call from_param("imx5",imx(5))
  call from_param("imx6",imx(6))
  call from_param("imx7",imx(7))
  call from_param("imx8",imx(8))
  call from_param("imx9",imx(9))
  call from_param("imx10",imx(10))

  if (exist_file("A_out")) then
        call to_history('n1',n1,"A_out")
        call to_history('d1',d1,"A_out")
        call to_history('o1',o1,"A_out")
        call to_history('n2',n3,"A_out")
        call to_history('d2',d3,"A_out")
        call to_history('o2',o3,"A_out")
  end if
  if (exist_file("B_out")) then
        call to_history('n1',n1,"B_out")
        call to_history('d1',d1,"B_out")
        call to_history('o1',o1,"B_out")
        call to_history('n2',n3,"B_out")
        call to_history('d2',d3,"B_out")
        call to_history('o2',o3,"B_out")
  end if
  if (exist_file("C_out")) then
        call to_history('n1',n1,"C_out")
        call to_history('d1',d1,"C_out")
        call to_history('o1',o1,"C_out")
        call to_history('n2',n3,"C_out")
        call to_history('d2',d3,"C_out")
        call to_history('o2',o3,"C_out")
  end if


!! allocate working space
  allocate(mask(n1,n2),mask1(n1*n2*n3),m(3*n1),data(n1*n2*n3),aa(3))
  allocate(A(n1,n3),B(n1,n3),C(n1,n3))
!!
A=0.;B=0.;C=0.
nAop=size(m)+size(aa)
pi = 4*atan(1.)
data=0.;mask=1.;m=0.;dado=0.

call sreed("in",data,n1*n2*n3*4)
if ((npe.ne.0).or.(npi.ne.0)) then
   call mute_init(n1,o1,d1,n2,o2,d2,n3,mask)
   call mute(npe,npi,emx,emt,imx,imt,mask1)
end if

if (nsolv.eq.1) then
	do i3=1,n3
	  m=0.
	  if (mod(i3,10).eq.0) write(0,*) "CMP  ", i3, " of ", n3
	  call avomodel_init(n1,n2,o2,d2,mask1((i3-1)*n1*n2+1:i3*n1*n2))
	  call solver_tiny(m=m, d=data((i3-1)*n1*n2+1:i3*n1*n2), Fop=avomodel_lop,stepper=cgstep,niter=niter)
	  call avomodel_close()
	  do i1=1,n1
	     A(i1,i3)=m(i1);B(i1,i3)=m(n1+i1);C(i1,i3)=m(2*n1+i1)
	  end do
	end do
elseif (nsolv.eq.2) then
	do i3=1,n3
	  m=0.
	  if (mod(i3,10).eq.0) write(0,*) "CMP  ", i3, " of ", n3
	  call avomodel_init(n1,n2,o2,d2,mask1((i3-1)*n1*n2+1:i3*n1*n2))
	  call solver_reg(m=m, d=data((i3-1)*n1*n2+1:i3*n1*n2), &
	      Fop=avomodel_lop,Aop=igrad1_lop,stepper=cgstep,nAop=nAop,niter=niter,eps=eps,&
	      verb=.false.)
	  call avomodel_close()
	  do i1=1,n1
	     A(i1,i3)=m(i1);B(i1,i3)=m(n1+i1);C(i1,i3)=m(2*n1+i1)
	  end do
	end do
elseif (nsolv.eq.3) then
	do i3=1,n3
	  m=0.
	  if (mod(i3,10).eq.0) write(0,*) "CMP  ", i3, " of ", n3
	  call avomodel_init(n1,n2,o2,d2,mask1((i3-1)*n1*n2+1:i3*n1*n2))
	  call solver_prc(m=m, d=data((i3-1)*n1*n2+1:i3*n1*n2), &
	      Fop=avomodel_lop,Sop=igrad1_lop,stepper=cgstep,nSop=nAop,niter=niter,eps=eps,&
	      verb=.false.)
	  call avomodel_close()
	  do i1=1,n1
	     A(i1,i3)=m(i1);B(i1,i3)=m(n1+i1);C(i1,i3)=m(2*n1+i1)
	  end do
	end do
end if

deallocate(data,mask1,m)
allocate(dado(n1*n2*n3))
  !!
  !! Creates simulated primaries
  !! 
do i3=1,n3
  do i2=1,n2
	theta=o2 + (i2-1)*d2
	do i1=1,n1
		dado((i3-1)*n1*n2+(i2-1)*n1+i1)=A(i1,i3)+B(i1,i3)*sind(theta)**2+C(i1,i3)*sind(theta)**2*tand(theta)**2
	end do
  end do
end do

call srite ("out",dado,n1*n2*n3*4)
call srite ("A_out",A,n1*n3*4)
call srite ("B_out",B,n1*n3*4)
call srite ("C_out",C,n1*n3*4)
call cgstep_close()

call date_and_time(date,tempo)
write(0,*) date,'  ',tempo

deallocate(dado,A,B,C)
call sep_end_prog()
end program AVO
