 module ext_SO_pw_mod

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

 implicit none

  public :: extr_init, sht_ext, rec_ext, imaging

  integer,private            			:: nxm,nz,nt,nwav
  real,private        				:: dxm,dz,dt,dwav
  real,private        				:: oxm,oz,ot,owav,kz2,kx2
  real,private					:: velw(:,:)

  type axis
     integer 					:: n
     real    					:: o, d
!     real, allocatable				:: all(:)
  end type axis

  public :: axis

 contains

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 subroutine extr_init(xm,z,t,wav,vel)

 type(axis),intent(in)			:: xm,z,t,wav
 real, allocatable			:: velw(:,:)

 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
 
 kx2=(dt/dxm)**2
 kz2=(dt/dz)**2
 
 allocate(velw(nz,nxm))
 velw=vel

 end subroutine extr_init
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine sht_ext(u)

real,intent(inout)			:: u(:,:,:)
integer					:: i,k
real					:: X,T,Z,alfa

!
! Loops
!
    do k=1,nxm
        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<(nxm-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/dxm*(u(i,k+1,2)-u(i,k,2)-(u(i,k+2,1)-u(i,k+1,1)))
            else if (k>=nxm-3) then
               u(i,k,3)=u(i,k,2)+u(i,k-1,2)-u(i,k-1,1)- &
                   alfa/dxm*(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

       end do
    end do
 end subroutine sht_ext
end module
