update rpn2_geoclaw, geoclaw_riemann_utils to initialize wave compone…

…nts 4:5 to 0 for Boussinesq code

src/geoclaw_riemann_utils.f

@@ -265,6 +265,8 @@ subroutine riemann_aug_JCP(maxiter,meqn,mwaves,hL,hR,huL,huR,
 
  enddo ! end iteration on Riemann problem
 
+ fw(:,:) = 0.d0 ! initialize before setting some waves
+
  do mw=1,mwaves
  sw(mw)=lambda(mw)
  fw(1,mw)=beta(mw)*r(2,mw)

src/rpn2_geoclaw.f

@@ -56,7 +56,7 @@ subroutine rpn2(ixy,maxm,meqn,mwaves,maux,mbc,mx,
  !local only
  integer m,i,mw,maxiter,mu,nv
  double precision wall(3)
- double precision fw(3,3)
+ double precision fw(meqn,3)
  double precision sw(3)
 
  double precision hR,hL,huR,huL,uR,uL,hvR,hvL,vR,vL,phiR,phiL,pL,pR
@@ -65,11 +65,21 @@ subroutine rpn2(ixy,maxm,meqn,mwaves,maux,mbc,mx,
  double precision hstar,hstartest,hstarHLL,sLtest,sRtest
  double precision tw,dxdc
 
+ logical :: debug
+
  logical rare1,rare2
 
  ! In case there is no pressure forcing
  pL = 0.d0
  pR = 0.d0
+ debug = .false.
+
+ ! initialize all components to 0
+ fw(:,:) = 0.d0
+ fwave(:,:,:) = 0.d0
+ s(:,:) = 0.d0
+ amdq(:,:) = 0.d0
+ apdq(:,:) = 0.d0
 
  !loop through Riemann problems at each grid cell
  do i=2-mbc,mx+mbc
@@ -80,13 +90,6 @@ subroutine rpn2(ixy,maxm,meqn,mwaves,maux,mbc,mx,
  write(*,*) 'Negative input: hl,hr,i=',qr(1,i-1),ql(1,i),i
  endif
 
- !Initialize Riemann problem for grid interface
- do mw=1,mwaves
- s(mw,i)=0.d0
- fwave(1,mw,i)=0.d0
- fwave(2,mw,i)=0.d0
- fwave(3,mw,i)=0.d0
- enddo
 
 c !set normal direction
  if (ixy.eq.1) then
@@ -213,7 +216,7 @@ subroutine rpn2(ixy,maxm,meqn,mwaves,maux,mbc,mx,
 
  maxiter = 1
 
- call riemann_aug_JCP(maxiter,3,3,hL,hR,huL,
+ call riemann_aug_JCP(maxiter,meqn,mwaves,hL,hR,huL,
  & huR,hvL,hvR,bL,bR,uL,uR,vL,vR,phiL,phiR,pL,pR,sE1,sE2,
  & drytol,g,rho,sw,fw)
 
@@ -273,8 +276,7 @@ subroutine rpn2(ixy,maxm,meqn,mwaves,maux,mbc,mx,
 
 
 c============= compute fluctuations=============================================
- amdq(1:3,:) = 0.d0
- apdq(1:3,:) = 0.d0
+
  do i=2-mbc,mx+mbc
  do mw=1,mwaves
  if (s(mw,i) < 0.d0) then
@@ -287,14 +289,17 @@ subroutine rpn2(ixy,maxm,meqn,mwaves,maux,mbc,mx,
  endif
  enddo
  enddo
-!-- do i=2-mbc,mx+mbc
-!-- do m=1,meqn
-!-- write(51,151) m,i,amdq(m,i),apdq(m,i)
-!-- write(51,152) fwave(m,1,i),fwave(m,2,i),fwave(m,3,i)
-!--151 format("++3 ampdq ",2i4,2e25.15)
-!--152 format("++3 fwave ",8x,3e25.15)
-!-- enddo
-!-- enddo
+
+ if (debug) then
+ do i=2-mbc,mx+mbc
+ do m=1,meqn
+ write(51,151) m,i,amdq(m,i),apdq(m,i)
+ write(51,152) fwave(m,1,i),fwave(m,2,i),fwave(m,3,i)
+ 151 format("+++ ampdq ",2i4,2e25.15)
+ 152 format("+++ fwave ",8x,3e25.15)
+ enddo
+ enddo
+ endif
 
  return
  end subroutine