 module sht_ext_mod

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

 implicit none

  public :: extr_init, sht_ext

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

  type axis
     integer :: n
     real    :: o, d, v !v=for the value of the variable (typical:v=(i-1)*d+o)
     real, dimension(:), allocatable:: all
  end type axis

  private
  public :: axis

 contains

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

 type(axis),intent(in)			:: xs,xh,xm,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
 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,*) " nwrk = ",nwrk," dwrk = ",dwrk," owrk = ",owrk," nvz = ",nz," dvz = ",dz," ovz = ",oz
 write(0,*) " nhpad = ", nhpad," nxm = "," ixw = ",ixw
 write(0,*) "**************************************************"



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

subroutine sht_ext(is,it,vel,u)

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


!
! Preliminary computations
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! To do: 01 - Different cmp, shot, receiver spacing
!!!!!!!!!!!!!!!!!!!!!!!!!!!!
owrk=(is-1)*ds+os+oh-nhpad*dxm
!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Taper Sponge-like: Bording, P - Finite difference modeling - nearly optimal sponge boundary conditions
!                   SEG Expanded Abstracts 23, 1921 (2004)
!
allocate(taper(nz,nwrk),velw(nz,nwrk))		!,u1(nz,nx,3))
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)

write(2,*) "#############passei01###############"
!
! Loops
!
!do j=2,nt-1

    do k=1,nwrk
write(2,*) "#############passei02###############",k
        do i=2,nz
write(2,*) "#############passei03###############",i
            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,it)+u(i,k-4,it))+ & 
                                0.02539682539683*(u(i,k+3,it)+u(i,k-3,it))- &
                                0.20000000000000*(u(i,k+2,it)+u(i,k-2,it))+ &
                                1.60000000000000*(u(i,k+1,it)+u(i,k-1,it))- &
                                2.84722222222222*(u(i,k,it)))
            end if

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

!            if (it==1) then
!	       T=2*u(i,k,it)
!	    else
               T=2*u(i,k,it)-u(i,k,it-1)
!	    end if

            u(i,k,it+1)=(Z+X+T)!*taper(i,k)

!
! Boundary conditions: Reynolds, A. - GEOPHYSICS. VOL. 43, NO. 6 (OCTOBER 1978)
!
            if (k<=4) then
               u(i,k,it+1)=u(i,k,it)+u(i,k+1,it)-u(i,k+1,it-1)+ &
                   alfa/dwrk*(u(i,k+1,it)-u(i,k,it)-(u(i,k+2,it-1)-u(i,k+1,it-1)))
               u(i,k,it+1)=taper(i,k)*u(i,k,it+1)
            else if (k>=nwrk-3) then
               u(i,k,it+1)=u(i,k,it)+u(i,k-1,it)-u(i,k-1,it-1)- &
                   alfa/dwrk*(u(i,k,it)-u(i,k-1,it)-(u(i,k-1,it-1)-u(i,k-2,it-1)))
               u(i,k,it+1)=taper(i,k)*u(i,k,it+1)
            end if
            
            if (i>=nz-3) then
               u(i,k,it+1)=u(i,k,it)+u(i-1,k,it)-u(i-1,k,it-1)- &
                   alfa/dz*(u(i,k,it)-u(i-1,k,it)-(u(i-1,k,it-1)-u(i-2,k,it-1)))   
               u(i,k,it+1)=taper(i,k)*u(i,k,it+1)
            end if
write(2,*) "#############passei04###############"
        end do
write(2,*) "#############passei05###############"
    end do
write(2,*) "#############passei06###############"
!end do
end subroutine
end module
