!! exploding reflector
!! acoustic modeling
!! constant density
!! optional boundary conditions
program ERFmod_
  use sep
  implicit none

  logical :: verb   !! verbose
  logical :: abc    !! absorbing boundary condition
  integer :: nz,nx,nt,nw
  integer :: iz,ix,it,iw
  real    :: dz,dx,dt,dw
  real    :: oz,ox,ot,ow
  real    :: energy
  integer :: bndx,bndz
  real    :: idx,idz,idt2
  real    :: c0,c1,c2
  integer :: iop,nop

  !! wavefield
  character(len=128)            :: wfld_
  real, dimension(:,:), pointer :: up,um,uo,ud

  !! source
  real, dimension(:),   pointer :: sou, wvl 

  !! velocity
  character(len=128)            :: vp_
  real, dimension(:,:), pointer :: vp, vp2

  !! density
  character(len=128)            :: ro_
  real, dimension(:,:), pointer :: ro

  !! reflectors
  character(len=128)            :: rr_
  real, dimension(:,:), pointer :: rr

  real    :: delp
  integer :: jsnap
  !----------------------------------------------------------------
  real, dimension(:), pointer :: bndz0,bndz1,bndx0,bndx1
  real, dimension(:), pointer :: xtaper,ztaper
  real                        :: xtapco,ztapco
  !----------------------------------------------------------------  
  nop = 2

  vp_ ='vp'
  ro_ ='ro'
  wfld_='wfld'
  rr_ ='rr'

  call sep_init()
  call from_param("verb",verb,.false.)
  call from_param("jsnap",jsnap,20)
  call from_param("abc",abc,.false.)

  call from_param("bndz",bndz,nop); if(bndz<nop) bndz=nop
  call from_param("bndx",bndx,nop); if(bndx<nop) bndx=nop

  call from_history("n1",nw)
  call from_param("nt",nt,nw)
  call from_param("dt",dt,.0001)

  call from_aux(vp_,"n1",nz,1 )
  call from_aux(vp_,"d1",dz,1.)
  call from_aux(vp_,"o1",oz,0.)

  call from_aux(vp_,"n2",nx,1 )
  call from_aux(vp_,"d2",dx,1.)
  call from_aux(vp_,"o2",ox,0.)

  !! taper coefficients
  call from_param("xtapco",xtapco,0.0002)
  call from_param("ztapco",ztapco,0.0002)
  !----------------------------------------------------------------
  !! output files
	call to_history("esize",4)

  call to_history("n1",nx)
  call to_history("d1",dx)
  call to_history("o1",ox)

  call to_history("n2",nt)
  call to_history("d2",dt)
  call to_history("o2",0.)

  call auxputch("n1",'i',nz,trim(wfld_))
  call auxputch("d1",'f',dz,trim(wfld_))
  call auxputch("o1",'f',oz,trim(wfld_))

  call auxputch("n2",'i',nx,trim(wfld_))
  call auxputch("d2",'f',dx,trim(wfld_))
  call auxputch("o2",'f',ox,trim(wfld_))

  call auxputch("n3",'i',nt/jsnap,trim(wfld_))
  call auxputch("d3",'f',dt*jsnap,trim(wfld_))
  call auxputch("o3",'f',      0.,trim(wfld_))

  call sep_close()
  !----------------------------------------------------------------
  !! expand domain
  nz = nz + 2*bndz
  nx = nx + 2*bndx

  allocate(up(nz,nx),um(nz,nx),uo(nz,nx),ud(nz,nx))
  allocate( ro(nz,nx))
  allocate( vp(nz,nx), vp2(nz,nx))
  allocate(sou(nt), wvl(nt))
  allocate( rr(nz,nx))
  !----------------------------------------------------------------
  !! material
  call sep_read( vp(1+bndz:nz-bndz,1+bndx:nx-bndx), vp_)
  do iz=1,bndz
     vp(   iz  ,:) = vp( 1+bndz,:)
     vp(nz-iz+1,:) = vp(nz-bndz,:)
  end do
  do ix=1,bndx
     vp(:,   ix  ) = vp(:, 1+bndx)
     vp(:,nx-ix+1) = vp(:,nx-bndx)
  end do

  call sep_read( ro(1+bndz:nz-bndz,1+bndx:nx-bndx), ro_)
  do iz=1,bndz
     ro(   iz  ,:) = ro( 1+bndz,:)
     ro(nz-iz+1,:) = ro(nz-bndz,:)
  end do
  do ix=1,bndx
     ro(:,   ix  ) = ro(:, 1+bndx)
     ro(:,nx-ix+1) = ro(:,nx-bndx)
  end do

  vp2 = vp*vp

  idx = 1./dx**2
  idz = 1./dz**2
  idt2= 1./dt**2
  !----------------------------------------------------------------
  !! source derivative (displacement to acceleration)
  sou=0.
  wvl=0.
  call sep_read(sou(1:nw))
  do iw=1+1,nw-1
     wvl(iw) = (-2.*sou(iw) + sou(iw+1) + sou(iw-1))
  end do
  wvl = wvl *idt2
  !----------------------------------------------------------------
  !! exploding reflectors
  call sep_read( rr(1+bndz:nz-bndz,1+bndx:nx-bndx), rr_)
  do iz=1,bndz
     rr(   iz  ,:) = 0.
     rr(nz-iz+1,:) = 0.
  end do
  do ix=1,bndx
     rr(:,   ix  ) = 0.
     rr(:,nx-ix+1) = 0.
  end do
  rr = rr / ro
  !----------------------------------------------------------------  
  if(abc) then
     !! one-way absorbing BC setup
     allocate(bndz0(nx),bndz1(nx))
     do ix=1,nx
        delp = vp(1+nop,ix) *dt/dz
        bndz0(ix) = (1.-delp)/(1.+delp)
        delp = vp(nz-nop,ix) *dt/dz
        bndz1(ix) = (1.-delp)/(1.+delp)
     end do

     allocate(bndx0(nz),bndx1(nz))
     do iz=1,nz
        delp = vp(iz,1+nop) *dt/dx
        bndx0(iz) = (1.-delp)/(1.+delp)
        delp = vp(iz,nx-nop) *dt/dx
        bndx1(iz) = (1.-delp)/(1.+delp)
     end do

     !! taper setup
     allocate(ztaper(bndz))
     do iz=1,bndz
        ztaper(iz) = exp(-ztapco * iz*iz)
     end do

     allocate(xtaper(bndx))
     do ix=1,bndx
        xtaper(ix) = exp(-xtapco * ix*ix)
     end do
  end if
  !----------------------------------------------------------------
  !! laplacian coefficients
  c0 = -30./12.
  c1 = +16./12.
  c2 = - 1./12.
  !----------------------------------------------------------------
  !! MAIN loop
  up=0;um=0;uo=0;ud=0;
  do it=1,nt
     !! 4th order laplacian
     do iz=1+nop,nz-nop
        do ix=1+nop,nx-nop
           ud(iz,ix) = &
           c0* uo(iz,  ix  ) * (idx + idz)        + &
           c1*(uo(iz  ,ix-1) + uo(iz  ,ix+1))*idx + &
           c2*(uo(iz  ,ix-2) + uo(iz  ,ix+2))*idx + &
           c1*(uo(iz-1,ix  ) + uo(iz+1,ix  ))*idz + &
           c2*(uo(iz-2,ix  ) + uo(iz+2,ix  ))*idz
        end do
     end do

     ud= ud *vp2

     !! exploding reflectors
     ud = ud - sou(it) * rr

     if(abc) then
        !! taper
        do iz=1,nz
           uo(iz,        1:bndx) = uo(iz,        1:bndx) * xtaper(bndx:1:-1) 
           uo(iz,nx-bndx+1:  nx) = uo(iz,nx-bndx+1:  nx) * xtaper(1:bndx:+1)
           um(iz,        1:bndx) = um(iz,        1:bndx) * xtaper(bndx:1:-1) 
           um(iz,nx-bndx+1:  nx) = um(iz,nx-bndx+1:  nx) * xtaper(1:bndx:+1)
           ud(iz,        1:bndx) = ud(iz,        1:bndx) * xtaper(bndx:1:-1) 
           ud(iz,nx-bndx+1:  nx) = ud(iz,nx-bndx+1:  nx) * xtaper(1:bndx:+1)
        end do
        do ix=1,nx
           uo(        1:bndz,ix) = uo(        1:bndz,ix) * ztaper(bndz:1:-1) 
           uo(nz-bndz+1:  nz,ix) = uo(nz-bndz+1:  nz,ix) * ztaper(1:bndz:+1)
           um(        1:bndz,ix) = um(        1:bndz,ix) * ztaper(bndz:1:-1) 
           um(nz-bndz+1:  nz,ix) = um(nz-bndz+1:  nz,ix) * ztaper(1:bndz:+1)
           ud(        1:bndz,ix) = ud(        1:bndz,ix) * ztaper(bndz:1:-1) 
           ud(nz-bndz+1:  nz,ix) = ud(nz-bndz+1:  nz,ix) * ztaper(1:bndz:+1)
        end do
     end if

     !! time step
     up = 2*uo - um + ud * dt**2
     um =   uo
     uo =   up

     if(abc) then
        !! absorbing BC
        do iop=1,nop
           uo( 1+nop-iop  ,:) = &
           um( 1+nop-iop+1,:) + bndz0 * (um( 1+nop-iop,:)-uo( 1+nop-iop+1,:))
           uo(nz-nop+iop  ,:) = &
           um(nz-nop+iop-1,:) + bndz1 * (um(nz-nop+iop,:)-uo(nz-nop+iop-1,:))
           uo(:, 1+nop-iop  ) = &
           um(:, 1+nop-iop+1) + bndx0 * (um(:, 1+nop-iop)-uo(:, 1+nop-iop+1))
           uo(:,nx-nop+iop  ) = &
           um(:,nx-nop+iop-1) + bndx1 * (um(:,nx-nop+iop)-uo(:,nx-nop+iop-1))
        end do
     end if

     if(mod(it,jsnap)==0) then
        call sep_write(uo(bndz+1:nz-bndz,bndx+1:nx-bndx),wfld_)
     end if
     call sep_write(uo(1+bndz,1+bndx:nx-bndx)) 
     if(verb) write(0,*) it,sum(uo)
  end do

  call exit(0)
end program ERFmod_








