# update2 # second order time update #% subroutine update2( dt2, u , bnd ) #include "commons" real dt2 WFIELD(u) WFIELDLAY(u) BOUND(bnd) BOUNDLAY(bnd)do ic = 1, ncomp { where(bnd==INTERIOR .or. bnd==FREE) u(3,ic,:,:,:)=2.*u(2,ic,:,:,:)-u(1,ic,:,:,:)+dt2*u(3,ic,:,:,:) u(1,ic,:,:,:)=u(2,ic,:,:,:) u(2,ic,:,:,:)=u(3,ic,:,:,:) endwhere }

return end

If we want to calculate to higher orders of time we can use the general
version of the time update of order n. That is used in subroutine
`update()` , which incorporates `update2()`
as the first step in the `1..norder` loop.

# update # time update using Taylor expansion of order n # that allows multiple source with different signature # # default will be one step of this routine per timestep %subroutine update(norder,u,nts,bnd,ncomp,dt2) integer norder, nts, ncomp real u(:,:,:,:,:), dt2 integer bnd(:,:,:) integer dtn,n,ic,nfac

nfac=1 dtn = 1. do n=1,norder { #the actual order is 2,4,6,8 even nfac = nfac*2*n dtn = dtn*dt2 sourceu = sourcef(:,:,:,n) # time differentiated sourcef call wave(u,3,..,sourceu,...) # apply -L^2+F do ic=1,ncomp { u(3,ic,:,:,:) = dt2*u(3,ic,:,:,:)*2./nfac } } do ic=1,ncomp { where(bnd==0 | bnd==1 ) u(3,ic,:,:,:) = 2.*u(2,ic,:,:,:)-u(1,ic,:,:,:)+u(3,ic,:,:,:) u(1,ic,:,:,:) = u(2,ic,:,:,:) u(2,ic,:,:,:) = u(3,ic,:,:,:) endwhere }

return end

The Chebychev time update repeatedly calculates .The key is that a recursion can be used to update the time step to high accuracy and thus the time step interval becomes less of an issue. It is very convenient, that it incorporates simpler schemes as special cases, just as loop limits.

11/17/1997