c ========================================================= subroutine rp1swt(maxmx,meqn,mwaves,mbc,mx,ql,qr,auxl,auxr, & wave,s,amdq,apdq) c ========================================================= c c # solve Riemann problems for the 1D shallow water equations c # (h)_t + (u h)_x = 0 c # (uh)_t + ( uuh + .5*gh^2 )_x = - g*h*H' c # using Roe's approximate Riemann solver with entropy fix for c # transonic rarefractions. c c # With modification of h to cancel source terms. c c # On input, ql contains the state vector at the left edge of each cell c # qr contains the state vector at the right edge of each cell c # On output, wave contains the waves, c # s the speeds, c # amdq the left-going flux difference A^- \Delta q c # apdq the right-going flux difference A^+ \Delta q c c # Note that the i'th Riemann problem has left state qr(i-1,:) c # and right state ql(i,:) c # From the basic clawpack routine step1, rp is called with ql = qr = q. c c Here meqn=mwaves=2 should be passed from the calling routine c implicit double precision (a-h,o-z) dimension ql(1-mbc:maxmx+mbc, meqn) dimension qr(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 auxl(1-mbc:maxmx+mbc, 1) dimension auxr(1-mbc:maxmx+mbc, 1) c c # local storage c --------------- dimension delta(2) parameter (maxmrp = 2002) dimension hl(-1:maxmrp),hr(-1:maxmrp) common /param/ grav c do 10 i=1-mbc,mx+mbc hu = ql(i,2) DH = auxl(i+1,1) - auxl(i,1) c c # delta h for the case hu=0: delh = -0.5d0*DH c c # Newton iteration to improve delh: do 5 iter=1,5 hp = ql(i,1) + delh hm = ql(i,1) - delh F = hu*hu*(1.d0/hp - 1.d0/hm) + 0.5d0*grav*(hp*hp - hm*hm) & + grav*ql(i,1)*DH Fprime = -hu*hu*(1.d0/hp**2 + 1.d0/hm**2) + grav*(hp+hm) dnewton = F/Fprime delh = delh - dnewton if (dabs(dnewton).lt.1d-6) go to 8 5 continue write(6,*) 'nonconvergence of newton in rp1swt' write(6,*) ' dnewton =',dnewton c c 8 continue hl(i) = ql(i,1) - delh hr(i) = ql(i,1) + delh 10 continue do 30 i=2-mbc,mx+mbc c c # compute Roe-averaged quantities: ubar = (qr(i-1,2)/dsqrt(hr(i-1)) + ql(i,2)/dsqrt(hl(i)))/ . ( dsqrt(hr(i-1)) + dsqrt(hl(i)) ) cbar=dsqrt(0.5d0*grav*(hr(i-1) + hl(i))) c # delta(1)=h(i)-h(i-1) and delta(2)=hu(i)-hu(i-1) delta(1) = hl(i) - hr(i-1) delta(2) = ql(i,2) - qr(i-1,2) c # compute coeffs in the evector expansion of delta(1),delta(2) a1 = 0.5d0*(-delta(2) + (ubar + cbar) * delta(1))/cbar a2 = 0.5d0*( delta(2) - (ubar - cbar) * delta(1))/cbar c # finally, compute the waves. wave(i,1,1) = a1 wave(i,2,1) = a1*(ubar - cbar) s(i,1) = ubar - cbar wave(i,1,2) = a2 wave(i,2,2) = a2*(ubar + cbar) s(i,2) = ubar + cbar 30 continue c # no entropy fix (need to add still!) c ---------------------------------------------- c # amdq = SUM s*wave over left-going waves c # apdq = SUM s*wave over right-going waves do 100 m=1,2 do 100 i=2-mbc, mx+mbc amdq(i,m) = 0.d0 apdq(i,m) = 0.d0 do 90 mw=1,mwaves if (s(i,mw) .lt. 0.d0) then amdq(i,m) = amdq(i,m) + s(i,mw)*wave(i,m,mw) else apdq(i,m) = apdq(i,m) + s(i,mw)*wave(i,m,mw) endif 90 continue 100 continue c return end