Merge pull request #349 from donnaaboise/fix-clawpack4-limiters

Minor updates to 4.x limiter;  flux2

src/solvers/fc2d_clawpack4.6/CMakeLists.txt

@@ -6,7 +6,6 @@ add_library(clawpack4.6_f OBJECT
   fortran_source/clawpack46_block.f
   fortran_source/clawpack46_bc2_default.f
   fortran_source/clawpack46_flux2.f
-  fortran_source/clawpack46_flux2fw.f
   fortran_source/clawpack46_step2.f
   fortran_source/clawpack46_step2_wrap.f
 )

src/solvers/fc2d_clawpack4.6/Makefile.am

@@ -16,7 +16,6 @@ libclawpack4_6_compiled_sources = \
 	src/solvers/fc2d_clawpack4.6/fortran_source/clawpack46_block.f \
 	src/solvers/fc2d_clawpack4.6/fortran_source/clawpack46_bc2_default.f \
 	src/solvers/fc2d_clawpack4.6/fortran_source/clawpack46_flux2.f \
-	src/solvers/fc2d_clawpack4.6/fortran_source/clawpack46_flux2fw.f \
 	src/solvers/fc2d_clawpack4.6/fortran_source/clawpack46_step2.f \
 	src/solvers/fc2d_clawpack4.6/fortran_source/clawpack46_step2_wrap.f
 

src/solvers/fc2d_clawpack4.6/fortran_source/clawpack46_flux2.f

@@ -56,7 +56,7 @@ subroutine clawpack46_flux2(ixy,maxm,meqn,maux,mbc,mx,
 c
       implicit none
 
-      integer ixy, maxm, meqn, mbc, mx, maux, use_fwaves, ierror
+      integer ixy, maxm, meqn, mbc, mx, maux, use_fwaves
       integer mwaves, mcapa, method(7), mthlim(mwaves)
 
       external rpn2, rpt2
@@ -113,8 +113,9 @@ subroutine clawpack46_flux2(ixy,maxm,meqn,maux,mbc,mx,
 c     # solve Riemann problem at each interface and compute Godunov updates
 c     ---------------------------------------------------------------------
 c
+c     !! # We pass in maux, even though the typdef does not require it.      
       call rpn2(ixy,maxm,meqn,mwaves,mbc,mx,q1d,q1d,
-     &          aux2,aux2,wave,s,amdq,apdq)
+     &          aux2,aux2,wave,s,amdq,apdq,maux)
 
 c
 c     # Set fadd for the donor-cell upwind method (Godunov)
@@ -197,9 +198,11 @@ subroutine clawpack46_flux2(ixy,maxm,meqn,maux,mbc,mx,
 
 
 c         # split the left-going flux difference into down-going and up-going:
+c         # !! Pass in maux, even though it is not in typedef. This works only
+c         # !! because maux is last argument
           ilr = 1
           call rpt2(ixy,maxm,meqn,mwaves,mbc,mx, q1d,q1d,
-     &          aux1,aux2,aux3, ilr,amdq,bmasdq,bpasdq)
+     &          aux1,aux2,aux3, ilr,amdq,bmasdq,bpasdq,maux)
 
 c         # modify flux below and above by B^- A^- Delta q and  B^+ A^- Delta q:
           do m=1,meqn
@@ -218,7 +221,7 @@ subroutine clawpack46_flux2(ixy,maxm,meqn,maux,mbc,mx,
 c              -going:
           ilr = 2
           call rpt2(ixy,maxm,meqn,mwaves,mbc,mx, q1d,q1d,
-     &               aux1,aux2,aux3, ilr,apdq,bmasdq,bpasdq)
+     &               aux1,aux2,aux3, ilr,apdq,bmasdq,bpasdq,maux)
 
 c         # modify flux below and above by B^- A^+ Delta q and  B^+ A^+ Delta q:
           do m=1,meqn

src/solvers/fc2d_clawpack4.6/fortran_source/clawpack46_inlinelimiter.f

@@ -30,31 +30,40 @@ subroutine clawpack46_inlinelimiter(maxm,meqn,mwaves,mbc,
 c     # given by dotl/wnorm2 and dotr/wnorm2, where wnorm2 is the 2-norm
 c     # of wave.
 c
-      implicit double precision (a-h,o-z)
-      dimension mthlim(mwaves)
-      dimension wave(1-mbc:maxm+mbc,meqn,mwaves)
-      dimension    s(1-mbc:maxm+mbc,mwaves)
+      implicit none
+
+      integer maxm, meqn,mwaves,mbc, mx
+      integer mthlim(mwaves)
+      double precision wave(1-mbc:maxm+mbc,meqn,mwaves)
+      double precision    s(1-mbc:maxm+mbc,mwaves)
+
+      integer mw, i, m
+      double precision wnorm2, dotl, dotr
+      double precision r, wlimitr, th, c
 c
 c
-      do 200 mw=1,mwaves
-         if (mthlim(mw) .eq. 0) go to 200
+      do mw=1,mwaves
+         if (mthlim(mw) .eq. 0) then
+            cycle
+         end if
          dotr = 0.d0
-         do 190 i = 0, mx+1
+         do i = 0, mx+1
             wnorm2 = 0.d0
             dotl = dotr
             dotr = 0.d0
-            do 5 m=1,meqn
+            do m=1,meqn
                wnorm2 = wnorm2 + wave(i,m,mw)**2
                dotr = dotr + wave(i,m,mw)*wave(i+1,m,mw)
-    5          continue
-            if (i .eq. 0) go to 190
-            if (wnorm2 .eq. 0.d0) go to 190
+            end do
+            if (i .eq. 0 .or. wnorm2 .eq. 0) then
+               cycle
+            endif
 c
             if (s(i,mw) .gt. 0.d0) then
                 r = dotl / wnorm2
-              else
+            else
                 r = dotr / wnorm2
-              endif
+            endif
 c
             go to (10,20,30,40,50,60) mthlim(mw)
 c
@@ -98,7 +107,7 @@ subroutine clawpack46_inlinelimiter(maxm,meqn,mwaves,mbc,
 c           ------------------------------
 c           # Generalized minmod
 c           ------------------------------
-            th = 1.3
+            th = 1.3d0
             wlimitr = dmax1(0.d0, dmin1(th*r,(1+r)/2.d0,th));
             go to 170
 c
@@ -107,12 +116,12 @@ subroutine clawpack46_inlinelimiter(maxm,meqn,mwaves,mbc,
 c
 c           # apply limiter to waves:
 c
-            do 180 m=1,meqn
+            do m=1,meqn
                wave(i,m,mw) = wlimitr * wave(i,m,mw)
-  180          continue
+            end do
+
+         end do
+      end do
 
-  190       continue
-  200    continue
-c
       return
       end