c c c ===================================================== subroutine rpn2swt(ixy,maxm,meqn,mwaves,mbc,mx,ql,qr,auxl,auxr, & wave,s,amdq,apdq) c ===================================================== c c # Roe-solver for the 2D shallow water equations with bottom topography. c # This implements the quasi-steady wave-propagation method for c # handling the source terms, as described in the paper c # ftp://amath.washington.edu/pub/rjl/papers/qsteady.ps.gz c c # This version also enforces zero-order extrapolation boundary c # conditions by setting all wave-strengths to zero at the boundaries. c c # solve Riemann problems along one slice of data. 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 c # This data is along a slice in the x-direction if ixy=1 c # or the y-direction if ixy=2. c # On output, wave contains the waves, s the speeds, c # and amdq, apdq the decomposition of the flux difference c # f(qr(i-1)) - f(ql(i)) c # into leftgoing and rightgoing parts respectively. c # With the Roe solver we have c # amdq = A^- \Delta q and apdq = A^+ \Delta q c # where A is the Roe matrix. An entropy fix can also be incorporated c # into the flux differences. 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 routines, this routine is called with ql = qr c c # The bottom topography H(x,y) is stored in the aux array. c # If (x(i),y(j)) represents cell centers. Then the aux array should c # be set so that c # aux(i,j,1) = H(x(i)-dx/2, y(j)) (Value of H at middle of left edge) c # aux(i,j,2) = H(x(i), y(j)-dy/2) (Value of H at middle of bottom edge) c c implicit double precision (a-h,o-z) c dimension wave(1-mbc:maxm+mbc, meqn, mwaves) dimension s(1-mbc:maxm+mbc, mwaves) dimension ql(1-mbc:maxm+mbc, meqn) dimension qr(1-mbc:maxm+mbc, meqn) dimension auxl(1-mbc:maxm+mbc, *) dimension auxr(1-mbc:maxm+mbc, *) dimension apdq(1-mbc:maxm+mbc, meqn) dimension amdq(1-mbc:maxm+mbc, meqn) c c local arrays -- common block comroe is passed to rpt2sh c ------------ parameter (maxm2 = 603) !# assumes at most 600x600 grid with mbc=3 dimension delta(3) logical efix common /param/ grav common /comroe/ u(-2:maxm2),v(-2:maxm2),a(-2:maxm2),hl(-2:maxm2), & hr(-2:maxm2) c data efix /.true./ !# use entropy fix for transonic rarefactions c if (-2.gt.1-mbc .or. maxm2 .lt. maxm+mbc) then write(6,*) 'Check dimensions of local arrays in rpn2' stop endif c c # set mu to point to the component of the system that corresponds c # to momentum in the direction of this slice, mv to the orthogonal c # momentum: c if (ixy.eq.1) then mu = 2 mv = 3 else mu = 3 mv = 2 endif c c # note that notation for u and v reflects assumption that the c # Riemann problems are in the x-direction with u in the normal c # direction and v in the orthogonal direcion, but with the above c # definitions of mu and mv the routine also works with ixy=2 c # and returns, waves, speeds, and flux differences for the c # Riemann problems u_t + g(u)_y = 0 in the y-direction. c c c # The constant depth h = ql(i,1) in the i'th cell is replaced by c # two values hl and hr with a jump discontinuity, chosen in such c # a way that the flux difference across this jump balances the c # source term arising from the bottom topography. c c # assumes ql = qr on input c c do 10 i=1-mbc,mx+mbc hu = ql(i,mu) DH = auxl(i+1,ixy) - auxl(i,ixy) 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 8 continue hl(i) = ql(i,1) - delh hr(i) = ql(i,1) + delh 10 continue c c # compute the Roe-averaged variables needed in the Roe solver. c # These are stored in the common block comroe since they are c # later used in routine rpt2sh to do the transverse wave splitting. c do 40 i = 2-mbc, mx+mbc hsqrtl = dsqrt(hr(i-1)) hsqrtr = dsqrt(hl(i)) hsq2 = hsqrtl + hsqrtr u(i) = (qr(i-1,mu)/hsqrtl + ql(i,mu)/hsqrtr) / hsq2 v(i) = (qr(i-1,mv)/hsqrtl + ql(i,mv)/hsqrtr) / hsq2 a(i) = dsqrt(grav*0.5d0*(hr(i-1)+hl(i))) 40 continue c c c # now split the jump in q at each interface into waves c c # find a1 thru a3, the coefficients of the 3 eigenvectors: do 50 i = 2-mbc, mx+mbc delta(1) = hl(i) - hr(i-1) delta(2) = ql(i,mu) - qr(i-1,mu) delta(3) = ql(i,mv) - qr(i-1,mv) a1 = ((u(i)+a(i))*delta(1) - delta(2))*(0.50d0/a(i)) a2 = -v(i)*delta(1) + delta(3) a3 = (-(u(i)-a(i))*delta(1) + delta(2))*(0.50d0/a(i)) c c # Compute the waves. c wave(i,1,1) = a1 wave(i,mu,1) = a1*(u(i)-a(i)) wave(i,mv,1) = a1*v(i) s(i,1) = u(i)-a(i) c wave(i,1,2) = 0.0d0 wave(i,mu,2) = 0.0d0 wave(i,mv,2) = a2 s(i,2) = u(i) c wave(i,1,3) = a3 wave(i,mu,3) = a3*(u(i)+a(i)) wave(i,mv,3) = a3*v(i) s(i,3) = u(i)+a(i) 50 continue c c # for extrapolation boundary conditions, want zero-strength waves at the c # boundaries. This section must be changed to use other boundary c # conditions. c do 51 i=2-mbc,1 do mw=1,3 s(i,mw) = 0.d0 do m=1,3 wave(i,m,mw) = 0.d0 enddo enddo 51 continue c do 52 i=mx+1,mx+mbc do mw=1,3 s(i,mw) = 0.d0 do m=1,3 wave(i,m,mw) = 0.d0 enddo enddo 52 continue c c c # compute flux differences amdq and apdq. c --------------------------------------- c if (efix) go to 110 c c # no entropy fix c ---------------- c c # amdq = SUM s*wave over left-going waves c # apdq = SUM s*wave over right-going waves c do 100 m=1,3 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 go to 900 c c----------------------------------------------------- c 110 continue c c # With entropy fix c ------------------ c c # compute flux differences amdq and apdq. c # First compute amdq as sum of s*wave for left going waves. c # Incorporate entropy fix by adding a modified fraction of wave c # if s should change sign. c do 200 i=2-mbc,mx+mbc c check 1-wave him1 = hr(i-1) s0 = qr(i-1,mu)/him1 - dsqrt(grav*him1) c check for fully supersonic case : if (s0.gt.0.0d0.and.s(i,1).gt.0.0d0) then do 60 m=1,3 amdq(i,m)=0.0d0 60 continue goto 200 endif c h1 = hr(i-1)+wave(i,1,1) hu1= qr(i-1,mu)+wave(i,mu,1) s1 = hu1/h1 - dsqrt(grav*h1) !speed just to right of 1-wave if (s0.lt.0.0d0.and.s1.gt.0.0d0) then c transonic rarefaction in 1-wave sfract = s0*((s1-s(i,1))/(s1-s0)) else if (s(i,1).lt.0.0d0) then c 1-wave is leftgoing sfract = s(i,1) else c 1-wave is rightgoing sfract = 0.0d0 endif do 120 m=1,3 amdq(i,m) = sfract*wave(i,m,1) 120 continue c check 2-wave if (s(i,2).gt.0.0d0) go to 200 !2 and 3 waves are right do 140 m=1,3 amdq(i,m) = amdq(i,m) + s(i,2)*wave(i,m,2) 140 continue c c check 3-wave c hi = hl(i) s03 = ql(i,mu)/hi + dsqrt(grav*hi) h3=hl(i)-wave(i,1,3) hu3=ql(i,mu)-wave(i,mu,3) s3=hu3/h3 + dsqrt(grav*h3) if (s3.lt.0.0d0.and.s03.gt.0.0d0) then c transonic rarefaction in 2-wave sfract = s3*((s03-s(i,3))/(s03-s3)) else if (s(i,3).lt.0.0d0) then c 3-wave is leftgoing sfract = s(i,3) else c 3-wave is rightgoing goto 200 endif do 160 m=1,3 amdq(i,m) = amdq(i,m) + sfract*wave(i,m,3) 160 continue 200 CONTINUE c c compute rightgoing flux differences : c do 220 m=1,3 do 220 i = 2-mbc,mx+mbc df = 0.0d0 do 210 mw=1,mwaves df = df + s(i,mw)*wave(i,m,mw) 210 continue apdq(i,m)=df-amdq(i,m) 220 continue c c 900 Continue return end