 program rtm2

 use sep
 use ext_mod2
! use sep3d_struct_mod
 use igrad1

 implicit none

 logical					:: diagnostic
 real, allocatable				:: source(:),data(:,:),vel(:,:),us(:,:,:)
 real, allocatable				:: dsource(:),ddata(:,:),ur(:,:,:),img(:,:),image(:,:)
 real, allocatable				:: ws(:,:),wr(:,:)
 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

if (not(diagnostic)) then
 call to_history("n1",nz); call to_history("d1",dz); call to_history("o1",oz)
 call to_history("n2",nxm); call to_history("d2",dxm); call to_history("o2",oxm)  
 call to_history("n3",1); call to_history("d3",1); call to_history("o3",1)  
end if

 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),ddata(nt,nxh))
  allocate(image(nz,nxm))
  allocate(data(nt,nxh))

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

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Diagnostic data
!
if (diagnostic) then
 if(exist_file("green")) then
    call to_history("n4",nxs,"green"); call to_history("d4",dxs,"green"); call to_history("o4",oxs,"green")
    call to_history("n3",nto,"green"); call to_history("d3",dt,"green"); call to_history("o3",oto,"green")
    call to_history("n2",nwrk,"green"); call to_history("d2",dwrk,"green"); 
    call to_history("n1",nz,"green"); call to_history("d1",dz,"green"); call to_history("o1",oz,"green")
end if

 if(exist_file("wrec")) then
    call to_history("n4",nxs,"wrec"); call to_history("d4",dxs,"wrec"); call to_history("o4",oxs,"wrec")
    call to_history("n3",nto,"wrec"); call to_history("d3",dt,"wrec"); call to_history("o3",-(nto-1)*dto,"wrec")
    call to_history("n2",nwrk,"wrec"); call to_history("d2",dwrk,"wrec")
    call to_history("n1",nz,"wrec"); call to_history("d1",dz,"wrec"); call to_history("o1",oz,"wrec")
 end if 
end if
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  allocate(us(nz,nwrk,3),ur(nz,nwrk,3),img(nz,nxm))
  allocate(ws(nz,nwrk),wr(nz,nwrk))

  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.;ur=0.;img=0.

      ddata=0;data=0;it1=0;it2=0;ih=0;
      call sreed("in", data, 4*nto*nxh)
!
! 1st derivative of the shot - accounts for phase
!
      do ih=1,nxh
         stat = igrad1_lop (.true.,.false.,ddata(:,ih),data(:,ih))
      end do

      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,1:2)=us(1,ixw,1:2)+dsource(it1-1:it1)
	 us(:,:,3)=0.
	 call sht_ext(vel,us)
	 if(exist_file("green")) then
 	     if (it1>=nwav2) then
		call to_history("o2",owrk,"green") 
		call srite("green",us(:,:,3),4*nz*nwrk)
	     end if
	 end if
	 if(it1 < nt-1)	 us(:,:,1:2)=us(:,:,2:3)
     end do

!
! Back-propagates receiver and source wavefield
!
     if(exist_file("wrec")) then
	call to_history("o2",owrk,"wrec")  
     end if


if(exist_file("wrec")) then
call srite("wrec",ur(:,:,1),4*nz*nwrk)
call srite("wrec",ur(:,:,1),4*nz*nwrk)
end if
     do it2=nt-nwav2,3,-1

	if (mod(it2,250).eq.0) write(0,*) "receivers from shot = ", is," from ",nxs," time = ",it2," from ",nt

	ur(1,nhpad+1:nhpad+nxh,2:3) = ur(1,nhpad+1:nhpad+nxh,2:3) + transpose(ddata(it2-1:it2,:))

        call rec_ext(it2,vel,ur)
        if(exist_file("wrec")) then
		call srite("wrec",ur(:,:,1),4*nz*nwrk)
	end if

	if(not(diagnostic)) then
	us(:,:,1)=0.
	call rec_ext(it2,vel,us)

!
! Applies the imaging condition
!

	if(mod(it2+2,4).eq.1) then
	   ws=0.;wr=0.
	   ws=us(:,:,1);wr=ur(:,:,1)
           call imaging(ws,wr,img)
           image(:,ixm:ixm+nwrk-1)=image(:,ixm:ixm+nwrk-1)+img
	endif

	ur(:,:,2:3)=ur(:,:,1:2);ur(:,:,1)=0.
	us(:,:,2:3)=us(:,:,1:2);us(:,:,1)=0.

	else
	ur(:,:,2:3)=ur(:,:,1:2);ur(:,:,1)=0.
	end if
     end do
  end do

  if(not(diagnostic)) call srite("out",image,4*nz*nxm)

end program
