 module ext_SO_mod

!  Extrapolacao de onda acustica em meio 2D
!  por diferencas finitas 8a ordem em z,x

 implicit none

  public :: extr_init, extr_init2, sht_ext, rec_ext, imaging

  integer,private            			:: ns,nh,nxm,nxso,nz,nt,nwav,nwrk,nhpad,ixw,ixm
  real,private        				:: ds,dh,dxm,dxso,dz,dt,dwav,dwrk
  real,private        				:: os,oh,oxm,oxso,oz,ot,owav,owrk,kz2,kx2

  type axis
     integer 					:: n
     real    					:: o, d
  end type axis

  private
  public :: axis

 contains

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 subroutine extr_init(ixw,xs,xh,xwrk,xm,xso,z,t,wav,nhpad)

 type(axis),intent(in)			:: xs,xh,xm,xso,z,t,wav
 type(axis)				:: xwrk
 integer,intent(in)			:: nhpad
 integer				:: ixw
 
 ns=xs%n;ds=xs%d;os=xs%o
 nh=xh%n;dh=xh%d;oh=xh%o
 nxm=xm%n;dxm=xm%d;oxm=xm%o
 nxso=xso%n;dxso=xso%d;oxso=xso%o
 nz=z%n;dz=z%d;oz=z%o
 nt=t%n;dt=t%d;ot=t%o
 nwav=wav%n;dwav=wav%d;owav=wav%o
 nwrk=xwrk%n;dwrk=xwrk%d;owrk=xwrk%o
 
 kx2=(t%d/xwrk%d)**2
 kz2=(t%d/z%d)**2


 write(0,*) "**************************************************"
 write(0,*) "***************** Parameters *********************"
 write(0,*) " "
 write(0,*) " ns = ",ns, " ds = ",ds," os = ",os
 write(0,*) " nt = ",nt, " dt = ",dt," nh = ",nh," dh = ",dh," oh = ",oh
 write(0,*) " nwav = ",nwav, " dwav = ",dwav," owav = ",owav
 write(0,*) " nhpad = ", nhpad," nxm = ",nxm," ixw = ",ixw
 write(0,*) " nxso = ",nxso, " dxso = ",dxso," oxso = ",oxso
 write(0,*) " nwrk = ",nwrk," dwrk = ",dwrk," owrk = ",owrk," nvz = ",nz," dvz = ",dz," ovz = ",oz
 write(0,*) "**************************************************"


 end subroutine extr_init
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

 subroutine extr_init2(owrk_in,ixm_in)

 integer,intent(in)			:: ixm_in
 real					:: owrk_in

 owrk=owrk_in;ixm=ixm_in

 end subroutine extr_init2

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

subroutine sht_ext(vel,u)

real,intent(inout)			:: u(:,:,:)
real,intent(in)				:: vel(:,:)
real, allocatable			:: velw(:,:)!,taper(:,:)
integer					:: i,k,   ni
real					:: X,T,Z,alfa


!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Taper Sponge-like: Bording, P - Finite difference modeling - nearly optimal sponge boundary conditions
!                   SEG Expanded Abstracts 23, 1921 (2004)
!
allocate(velw(nz,nwrk))
velw=0;
!allocate(taper(nz,nwrk),velw(nz,nwrk))
!taper=1.
!do i=nz-9,nz
!    taper(i,:)=exp(((-1)*(i-nz+9)*0.008)**2)
!end do
!do i=1,10
!    taper(:,i)=exp(((-1)*(10-i)*0.004)**2)
!end do
!do i=nwrk-9,nwrk
!    taper(:,i)=exp(((-1)*(i-nwrk+9)*0.004)**2)
!end do


!
! Indexing velocity
!

ni=int((owrk-oxm)/dwrk)+1
velw(:,:)=vel(:,ni:ni+nwrk-1)

!
! Loops
!

!write(0,*) "nwrk",nwrk
    do k=1,nwrk
        do i=2,nz
            Z=0;X=0;T=0
            alfa=velw(i,k)*dt
!
! 8th order extrapolator in space: Fornberg, B - Mathematics of Computation, Vol. 51, No. 184. (Oct., 1988)
!
            if (k>4.and.k<(nwrk-3)) then
                X=kx2*velw(i,k)**2*(-0.00178571428571*(u(i,k+4,2)+u(i,k-4,2))+ & 
                                0.02539682539683*(u(i,k+3,2)+u(i,k-3,2))- &
                                0.20000000000000*(u(i,k+2,2)+u(i,k-2,2))+ &
                                1.60000000000000*(u(i,k+1,2)+u(i,k-1,2))- &
                                2.84722222222222*(u(i,k,2)))
            end if

            if (i==2) then
                Z=kz2*velw(i,k)**2*(u(i+1,k,2)+u(i-1,k,2)-2*u(i,k,2))
            else if (i==3.or.i==4) then
                Z=kz2/12*velw(i,k)**2*(-u(i+2,k,2)-u(i-2,k,2)+16*u(i+1,k,2)+16*u(i-1,k,2)-30*u(i,k,2))
            else if (i==4) then
                Z=kz2/12*velw(i,k)**2*(-u(i+2,k,2)-u(i-2,k,2)+16*u(i+1,k,2)+16*u(i-1,k,2)-30*u(i,k,2))
            else if (i>4.and.i<(nz-3)) then
                Z=kz2*velw(i,k)**2*(-0.00178571428571*(u(i+4,k,2)+u(i-4,k,2))+ & 
                                0.02539682539683*(u(i+3,k,2)+u(i-3,k,2))- &
                                0.20000000000000*(u(i+2,k,2)+u(i-2,k,2))+ &
                                1.60000000000000*(u(i+1,k,2)+u(i-1,k,2))- &
                                2.84722222222222*(u(i,k,2)))
            end if

            T=2*u(i,k,2)-u(i,k,1)

            u(i,k,3)=(Z+X+T)
!
! Boundary conditions: Reynolds, A. - GEOPHYSICS. VOL. 43, NO. 6 (OCTOBER 1978)
!
            if (k<=4) then
               u(i,k,3)=u(i,k,2)+u(i,k+1,2)-u(i,k+1,1)+ &
                   alfa/dwrk*(u(i,k+1,2)-u(i,k,2)-(u(i,k+2,1)-u(i,k+1,1)))
            else if (k>=nwrk-3) then
               u(i,k,3)=u(i,k,2)+u(i,k-1,2)-u(i,k-1,1)- &
                   alfa/dwrk*(u(i,k,2)-u(i,k-1,2)-(u(i,k-1,1)-u(i,k-2,1)))
            end if
            
            if (i>=nz-3) then
               u(i,k,3)=u(i,k,2)+u(i-1,k,2)-u(i-1,k,1)- &
                   alfa/dz*(u(i,k,2)-u(i-1,k,2)-(u(i-1,k,1)-u(i-2,k,1)))   
            end if

!
! Boundary condition at z=0
!           
!            if (i<4 .and. i>4 .and. i<(nz-3)) then
!		u(i,k,3)=kx2*velw(i,k)**2*(-0.00178571428571*(u(i,k+4,2)+u(i,k-4,2))+ & 
!                                0.02539682539683*(u(i,k+3,2)+u(i,k-3,2))- &
!                                0.20000000000000*(u(i,k+2,2)+u(i,k-2,2))+ &
!                                1.60000000000000*(u(i,k+1,2)+u(i,k-1,2))- &
!                                2.84722222222222*(u(i,k,2)))		+ &
!				alfa**2/dz*(u(i+1,k,2)-u(i,k,2))	+ &
!				2*u(i,k,2) - u(i,k,1)

!!               u(i,k,3)=u(i,k,2)+u(i+1,k,2)-u(i+1,k,1)+ &
!!                   alfa/dz*(u(i+1,k,2)-u(i,k,2)-(u(i+2,k,1)-u(i+1,k,1)))   
!            end if


        end do
    end do
 deallocate(velw)
end subroutine sht_ext

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

subroutine rec_ext(it,vel,u)

real,intent(inout)			:: u(:,:,:)
real,intent(in)				:: vel(:,:)
integer,intent(in)			:: it
real, allocatable			:: velw(:,:)!,taper(:,:)	
integer					:: i,k,   ni
real					:: X,T,Z,alfa

!
! Taper Sponge-like: Bording, P - Finite difference modeling - nearly optimal sponge boundary conditions
!                   SEG Expanded Abstracts 23, 1921 (2004)
!
allocate(velw(nz,nwrk))
velw=0.
!allocate(taper(nz,nwrk),velw(nz,nwrk))
!taper=1.
!do i=nz-9,nz
!    taper(i,:)=exp(((-1)*(i-nz+9)*0.008)**2)
!end do
!do i=1,10
!    taper(:,i)=exp(((-1)*(10-i)*0.004)**2)
!end do
!do i=nwrk-9,nwrk
!    taper(:,i)=exp(((-1)*(i-nwrk+9)*0.004)**2)
!end do

!
! Indexing velocity
!

ni=int((owrk-oxm)/dwrk)+1
velw(:,:)=vel(:,ni:ni+nwrk-1)

!
! Loops
!
    do k=1,nwrk
        do i=1,nz
            Z=0;X=0;T=0
            alfa=velw(i,k)*dt
            
            if (it==nt) then
                T=u(i,k,2)
            else
                T=u(i,k,3)-2*u(i,k,2)
            end if
            
            if (k>4 .and. k<(nwrk-3)) then
                X=dwrk**(-2)*(-0.00178571428571*(u(i,k+4,2)+u(i,k-4,2))+ & 
                            0.02539682539683*(u(i,k+3,2)+u(i,k-3,2))- &
                            0.20000000000000*(u(i,k+2,2)+u(i,k-2,2))+ &
                            1.60000000000000*(u(i,k+1,2)+u(i,k-1,2))- &
                            2.84722222222222*(u(i,k,2)))
            end if
                        
            if (i==1) then
                Z=dz**(-2)*(u(i+1,k,2)-2*u(i,k,2))
            else if (i==2) then
                Z=dz**(-2)*(u(i+1,k,2)+u(i-1,k,2)-2*u(i,k,2))
            else if (i==3 .or. i==4) then
                Z=dz**(-2)/12*(-u(i+2,k,2)-u(i-2,k,2)+16*u(i+1,k,2)+16*u(i-1,k,2)-30*u(i,k,2))
            else if (i>4 .and. i<(nz-3)) then
                Z=dz**(-2)*(-0.00178571428571*(u(i+4,k,2)+u(i-4,k,2))+ & 
                            0.02539682539683*(u(i+3,k,2)+u(i-3,k,2))- &
                            0.20000000000000*(u(i+2,k,2)+u(i-2,k,2))+ &
                            1.60000000000000*(u(i+1,k,2)+u(i-1,k,2))- &
                            2.84722222222222*(u(i,k,2)))
            end if
            
            u(i,k,1)=-T+(Z+X)*(velw(i,k)*dt)**(2)
            
!
! Boundary conditions: Reynolds, A. - GEOPHYSICS. VOL. 43, NO. 6 (OCTOBER 1978)
!
            if (k<=4) then
               u(i,k,1)=u(i,k,2)+u(i,k+1,2)-u(i,k+1,3)+ &
                   alfa/dwrk*(u(i,k+1,2)-u(i,k,2)-(u(i,k+2,3)-u(i,k+1,3)))
            else if (k>=nwrk-3) then
               u(i,k,1)=u(i,k,2)+u(i,k-1,2)-u(i,k-1,3)- &
                   alfa/dwrk*(u(i,k,2)-u(i,k-1,2)-(u(i,k-1,3)-u(i,k-2,3)))
            end if
            
            if (i>=nz-3) then
               u(i,k,1)=u(i,k,2)+u(i-1,k,2)-u(i-1,k,3)- &
                   alfa/dz*(u(i,k,2)-u(i-1,k,2)-(u(i-1,k,3)-u(i-2,k,3)))   
            end if

            end do
        end do
 deallocate(velw)
end subroutine rec_ext
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine imaging(us,ur,imgx,imgz)

real,intent(out),optional		:: imgz(:,:,:)
real,intent(out)			:: imgx(:,:,:)
real,intent(in)				:: us(:,:),ur(:,:)
integer					:: i,k,j,m,ihso,ivso

imgx=0.

if (present(imgz)) then
   imgz=0.
   do j=1,nxso
      ihso=nint(((j-1)*dxso+oxso)/dwrk)
      ivso=nint(((j-1)*dxso+oxso)/dz)
      do k=1,nwrk
          do i=1,nz
	     if ((i-ivso<1).or.(i-ivso>nz).or.(i+ivso<1).or.(i+ivso>nz)) then
	       cycle
	     else
	       imgz(i,k,j)=imgz(i,k,j)+us(i-ivso,k)*ur(i+ivso,k)
	     end if
	     if ((k+ihso<1).or.(k+ihso>nwrk).or.(k-ihso<1).or.(k-ihso>nwrk)) then
	       cycle
	     else
	       imgx(i,k,j)=imgx(i,k,j)+us(i,k-ihso)*ur(i,k+ihso)
	     end if
          end do
       end do
   end do
else
   do j=1,nxso
      ihso=nint(((j-1)*dxso+oxso)/(dwrk))
      do k=1,nwrk
          do i=1,nz
	     if ((k-ihso<1).or.(k-ihso>nwrk).or.(k+ihso<1).or.(k+ihso>nwrk)) then
	        cycle
	     else 
	        imgx(i,k,j)=imgx(i,k,j)+us(i,k-ihso)*ur(i,k+ihso)
	     end if
         end do
     end do
  end do
end if
end subroutine
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

end module
