 module ext_mod2

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

 implicit none

  public :: extr_init2, sht_ext, rec_ext, imaging,extr_init_o,extr_init_d,extr_init_n,extr_init_i

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

 contains

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 subroutine extr_init_o(oxs_in,oxh_in,oxm_in,oz_in,ot_in,owav_in)
 real			::oxs_in,oxh_in,oxm_in,oz_in,ot_in,owav_in

 os=oxs_in;oh=oxh_in;oxm=oxm_in
 oz=oz_in;ot=ot_in;owav=owav_in

 end subroutine extr_init_o
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 subroutine extr_init_d(dxs_in,dxh_in,dxm_in,dz_in,dt_in,dwav_in,dwrk_in)
 real			::dxs_in,dxh_in,dxm_in,dz_in,dt_in,dwav_in,dwrk_in

 ds=dxs_in;dh=dxh_in;dxm=dxm_in
 dz=dz_in;dt=dt_in;dwav=dwav_in
 dwrk=dwrk_in

 kx2=(dt/dwrk)**2
 kz2=(dt/dz)**2

 end subroutine extr_init_d
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 subroutine extr_init_n(nxs_in,nxh_in,nxm_in,nz_in,nt_in,nwav_in,nwrk_in)
 integer			::nxs_in,nxh_in,nxm_in,nz_in,nt_in,nwav_in,nwrk_in

 ns=nxs_in;nh=nxh_in;nxm=nxm_in
 nz=nz_in;nt=nt_in;nwav=nwav_in
 nwrk=nwrk_in

 end subroutine extr_init_n
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 subroutine extr_init_i(ixw_in,nhpad_in)

 integer,intent(in)			:: ixw_in,nhpad_in

 ixw=ixw_in;nhpad=nhpad_in
 end subroutine extr_init_i
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 subroutine extr_init2(owrk_in,ixm_in)

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

 owrk=owrk_in;ixm=ixm_in
 write(0,*) owrk
 end subroutine extr_init2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine sht_ext(vel,u)

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


!
! Preliminary computations
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! To do: 01 - Different cmp, shot, receiver spacing
!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Taper Sponge-like: Bording, P - Finite difference modeling - nearly optimal sponge boundary conditions
!                   SEG Expanded Abstracts 23, 1921 (2004)
!
allocate(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=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.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)!*taper(i,k)

!
! 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)))
!               u(i,k,3)=taper(i,k)*u(i,k,3)
            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)))
!               u(i,k,3)=taper(i,k)*u(i,k,3)
            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)))   
!               u(i,k,3)=taper(i,k)*u(i,k,3)
            end if

        end do
    end do

end subroutine sht_ext

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

subroutine rec_ext(it,vel,u)

real,intent(inout)			:: u(:,:,:)
real,intent(in)				:: vel(:,:)
integer,intent(in)			:: it!,is
real, allocatable			:: velw(:,:)	
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))
!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)))
!               u(i,k,3)=taper(i,k)*u(i,k,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)))
!               u(i,k,3)=taper(i,k)*u(i,k,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)))   
!               u(i,k,3)=taper(i,k)*u(i,k,3)
            end if
            
            end do
        end do

end subroutine rec_ext
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine imaging(us,ur,img)

real					:: img(:,:)
real					:: us(:,:),ur(:,:)
integer					:: i,k

img=0.
do k=1,nwrk
    do i=1,nz
	img(i,k)=us(i,k)*ur(i,k)
     end do
end do

end subroutine
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
end module
