 program rtm2

 use sep
 use ext_mod2
 use igrad1

 implicit none

 logical					:: diagnostic
 real, allocatable				:: source(:),data(:,:),vel(:,:),us(:,:,:)
 real, allocatable				:: dsource(:)
 real						:: xmax,zmax
 real						:: oto,dto
 integer					:: is,it1,it2,ix,iz,ixw,iwav,stat,ih,ixm,im,i3,iwrk
 integer					:: nhpad,nto,nwav2
 real						:: oxs,oxh,oxm,oz,ot,owav,ovx,ovz,owrk
 real						:: dxs,dxh,dxm,dz,dt,dwav,dvx,dvz,dwrk
 integer					:: nxs,nxh,nxm,nz,nt,nwav,nvx,nvz,nwrk

! type(axis3)					:: xs,xh,xm,z,t,wav,vx,vz
! type(axis2)					:: xwrk

! type(sep3d)					:: axgr
 

 call sep_init()

 call from_history("n1",nt); call from_history("d1",dt); call from_history("o1",ot)
 call from_history("n2",nxh); call from_history("d2",dxh); call from_history("o2",oxh)  
 call from_history("n3",nxs); call from_history("d3",dxs); call from_history("o3",oxs)

 call from_aux("source","n1",nwav); call from_aux("source","d1",dwav); call from_aux("source","o1",owav)

 call from_aux("vel","n1",nvz); call from_aux("vel","d1",dvz); call from_aux("vel","o1",ovz)
 call from_aux("vel","n2",nvx); call from_aux("vel","d2",dvx); call from_aux("vel","o2",ovx)

 call from_param("xmin", oxm, ovx)
 call from_param("xmax", xmax, ovx+(nvx-1)*dvx)
 call from_param("nhpad", nhpad, 100)
 call from_param("zmax", zmax)
 call from_param("diagnostic",diagnostic,.false.)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! To do: 01 - Different cmp, shot, receiver spacing
!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 nz=int((zmax-ovz)/dvz)+1
 oz=ovz
 dz=dvz
 nxm=int((xmax-oxm)/dvx)+1
 nwrk=2*nhpad+nxh
 dwrk=dvx
 dxm=dvx
 ixw=nhpad-int(oxh/dxh)+1

 call to_history("n1",nt); call to_history("d1",dt); call to_history("o1",ot)
 call to_history("n2",nxh); call to_history("d2",dxh); call to_history("o2",oxh)  
 call to_history("n3",nxs); call to_history("d3",dxs); call to_history("o3",oxs)  

 write(0,*) "**************************************************"
 write(0,*) "***************** Parameters *********************"
 write(0,*) " "
 write(0,*) " ns = ",nxs, " ds = ",dxs," os = ",oxs
 write(0,*) " nt = ",nt, " dt = ",dt," nh = ",nxh," dh = ",dxh," oh = ",oxh
 write(0,*) " nwav = ",nwav, " dwav = ",dwav," owav = ",owav
 write(0,*) " nvx = ",nvx," dvx = ",dvx," ovx = ",ovx," nvz = ",nz," dvz = ",dz," ovz = ",oz
 write(0,*) " xmin = ",oxm," xmax = ",xmax," nhpad = ", nhpad
 write(0,*) " nxm = ",nxm," nwrk = ",nwrk," dwrk = ",dwrk," ixw = ",ixw
 write(0,*) "**************************************************"
 write(0,*) "**************************************************"


!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  if (nxm>nvx) call seperr("Image domain is bigger than the velocity model in x - pad velocity (use extent = 1)")

  allocate(source(nwav),vel(nz,nvx))
  allocate(dsource(nwav))

  nto=nt;oto=ot;dto=dt
  nwav2=int(abs(owav)/dwav)+1
  nt=nt+nwav2
  ot=owav

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Diagnostic data
!
  allocate(us(nz,nwrk,nt))

  call sreed("source", source, 4*nwav)
  dsource=0
  stat = igrad1_lop (.true.,.false.,dsource,source)
  deallocate(source)

  call sreed("vel", vel, 4*nz*nvx)

  call extr_init_o (oxs,oxh,oxm,oz,ot,owav)
  call extr_init_d (dxs,dxh,dxm,dz,dt,dwav,dwrk)
  call extr_init_n (nxs,nxh,nxm,nz,nt,nwav,nwrk)
  call extr_init_i (ixw,nhpad)
!
! Shot loop
!
  do is=1,nxs
      us=0.
      owrk=0;ixm=0;
      owrk=(is-1)*dxs+oxs+oxh-nhpad*dwrk
      ixm=floor((owrk-oxm)/dxm)+1
      call extr_init2(owrk,ixm)

      write(0,*) "    OWRK          OIMAGE          DX               IXM" 
      write(0,*) owrk,oxm,dxm,ixm

!
! Propagates source wavefield
!

     do it1=2,nt-1

	 if (mod(it1,250).eq.0) write(0,*) "shot = ", is," from ",nxs," time = ",it1," from ",nt
      	 if(it1 <= nwav) us(1,ixw,it1)=us(1,ixw,it1)+dsource(it1-1:it1)
	 us(:,:,3)=0.
!	write(0,*) "***passei01***"
	 call sht_ext(vel,us)
!	write(0,*) "***passei02***", it1
!	 if (it1>=nwav2) then
!	     call srite("out",us(1,:,3),4*nxh)
!	write(0,*) "***passei03***"
!	 end if

	 us(:,:,1:2)=us(:,:,2:3)
     end do
     call srite("out",us(1,:,nwav2:nt-nwav2-1),4*nxh*nto)
  end do

end program
