Merge pull request #111 from rjleveque/geoclaw_riemann_utils_update

Update to JCP GeoClaw Riemann solver

src/geoclaw_riemann_utils.f

@@ -36,6 +36,7 @@ subroutine riemann_aug_JCP(maxiter,meqn,mwaves,hL,hR,huL,huR,
       double precision delh,delhu,delphi,delb,delnorm
       double precision rare1st,rare2st,sdelta,raremin,raremax
       double precision criticaltol,convergencetol,raretol
+      double precision criticaltol_2, hustar_interface
       double precision s1s2bar,s1s2tilde,hbar,hLstar,hRstar,hustar
       double precision huRstar,huLstar,uRstar,uLstar,hstarHLL
       double precision deldelh,deldelphi
@@ -109,7 +110,10 @@ subroutine riemann_aug_JCP(maxiter,meqn,mwaves,hL,hR,huL,huR,
 
 
 c     !determine the steady state wave -------------------
-      criticaltol = 1.d-6
+      !criticaltol = 1.d-6
+      ! MODIFIED:
+      criticaltol = max(drytol*g, 1d-6)
+      criticaltol_2 = sqrt(criticaltol)
       deldelh = -delb
       deldelphi = -g*0.5d0*(hR+hL)*delb
 
@@ -145,15 +149,25 @@ subroutine riemann_aug_JCP(maxiter,meqn,mwaves,hL,hR,huL,huR,
          s1s2tilde= max(0.d0,uLstar*uRstar) - g*hbar
 
 c        !find if sonic problem
-         sonic=.false.
-         if (abs(s1s2bar).le.criticaltol) sonic=.true.
-         if (s1s2bar*s1s2tilde.le.criticaltol) sonic=.true.
-         if (s1s2bar*sE1*sE2.le.criticaltol) sonic = .true.
-         if (min(abs(sE1),abs(sE2)).lt.criticaltol) sonic=.true.
-         if (sE1.lt.0.d0.and.s1m.gt.0.d0) sonic = .true.
-         if (sE2.gt.0.d0.and.s2m.lt.0.d0) sonic = .true.
-         if ((uL+sqrt(g*hL))*(uR+sqrt(g*hR)).lt.0.d0) sonic=.true.
-         if ((uL-sqrt(g*hL))*(uR-sqrt(g*hR)).lt.0.d0) sonic=.true.
+         ! MODIFIED from 5.3.1 version
+         sonic = .false.
+         if (abs(s1s2bar) <= criticaltol) then
+            sonic = .true.
+         else if (s1s2bar*s1s2tilde <= criticaltol**2) then
+            sonic = .true.
+         else if (s1s2bar*sE1*sE2 <= criticaltol**2) then
+            sonic = .true.
+         else if (min(abs(sE1),abs(sE2)) < criticaltol_2) then
+            sonic = .true.
+         else if (sE1 <  criticaltol_2 .and. s1m > -criticaltol_2) then
+            sonic = .true.
+         else if (sE2 > -criticaltol_2 .and. s2m <  criticaltol_2) then
+            sonic = .true.
+         else if ((uL+dsqrt(g*hL))*(uR+dsqrt(g*hR)) < 0.d0) then
+            sonic = .true.
+         else if ((uL- dsqrt(g*hL))*(uR- dsqrt(g*hR)) < 0.d0) then
+            sonic = .true.
+         end if
 
 c        !find jump in h, deldelh
          if (sonic) then
@@ -254,9 +268,18 @@ subroutine riemann_aug_JCP(maxiter,meqn,mwaves,hL,hR,huL,huR,
          fw(3,mw)=beta(mw)*r(2,mw)
       enddo
       !find transverse components (ie huv jumps).
+      ! MODIFIED from 5.3.1 version
       fw(3,1)=fw(3,1)*vL
       fw(3,3)=fw(3,3)*vR
-      fw(3,2)= hR*uR*vR - hL*uL*vL - fw(3,1)- fw(3,3)
+      fw(3,2)= 0.d0
+ 
+      hustar_interface = huL + fw(1,1)   ! = huR - fw(1,3)
+      if (hustar_interface <= 0.0d0) then
+          fw(3,1) = fw(3,1) + (hR*uR*vR - hL*uL*vL - fw(3,1)- fw(3,3))
+        else
+          fw(3,3) = fw(3,3) + (hR*uR*vR - hL*uL*vL - fw(3,1)- fw(3,3))
+        end if
+
 
       return
 
@@ -623,4 +646,3 @@ subroutine riemanntype(hL,hR,uL,uR,hm,s1m,s2m,rare1,rare2,
       return
 
       end ! subroutine riemanntype----------------------------------------------------------------
-