 program rtm_SO_pw

 use sep
 use ext_SO_mod
 use sep3d_struct_mod
 use igrad1

 implicit none

 logical					:: forward,vodcig
 real, allocatable				:: source(:),data(:,:),vel(:,:),velw(:,:),us(:,:,:),ws(:,:,:)
 real, allocatable				:: dsource(:),ddata(:,:),ur(:,:,:),dsource2d(:,:)
 real, allocatable				:: imgx(:,:,:),imagex(:,:,:),imgz(:,:,:),imagez(:,:,:)
 real						:: xmax,zmax,oto
 integer					:: is,it,ix,iz,ixw,iwav,stat,ih,ixm,i
 integer					:: nhpad,nto,nwav2

 type(axis)					:: xs,xh,xm,xwrk,xso,z,t,wav,vx,vz
 type(sep3d)					:: axgr
 

 call sep_init()
 call sep_begin_prog()

 call from_history("n1",t%n); call from_history("d1",t%d); call from_history("o1",t%o)
 call from_history("n2",xh%n); call from_history("d2",xh%d); call from_history("o2",xh%o)  
 call from_history("n3",xs%n); call from_history("d3",xs%d); call from_history("o3",xs%o)

 call from_aux("source","n1",wav%n); call from_aux("source","d1",wav%d); call from_aux("source","o1",wav%o)

 call from_aux("vel","n1",vz%n); call from_aux("vel","d1",vz%d); call from_aux("vel","o1",vz%o)
 call from_aux("vel","n2",vx%n); call from_aux("vel","d2",vx%d); call from_aux("vel","o2",vx%o)

 call from_param("xmin", xm%o, vx%o)
 call from_param("xmax", xmax, vx%o+(vx%n-1)*vx%d)
 call from_param("nhpad", nhpad, 100)
 call from_param("nsof",xso%n,1); call from_param("dsof",xso%d,0.); call from_param("osof",xso%o,0.);
 call from_param("zmax", zmax)
 call from_param("forward",forward,.false.)
 call from_param("vodcig",vodcig,.false.)


 z%n=int((zmax-vz%o)/vz%d)+1
 z%o=vz%o
 z%d=vz%d
 xm%n=int((xmax-xm%o)/vx%d)+1
 xwrk%n=2*nhpad+xh%n
 xwrk%d=vx%d
 xm%d=vx%d
 ixw=nhpad-int(xh%o/xh%d)+1

 call to_history("n1",z%n); call to_history("d1",z%d); call to_history("o1",z%o)
 call to_history("n2",xm%n); call to_history("d2",xm%d); call to_history("o2",xm%o)  
 call to_history("n3",xso%n); call to_history("d3",xso%d); call to_history("o3",xso%o)  
 call to_history("n4",1); call to_history("d4",1); call to_history("o4",1)

 if(exist_file("wrec")) then
    call to_history("n3",t%n,"wrec"); call to_history("d3",t%d,"wrec"); call to_history("o3",t%o,"wrec")
    call to_history("n2",xwrk%n,"wrec"); call to_history("d2",xwrk%d,"wrec"); call to_history("o2",xwrk%o,"wrec")  
    call to_history("n1",z%n,"wrec"); call to_history("d1",z%d,"wrec"); call to_history("o1",z%o,"wrec")
 end if

 if(exist_file("vertcig")) then
    call to_history("n1",z%n,"vertcig"); call to_history("d1",z%d,"vertcig"); call to_history("o1",z%o,"vertcig")
    call to_history("n2",xm%n,"vertcig"); call to_history("d2",xm%d,"vertcig"); call to_history("o2",xm%o,"vertcig")  
    call to_history("n3",xso%n,"vertcig"); call to_history("d3",xso%d,"vertcig"); call to_history("o3",xso%o,"vertcig")  
    call to_history("n4",1); call to_history("d4",1); call to_history("o4",1)
 end if


 write(0,*) "**************************************************"
 write(0,*) "***************** Parameters *********************"
 write(0,*) " "
 write(0,*) " ns = ",xs%n, " ds = ",xs%d," os = ",xs%o
 write(0,*) " nt = ",t%n, " dt = ",t%d," nh = ",xh%n," dh = ",xh%d," oh = ",xh%o
 write(0,*) " nwav = ",wav%n, " dwav = ",wav%d," owav = ",wav%o
 write(0,*) " nvx = ",vx%n," dvx = ",vx%d," ovx = ",vx%o," nvz = ",z%n," dvz = ",z%d," ovz = ",z%o
 write(0,*) " xmin = ",xm%o," xmax = ",xmax," nhpad = ", nhpad
 write(0,*) " nxm = ",xm%n," nwrk = ",xwrk%n," dwrk = ",xwrk%d," ixw = ",ixw
 write(0,*) "**************************************************"
 write(0,*) "**************************************************"


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

  if (xm%n>vx%n) call seperr("Image domain is bigger than the velocity model in x - pad velocity (use extent = 1)")

  allocate(source(wav%n),vel(vz%n,vx%n),velw(z%n,vx%n),data(t%n,xh%n))
  allocate(dsource(wav%n),ddata(t%n,xh%n))
  allocate(imagex(z%n,xm%n,xso%n),imgx(z%n,xm%n,xso%n))

  if(exist_file("vertcig")) then
     allocate(imagez(z%n,xm%n,xso%n),imgz(z%n,xm%n,xso%n))
  end if

  allocate(us(z%n,xwrk%n,3),ur(z%n,xwrk%n,3))
  allocate(ws(z%n,xwrk%n,1))

  nto=t%n;oto=t%o
  nwav2=int(abs(wav%o)/wav%d)+1
  t%n=t%n+nwav2
  t%o=wav%o

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

  call sreed("vel", vel, 4*vz%n*vx%n)
  velw=vel(1:z%n,:)
  deallocate(vel)

  call extr_init(ixw,xs,xh,xwrk,xm,xso,z,t,wav,nhpad)
!
!
!
 allocate(dsource2d(xwrk%n,wav%n))
    do i=1,xwrk%n
       dsource2d(i,:)=dsource
    end do
!
! Shot loop
!
  do is=1,xs%n
      us=0.;ddata=0.;imgx=0.;ws=0.
      if (vodcig) imgz=0.
      data=0.
      call sreed("in", data, 4*nto*xh%n)
!
! 1st derivative of the shot - accounts for phase
!
      do ih=1,xh%n
         stat = igrad1_lop (.true.,.false.,ddata(:,ih),data(:,ih))
      end do

      xwrk%o=(is-1)*xs%d+xs%o+xh%o-nhpad*xwrk%d
      ixm=floor((xwrk%o-xm%o)/xm%d)+1
      call extr_init2(xwrk%o,ixm)

!
! Propagates source wavefield
!

     call init_sep3d(axgr, "OUTPUT","FLOAT", "REGULAR", ctag="green")
     call axis_allocate(axgr,3)
     axgr%o(1)=z%o; axgr%n(1)=z%n; axgr%d(1)=z%d
     axgr%o(2)=xwrk%o; axgr%n(2)=xwrk%n; axgr%d(2)=xwrk%d
     axgr%o(3)=oto; axgr%n(3)=nto; axgr%d(3)=t%d
     call sep3d_set_sep3d(axgr)
     call sep3d_write_description("green",axgr)

     do it=2,t%n-1
	 if (mod(it,250).eq.0) write(0,*) "shot = ", is," from ",xs%n," time = ",it," from ",t%n
	 if(it <= wav%n) us(1,:,1:2)=us(1,:,1:2)+dsource2d(:,it-1:it)
	 call sht_ext(velw,us)
	 if (it>=nwav2) call srite("green",us(:,:,3),4*z%n*xwrk%n)
	 us(:,:,1:2)=us(:,:,2:3);us(:,:,3)=0
     end do

     if(forward) stop

     call sep3df_delete(axgr,"green")
     call init_sep3d("green",axgr,"SCRATCH")

!
! Back-propagates receivers
!
!      ur=0

!      do it=t%n-nwav2,3,-1
!	 if (mod(it,250).eq.0) write(0,*) "receivers from shot = ", is," from ",xs%n," time = ",it," from ",t%n
!	 ur(1,nhpad+1:nhpad+xh%n,2:3) = ur(1,nhpad+1:nhpad+xh%n,2:3) + transpose(ddata(it-1:it,:))
!        call rec_ext(it,velw,ur)
!        if(exist_file("wrec")) then
!		if (it>=nwav2) then
!		   call srite("wrec",ur(:,:,1),4*z%n*xwrk%n)
!		end if
!	end if

!
! Applies the imaging condition
!

!	if(mod(it+2,4).eq.1) then
!           call sreed_window("green",3,(/z%n,xwrk%n,nto/),(/z%n,xwrk%n,1/),(/0,0,it-1/),(/1,1,1/),4,ws)
!	   if (vodcig) then
!              call imaging(ws(:,:,1),ur(:,:,1),imgx,imgz)
!              imagex(:,ixm:ixm+xwrk%n-1,:)=imagex(:,ixm:ixm+xwrk%n-1,:)+imgx
!              imagez(:,ixm:ixm+xwrk%n-1,:)=imagez(:,ixm:ixm+xwrk%n-1,:)+imgz
!	   else
!              call imaging(ws(:,:,1),ur(:,:,1),imgx)
!              imagex(:,ixm:ixm+xwrk%n-1,:)=imagex(:,ixm:ixm+xwrk%n-1,:)+imgx
!	   end if
!	end if
!   	ur(:,:,2:3)=ur(:,:,1:2);ur(:,:,1)=0

!      end do

      call sep3df_delete(axgr,"green")
!      if(exist_file("wrec")) then
!         do it=1,nwav2
!             call srite("wrec",ur(:,:,1),4*z%n*xwrk%n)
!         end do
!      end if

  end do

!  call srite("out",imagex,4*z%n*xm%n*xso%n)

!  if(exist_file("vertcig")) then 
!     call srite("vertcig",imagez,4*z%n*xm%n*xso%n)
!     deallocate(imagez,imgz)
!  end if

  deallocate(velw,data)
  deallocate(dsource,ddata)
  deallocate(imagex,imgx)
  deallocate(us,ur)
  deallocate(ws)

 call sep_end_prog()
end program
