subroutine copyq1(maxmx,meqn,mbc,mx,q1,q2)

copy the contents of q1 into q2 c

implicit double precision (a-h,o-z)
dimension q1(1-mbc:maxmx+mbc, meqn)
dimension q2(1-mbc:maxmx+mbc, meqn)

do 10 i = 1-mbc, mx+mbc

10 continue



    subroutine claw1(maxmx,meqn,mwaves,mbc,mx, q,aux,dx,tstart,tend,dtv,cflv,nv,method,mthlim, work,mwork,info,bc1,rp1,src1)


Solves a hyperbolic system of conservation laws in one space dimension of the general form
capa * q_t + A q_x = psi
The "capacity function" capa(x) and source term psi are optional (see below). For a more complete description see the documentation in the directory claw/doc, especially and Sample driver programs and user-supplied subroutines are available  see the the directory claw/clawpack/1d/example for one example, and codes in claw/applications for more extensive examples. The user must supply the following subroutines: cbc1, rp1 subroutines specifying the boundary conditions and Riemann solver. These are described in greater detail below. In addition, if the equation contains source terms psi, then the user  must provide: src1 subroutine that solves capa * q_t = psi   over a single time step. These routines must be declared EXTERNAL in the main program. For description of the calling sequences, see below.

The user must supply the following subroutines: bc1, rp1 subroutines specifying the boundary conditions and Riemann solver. These are described in greater detail below. In addition, if the equation contains source terms psi, then the user must provide: src1 subroutine that solves capa * q_t = psi over a single time step. These routines must be declared EXTERNAL in the main program. For description of the calling sequences, see below.

Description of parameters

maxmx : maximum number of interior grid points in x, and is used in declarations of the array q.

meqn    : number of equations in the system of conservation laws.

mwaves : number of waves that result from the solution of each Riemann problem. Often mwaves = meqn but for some problems these may be different.

mbc      :  number of "ghost cells" that must be added on to each side of the domain to handle boundary conditions. The cells actually in the physical domain are labelled from 1 to mx in x. The arrays are dimensioned actually indexed from 1-mbc to mx+mbc. For the methods currently implemented, mbc = 2 should be used. If the user implements another method that has a larger stencil and hence requires more ghost cells, a larger value of mbc could be used. q is extended from the physical domain to the ghost cells by the user-supplied routine bc1.

mx     :   number of grid cells in the x-direction, in the physical domain. In addition there are mbc grid cells along each edge of the grid that are used for boundary conditions. Must have mx .le. maxmx q(1-mbc:maxmx+mbc, meqn)
On input: initial data at time tstart.
On output: final solution at time tend. q(i,m) = value of mth component in the i'th cell.
Values within the physical domain are in q(i,m) for i = 1,2,...,mx mbc extra cells on each end are needed for boundary conditions as specified in the routine bc1.

method(1:7)   :   array of values specifying the numerical method to use

method(1) = 0 if fixed size time steps are to be taken. In this case, dt = dtv(1) in all steps.
= 1 if variable time steps are to be used. In this case, dt = dtv(1) in the first step and thereafter the value cflv(2) is used to choose the next time step based on the maximum wave speed seen in the previous step. Note that since this value comes from the previous step, the Courant number will not in general be exactly equal to the desired value If the actual Courant number in the next step is greater than 1, then this step is redone with a smaller dt.

The recommended choice of methods for most problems is method(1) = 1, method(2) = 2.


aux(1-mbc:maxmx+mbc, maux)   :   array of auxiliary variables that are used in specifying the problem.
If method(7) = 0 then there are no auxiliary variables and aux can be a dummy variable.
If method(7) = maux > 0 then there are maux auxiliary variables and aux must be dimensioned as above.
Capacity functions are one particular form of auxiliary variable. These arise in some applications, e.g. variable coefficients in advection or acoustics problems. See Clawpack Note # 5 for examples.
If method(6) = 0 then there is no capacity function.
If method(6) = mcapa > 0 then there is a capacity function and capa(i), the "capacity" of the ith cell, is assumed to be stored in aux(i,mcapa). In this case we require method(7).ge.mcapa.

dx : grid spacing in x. (for a computation in ax <= x <= bx, set dx = (bx-ax)/mx.)

tstart : initial time.

tend desired final time (on input) ; actual time reached (on output).

dtv(1:5)   :   array of values related to the time step: (Note: method(1)=1 indicates variable size time steps)

cflv(1:4)  : array of values related to Courant number:

nv(1:2)  :   array of values related to the number of time steps:

nv(1) = unused if method(1) = 0
= maximum number of time steps allowed if method(1) = 1
nv(2) = number of time steps taken (on output).

mthlim(1:mwaves)   :   array of values specifying the flux limiter to be used in each wave family mw. Often the same value will be used for each value of mw, but in some cases it may be desirable to use different limiters. For example, for the Euler equations the superbee limiter might be used for the contact discontinuity (mw=2) while another limiter is used for the nonlinear waves. Several limiters are built in and others can be added by modifying the subroutine philim.

work(mwork)   :   double precision work array of length at least mwork

mwork      :     length of work array. Must be at least (maxmx + 2*mbc) * (2 + 3*meqn + mwaves + meqn*mwaves) If mwork is too small then the program returns with info = 4 and prints the necessary value of mwork to unit 6.

info   :   output value yielding error information:

User-supplied subroutines

bc1   :   subroutine that specifies the boundary conditions. This subroutine should extend the values of q from cells 1:mx to the mbc ghost cells along each edge of the domain.

The form of this subroutine is

rp1    :   user-supplied subroutine that implements the Riemann solver

The form of this subroutine is

On input, ql contains the state vector at the left edge of each cell , qr contains the state vector at the right edge of each cell , auxl contains auxiliary values at the left edge of each cell and auxr contains auxiliary values at the right edge of each cell.

Note that the ith Riemann problem has left state qr(i-1,:) and right state ql(i,:) In the standard clawpack routines, this Riemann solver is called with ql=qr=q along this slice. More flexibility is allowed in case the user wishes to implement another solution method that requires left and rate states at each interface. c If method(7)=maux > 0 then the auxiliary variables along this slice are passed in using auxl and auxr. Again, in the standard routines auxl=auxr=aux in the call to rp1.

On output, wave(i,m,mw) is the m'th component of the jump across wave number mw in the ith Riemann problem. s(i,mw) is the wave speed of wave number mw in the ith Riemann problem. amdq(i,m) = m'th component of A^- Delta q, apdq(i,m) = m'th component of A^+ Delta q, the decomposition of the flux difference f(qr(i-1)) - f(ql(i)) into leftgoing and rightgoing parts respectively.

It is assumed that each wave consists of a jump discontinuity propagating at a single speed, as results, for example, from a Roe approximate Riemann solver. An entropy fix can be included into the specification of amdq and apdq.

src1   :   subroutine for the source terms that solves the equation capa * q_t = psi over time dt.

If method(5)=0 then the equation does not contain a source term and this routine is never called. A dummy argument can be used with many compilers, or provide a dummy subroutine that does nothing (such a subroutine can be found in clawpack/1d/misc/src1xx.f)

The form of this subroutine is

If method(7)=0 or the auxiliary variables are not needed in this solver, then the latter dimension statement can be omitted, but aux should still appear in the argument list.

On input, q(i,m) contains the data for solving the source term equation. On output, q(i,m) should have been replaced by the solution to the source term equation after a step of length dt.

Beginning of claw1 code

implicit double precision (a-h,o-z)
dimension q(1-mbc:maxmx+mbc, meqn)
dimension aux(1-mbc:maxmx+mbc, *)
dimension work(mwork)
dimension mthlim(mwaves),method(7),dtv(5),cflv(4),nv(2)

info = 0
t = tstart
maxn = nv(1)
dt = dtv(1)
cflmax = 0.d0
dtmin = dt
dtmax = dt
nv(2) = 0

if (mx .gt. maxmx) then

info = 1
go to 900


if (method(1) .eq. 0) then fixed size time steps. Compute the number of steps:

maxn = (tend - tstart + 1d-10) / dt

if (dabs(maxn*dt - (tend-tstart)) .gt. 1d-8) then
dt doesn't divide time interval integer number of times
info = 2
go to 900


if (method(1).eq.1 .and. cflv(1).gt.1.d0) then

info = 3
go to 900



i0f = 1
i0wave = i0f + (maxmx + 2*mbc) * meqn
i0s = i0wave + (maxmx + 2*mbc) * meqn * mwaves
i0dtdx = i0s + (maxmx + 2*mbc) * mwaves
i0qwork = i0dtdx + (maxmx + 2*mbc)
i0amdq = i0qwork + (maxmx + 2*mbc) * meqn
i0apdq = i0amdq + (maxmx + 2*mbc) * meqn
i0dtdx = i0apdq + (maxmx + 2*mbc) * meqn
i0end = i0dtdx + (maxmx + 2*mbc) - 1

if (mwork .lt. i0end) then

write(6,*) 'mwork must be increased to ',i0end
info = 4
go to 900



if (maxn.eq.0)

go to 900

do 100 n=1,maxn

told = t

if (told+dt .gt. tend)
dt = tend - told

if (method(1).eq.1) then
save old q in case we need to retake step with smaller dt:
call copyq1(maxmx,meqn,mbc,mx,q,work(i0qwork))

40 continue
dt2 = dt / 2.d0
thalf = t + dt2 midpoint in time for Strang splitting
t = told + dt time at end of step

100 continue
900 continue

if (method(1).eq.1 .and. .and. nv(2) .eq. maxn) then
too many timesteps
info = 11

if (method(1).eq.0 .and. cflmax .gt. cflv(1)) then
Courant number too large with fixed dt
info = 12

tend = t
cflv(3) = cflmax
cflv(4) = cfl
dtv(3) = dtmin
dtv(4) = dtmax
dtv(5) = dt


   subroutine step1(maxmx,meqn,mwaves,mbc,mx,q,aux,dx,dt, method,mthlim,cfl,f,wave,s,amdq,apdq,dtdx,rp1)


Take one time step, updating q.

method(1) = 1 ==> Godunov method
method(1) = 2 ==> Slope limiter
method mthlim(p) controls what limiter is used in the pth family


amdq(1-mbc:maxmx+mbc, meqn)   :   left-going flux-differences
apdq(1-mbc:maxmx+mbc, meqn)    :   right-going flux-differences e.g.
amdq(i,m)     :  m'th component of A^- \Delta q from ith Riemann problem (between cells i-1 and i).

wave(1-mbc:maxmx+mbc, meqn, mwaves)     :    waves from solution of Riemann problems
wave(i,m,mw)      :     mth component of jump in q across wave in family mw in Riemann problem between states i-1 and i.

s(1-mbc:maxmx+mbc, mwaves) :  wave speeds, s(i,mw)
                                                               :  speed of wave in family mw in Riemann problem between states i-1 and i.

f(1-mbc:maxmx+mbc, meqn)      :   correction fluxes for second order method f(i,m)
                                                               :   mth component of flux at left edge of ith cell

Beginning of the code

implicit double precision (a-h,o-z)
dimension q(1-mbc:maxmx+mbc, meqn)
dimension aux(1-mbc:maxmx+mbc, *)
dimension f(1-mbc:maxmx+mbc, meqn)
dimension s(1-mbc:maxmx+mbc, mwaves)
dimension wave(1-mbc:maxmx+mbc, meqn, mwaves)
dimension amdq(1-mbc:maxmx+mbc, meqn)
dimension apdq(1-mbc:maxmx+mbc, meqn)
dimension dtdx(1-mbc:maxmx+mbc)
dimension method(7),mthlim(mwaves)
logical limit

limit = .false.

do 5 mw=1,mwaves

if (mthlim(mw) .gt. 0) limit = .true.

5 continue

mcapa = method(6)

do 10 i=1-mbc,mx+mbc

if ( then

if (aux(i,mcapa) .le. 0.d0) then
write(6,*) 'Error -- capa must be positive'

dtdx(i) = dt / (dx*aux(i,mcapa))


dtdx(i) = dt/dx


10 continue


call rp1(maxmx,meqn,mwaves,mbc,mx,q,q,aux,aux,wave,s,amdq,apdq)


Modify q for Godunov update:
Note this may not correspond to a conservative flux-differencing for equations not in conservation form. It is conservative if amdq + apdq = f(q(i)) - f(q(i-1)).

do 40 i=1,mx+1

do 40 m=1,meqn
q(i,m) = q(i,m) - dtdx(i)*apdq(i,m)
q(i-1,m) = q(i-1,m) - dtdx(i-1)*amdq(i,m)

40 continue


Compute maximum wave speed:
cfl = 0.d0
do 50 i=1,mx+1
cfl = dmax1(cfl, dtdx(i)*dabs(s(i,1)), dtdx(i)*dabs(s(i,mwaves)))
50 continue

if (method(2) .eq. 1)

go to 900

Compute correction fluxes for second order q_{xx} terms:
do 100 m = 1, meqn
do 100 i = 1-mbc, mx+mbc
f(i,m) = 0.d0
100 continue

Apply limiter to waves: if (limit) call limiter(maxmx,meqn,mwaves,mbc,mx,wave,s,mthlim)
do 120 i=1,mx+1
do 120 m=1,meqn
do 110 mw=1,mwaves
dtdxave = 0.5d0 * (dtdx(i-1) + dtdx(i))
f(i,m) = f(i,m) + 0.5d0 * dabs(s(i,mw)) * (1.d0 - dabs(s(i,mw))*dtdxave) * wave(i,m,mw)


Third order corrections:
(still experimental... works well for smooth solutions with no limiters but not well with limiters so far)

if (method(2).lt.3)

go to 110
if (s(i,mw) .gt. 0.d0) then
dq2 = wave(i,m,mw) - wave(i-1,m,mw)
dq2 = wave(i+1,m,mw) - wave(i,m,mw)
f(i,m) = f(i,m) - s(i,mw)/6.d0 * (1.d0 - (s(i,mw)*dtdxave)**2) * dq2

110 continue
120 continue
140 continue


Update q by differencing correction fluxes
(Note: Godunov update has already been performed above)

do 150 m=1,meqn
do 150 i=1,mx
q(i,m) = q(i,m) - dtdx(i) * (f(i+1,m) - f(i,m))
150 continue
900 continue