!! acoustic modeling
!! variable density
!! no boundary condition
program aWFk_
  use sep
  implicit none

  logical :: verb
  integer :: nz,nx,nt,nw, nsz, nsx
  integer :: iz,ix,it,iw, isz, isx
  real    :: dz,dx,dt,dw, dsz, dsx
  real    :: oz,ox,ot,ow, osz, osx
  real    :: energy
  integer :: bndx,bndz
  real    :: idx,idz, idt2
  real    :: c1, c2, c3, c4,  c0
  real    :: c1x,c2x,c3x,c4x
  real    :: c1z,c2z,c3z,c4z
  integer :: iop,nop

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

  !! source
  real, dimension(:),   pointer :: sou, wvl
  integer :: izsou,ixsou

  !! modulus: K = r * v^2
  character(len=128)            :: kk_
  real, dimension(:,:), pointer :: kk, vp

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

  real    :: delp
  integer :: jsnap

  !----------------------------------------------------------------
  real, dimension(:,:), pointer :: tt
  real, dimension(:,:), pointer :: iroxp,irozp
  real, dimension(:,:), pointer :: iroxm,irozm
  !----------------------------------------------------------------
  real, dimension(:), pointer :: bndz0,bndz1,bndx0,bndx1
  real, dimension(:), pointer :: xtaper,ztaper
  real                        :: xtapco,ztapco
  !----------------------------------------------------------------  
  nop = 4

  kk_ ='kk'
  ro_ ='ro'
  wfld_='wfld'

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

  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(kk_,"n1",nz,1 ) 
  call from_aux(kk_,"d1",dz,1.)
  call from_aux(kk_,"o1",oz,0.)

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

  !! source location
  call from_param("nsz",nsz,1 )
  call from_param("dsz",dsz,1.)
  call from_param("osz",osz,0.)

  call from_param("nsx",nsx,1 )
  call from_param("dsx",dsx,1.)
  call from_param("osx",osx,0.)

  !! taper coefficients
  call from_param("xtapco",xtapco,0.0002)
  call from_param("ztapco",ztapco,0.0002)

  !----------------------------------------------------------------
  !! output files
  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 to_history("n3",nsx)
  call to_history("d3",dsx)
  call to_history("o3",osx)

  call to_history("n4",nsz)
  call to_history("d4",dsz)
  call to_history("o4",osz)

  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 auxputch("n4",'i',nsx,trim(wfld_))
  call auxputch("d4",'f',dsx,trim(wfld_))
  call auxputch("o4",'f',osx,trim(wfld_))

  call auxputch("n5",'i',nsz,trim(wfld_))
  call auxputch("d5",'f',dsz,trim(wfld_))
  call auxputch("o5",'f',osz,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( kk(nz,nx)); allocate( vp(nz,nx))
  allocate(sou(nt), wvl(nt))

  !----------------------------------------------------------------  
  !! material
  call sep_read( kk(1+bndz:nz-bndz,1+bndx:nx-bndx), kk_)
  do iz=1,bndz
     kk(   iz  ,:) = kk( 1+bndz,:)
     kk(nz-iz+1,:) = kk(nz-bndz,:)
  end do
  do ix=1,bndx
     kk(:,   ix  ) = kk(:, 1+bndx)
     kk(:,nx-ix+1) = kk(:,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

  vp = sqrt(kk/ro)

  idx = 1./dx
  idz = 1./dz
  idt2= 1./dt**2
  !----------------------------------------------------------------
  !! source derivative (displacement to acceleration)
  sou=0
  wvl=0.
  call sep_read(sou(1:nw))

!  do iw=1+1,nt-1
!     wvl(iw) = (-2.*sou(iw) + sou(iw+1) + sou(iw-1))
!  end do
!  wvl = wvl * idt2
  wvl = sou
  !----------------------------------------------------------------
  !! 1/ro
  allocate(iroxp(nz,nx),iroxm(nz,nx))
  allocate(irozp(nz,nx),irozm(nz,nx))
  allocate(  tt(nz,nx))

  do iz=1+1,nz-1
     irozp(iz,:) = 2./(ro(iz+1,:)+ro(iz  ,:))
     irozm(iz,:) = 2./(ro(iz  ,:)+ro(iz-1,:))
  end do
  irozp( 1,:) = irozp(   2,:)
  irozp(nz,:) = irozp(nz-1,:)
  irozm( 1,:) = irozm(   2,:)
  irozm(nz,:) = irozm(nz-1,:)

  do ix=1+1,nx-1
     iroxp(:,ix) = 2./(ro(:,ix+1)+ro(:,ix  ))
     iroxm(:,ix) = 2./(ro(:,ix  )+ro(:,ix-1))
  end do
  iroxp(:, 1) = iroxp(:,   2)
  iroxp(:,nx) = iroxp(:,nx-1)
  iroxm(:, 1) = iroxm(:,   2)
  iroxm(:,nx) = iroxm(:,nx-1)
  !----------------------------------------------------------------
  !! derivative coefficients
  c1= 1.      *1225./1024.
  c2=-1./15.  *1225./1024.
  c3= 1./125. *1225./1024.
  c4=-1./1715.*1225./1024.

  c1x = c1 * idx; c1z = c1 * idz
  c2x = c2 * idx; c2z = c2 * idz
  c3x = c3 * idx; c3z = c3 * idz
  c4x = c4 * idx; c4z = c4 * idz
  !----------------------------------------------------------------
  do isz=1,nsz;    izsou = bndz + 1+((osz+(isz-1)*dsz) - oz)/dz
     do isx=1,nsx; ixsou = bndx + 1+((osx+(isx-1)*dsx) - ox)/dx
        write(0,*) "sou ",(isz-1)*nsx+isx,"/",nsx*nsz," @",osz+(isz-1)*dsz,osx+(isx-1)*dsx

        sou = wvl / ro(izsou,ixsou)

        !! MAIN loop
        up=0;um=0;uo=0;ud=0;
        do it=1,nt

!!$           !! 1st order derivative
!!$           do iz=1+nop,nz-nop
!!$              tt(iz,:) = &
!!$              irozp(iz,:) * idz * (uo(iz+1,:) - uo(iz  ,:)) - &
!!$              irozm(iz,:) * idz * (uo(iz  ,:) - uo(iz-1,:)) 
!!$           end do
!!$           ud = tt * idz
!!$
!!$           do ix=1+nop,nx-nop
!!$              tt(:,ix) = &
!!$              iroxp(:,ix) * idx * (uo(:,ix+1) - uo(:,ix  )) - &
!!$              iroxm(:,ix) * idx * (uo(:,ix  ) - uo(:,ix-1))
!!$           end do
!!$           ud = ud + tt * idx
!!$

!!$           !! 1st order derivatives
!!$           do iz=1+nop-1,nz-nop
!!$              tt(iz,:) = (uo(iz+1,:) - uo(iz  ,:)) * idz
!!$           end do
!!$           tt = tt * irozp
!!$           do iz=1+nop,nz-nop
!!$              ud(iz,:) = (tt(iz  ,:) - tt(iz-1,:)) * idz
!!$           end do
!!$
!!$           do ix=1+nop-1,nx-nop
!!$              tt(:,ix) = (uo(:,ix+1) - uo(:,ix  )) * idx
!!$           end do
!!$           tt = tt * iroxp
!!$           do ix=1+nop,nx-nop
!!$              ud(:,ix) = &
!!$              ud(:,ix) + (tt(:,ix  ) - tt(:,ix-1)) * idx
!!$           end do

           !! 4th order derivatives (GOOD)
           do iz=nop,nz-nop ! forward derivative
              tt(iz,:) = &
              c1z * (uo(iz+1,:) - uo(iz-0,:)) + &
              c2z * (uo(iz+2,:) - uo(iz-1,:)) + &
              c3z * (uo(iz+3,:) - uo(iz-2,:)) + &
              c4z * (uo(iz+4,:) - uo(iz-3,:))
           end do
           tt = tt * iroxp
           do iz=1+nop,nz-nop ! backward derivative
              ud(iz,:) = &
              c1z * (tt(iz+0,:) - tt(iz-1,:)) + &
              c2z * (tt(iz+1,:) - tt(iz-2,:)) + &
              c3z * (tt(iz+2,:) - tt(iz-3,:)) + &
              c4z * (tt(iz+3,:) - tt(iz-4,:))
           end do

           do ix=nop,nx-nop ! forward derivative
              tt(:,ix) = &
              c1x * (uo(:,ix+1) - uo(:,ix-0)) + &
              c2x * (uo(:,ix+2) - uo(:,ix-1)) + &
              c3x * (uo(:,ix+3) - uo(:,ix-2)) + &
              c4x * (uo(:,ix+4) - uo(:,ix-3))
           end do
           tt = tt * irozp
           do ix=1+nop,nx-nop ! backward derivative
              ud(:,ix) = &
              ud(:,ix) + &
              c1x * (tt(:,ix+0) - tt(:,ix-1)) + &
              c2x * (tt(:,ix+1) - tt(:,ix-2)) + &
              c3x * (tt(:,ix+2) - tt(:,ix-3)) + &
              c4x * (tt(:,ix+3) - tt(:,ix-4))
           end do

           ud = ud * kk

           !! acceleration source
           ud(izsou,ixsou) = &
           ud(izsou,ixsou) - sou(it)

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

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

  call exit(0)
end program AWFk_








