The wave operator

# apply the simple wave equation operator
# repeatedly  useable for more complicated time updates
#
subroutine wave(u,strain,stiff,stress,idim,ndim,icomp,ncomp,coef,nlo,nhi_
                                                        ,abbin,ccin,truin,bnd)
implicit none
integer                                                    nlo,nhi, bnd(:,:,:)
real             u(:,:,:,:),   strain(:,:,:,:), stress(:,:,,:)
real                                                coef(nlo:nhi),ddx(3)
#mf$layout       u(:serial,:news,:news,:news)
#mf$layout       stress(:serial,:news,:news,:news)
#mf$layout       strain(:serial,:news,:news,:news)
#mf$layout       bnd(:news,:news,:news)
integer           idim(3),ndim,icomp(3),ncomp,  abbin(3,3)

call del(0,u(2,:,:,:,:),idim,ndim,icomp,ncomp,strain,abbin,coef,nlo,nhi,ddx,bnd)

call getstress(strain,stiff,stress,idim,ndim,icomp,ncomp,abbin,ccin,truin,bnd)

where(free)	#free surface constraint
  stress(abbin(1,3)) = 0.
end where

call del(1,u(2,:,:,:,:),idim,ndim,icomp,ncomp,strain,abbin,coef,nlo,nhi,ddx,bnd)
#add in sources

return
end

Combining the previous 3 steps we can compose the wave operator $\Box$which is the left most term in (1) and is given in subroutine wave() . This routine can be used if we calculate the time update in a more complicated fashion, like Chebychev recursion.