add code for debugging difference between dclaw 4 and dclaw 5

geoclaw/2d/lib/flux2fw_geo.f

@@ -151,6 +151,13 @@ subroutine flux2(ixy,maxm,meqn,maux,mbc,mx,
  faddp(i,m) = faddp(i,m) - apdq(i,m)
  faddm(i,m) = faddm(i,m) + amdq(i,m)
  enddo
+
+ if (ixy.eq.2.and.i.ge.290.and.i.le.297) then
+ write(56,*) "+++ flux2, ixy=2, amdq"
+ write(56,599) i, (amdq(i,mm), mm=1,meqn)
+ 599 format(i5,7e15.7)
+ endif
+
  if (relimit) then
  faddp(i,1) = faddp(i,1) + dxdc*q1d(i,mu)
  faddm(i,1) = faddp(i,1)
@@ -194,6 +201,8 @@ subroutine flux2(ixy,maxm,meqn,maux,mbc,mx,
  119 continue
  faddm(i,m) = faddm(i,m) + 0.5d0 * cqxx(i,m)
  faddp(i,m) = faddp(i,m) + 0.5d0 * cqxx(i,m)
+
+ write(56,*) "++++ second order, shouldn't be here"
  120 continue
 c
 c

geoclaw/2d/lib/step2_geo.f

@@ -55,6 +55,7 @@ subroutine step2(maxm,maxmx,maxmy,meqn,maux,mbc,mx,my,
  dimension work(mwork)
 
  logical relimit
+ integer ii,jj,mm
  common /comxyt/ dtcom,dxcom,dycom,tcom,icom,jcom
  include "call.i"
 c
@@ -164,7 +165,19 @@ subroutine step2(maxm,maxmx,maxmy,meqn,maux,mbc,mx,my,
  gp(i,j,m) = gp(i,j,m) + gaddp(i,m,1)
  gm(i,j+1,m) = gm(i,j+1,m) + gaddm(i,m,2)
  gp(i,j+1,m) = gp(i,j+1,m) + gaddp(i,m,2)
+
  25 continue
+ if (j.eq.289 .or. j.eq.290) then
+ write(55,*) 
+ write(55,*) "+++ x-sweep"
+ write(55,*) "+++ j", j
+ do ii = 1,mx+1
+ write(55,599) ii, (gaddm(ii,mm,1), mm=1,7)
+ 599 format(i5,7e15.7)
+ write(55,599) ii, (gaddm(ii,mm,2), mm=1,7)
+ enddo 
+ endif
+
  50 continue
 c
 c
@@ -223,6 +236,17 @@ subroutine step2(maxm,maxmx,maxmy,meqn,maux,mbc,mx,my,
  fm(i+1,j,m) = fm(i+1,j,m)+ gaddm(j,m,2)
  fp(i+1,j,m) = fp(i+1,j,m)+ gaddp(j,m,2)
  75 continue
+
+
+ if (i.eq.173 .or. i.eq.172) then
+ write(55,*) 
+ write(55,*) "+++ y-sweep"
+ write(55,*) "+++ i", i
+ do jj = 1,my+1
+ write(55,599) jj, (faddm(jj,mm), mm=1,7)
+ enddo
+ endif
+
  100 continue
 
 c # relimit correction fluxes if they drive a cell negative

geoclaw/2d/lib/valout_geo.f

@@ -111,12 +111,12 @@ subroutine valout (lst, lend, time, nvar, naux)
  endif
  surface = alloc(iadd(i,j,1)) + alloc(iaddaux(i,j,1))
  & - alloc(iadd(i,j,7))
- write(matunit1,109) (alloc(iadd(i,j,ivar)), ivar=1,nvar),
+ write(matunit1,109) i, j, (alloc(iadd(i,j,ivar)), ivar=1,1),
  & surface
  enddo
  write(matunit1,*) ' '
  enddo
- 109 format(9e26.16)
+ 109 format(2i5, 9e26.16)
 
 
  mptr = node(levelptr, mptr)
@@ -197,7 +197,7 @@ subroutine valout (lst, lend, time, nvar, naux)
  alloc(iaddaux(i,j,ivar)) = 0.d0
  endif
  enddo
- write(matunit3,109) (alloc(iaddaux(i,j,ivar)),
+ write(matunit3,109) i, j, (alloc(iaddaux(i,j,ivar)),
  & ivar=1,naux)
  enddo
  write(matunit3,*) ' '

geoclaw/2d/lib_dig/digclaw_mod.f90

@@ -487,7 +487,7 @@ subroutine auxeval(h,u,v,m,p,phi_bed,theta,kappa,S,rho,tanpsi,D,tau,sigbed,kperm
  ! RPJ reduced tau_solid by m/mcrit to reflect the loss of granular friction for low values of m.
  ! this regularization may be improved.
 
- tau = dmax1(0.d0,sigbed*tan(atan(mu_bf)+atan(tanpsi)))
+ tau = dmax1(0.d0,sigbed*tan(atan(mu_bf)+atan(tanpsi)))*((tanh(100.0*(m-0.05))+1.0)*0.5)
 
  ! viscosity (depth averaged is vh^(1/2)/2)
  ! in Rocha, Johnson, and Gray, they calculate viscosity.

geoclaw/2d/lib_dig/qinit_geo.f

@@ -43,13 +43,13 @@ subroutine qinit(maxmx,maxmy,meqn,mbc,mx,my,xlower,ylower,
 
  xintlow = dmax1(xlower,xlowqinit(mf))
  xinthi = dmin1(xhigher,xhiqinit(mf))
- istart = max(1-mbc,int(0.5 + (xintlow-xlower)/dx)-mbc)
- iend = min(mx+mbc,int(1.0 + (xinthi-xlower)/dx)+mbc)
+ istart = max(1-mbc,int(0.5d0 + (xintlow-xlower)/dx)-mbc)
+ iend = min(mx+mbc,int(1.0d0 + (xinthi-xlower)/dx)+mbc)
 
  yintlow = dmax1(ylower,ylowqinit(mf))
  yinthi = dmin1(yhigher,yhiqinit(mf))
- jstart = max(1-mbc,int(0.5 + (yintlow-ylower)/dy)-mbc)
- jend = min(my+mbc,int(1.0 + (yinthi-ylower)/dy)+mbc)
+ jstart = max(1-mbc,int(0.5d0 + (yintlow-ylower)/dy)-mbc)
+ jend = min(my+mbc,int(1.0d0 + (yinthi-ylower)/dy)+mbc)
 
  do i=istart,iend
  x = xlower + (i-0.5d0)*dx
@@ -125,7 +125,7 @@ subroutine qinit(maxmx,maxmy,meqn,mbc,mx,my,xlower,ylower,
  q(i,j,5) = q(i,j,1)*q(i,j,5)
  endif
  if (q(i,j,1).le.drytolerance) then
- do m = 2,meqn
+ do m = 1,meqn
  q(i,j,m) = 0.d0
  enddo
  endif

geoclaw/2d/lib_dig/riemannsolvers_geo.f

@@ -68,12 +68,12 @@ subroutine riemann_dig2_aug_sswave_ez(ixy,meqn,mwaves,hL,hR,
  pmL = chiL
  pmR = chiR
  pm = 0.5*(chiL + chiR)
- pm = min(1.0,pm)
- pm = max(0.0,pm)
- if (alpha_seg==1.0) then
- seg = 0.0
+ pm = min(1.0d0,pm)
+ pm = max(0.0d0,pm)
+ if (alpha_seg==1.0d0) then
+ seg = 0.0d0
  else
- seg = 1.0
+ seg = 1.0d0
  endif
  !pmtanh01 = seg*(0.5*(tanh(20.0*(pm-0.80))+1.0))
  !pmtanh01 = seg*(0.5*(tanh(40.0*(pm-0.90))+1.0))
@@ -94,17 +94,17 @@ subroutine riemann_dig2_aug_sswave_ez(ixy,meqn,mwaves,hL,hR,
  & kappa,SN,rhoL,tanpsi,D,tauL,sigbed,kperm,compress,pmL)
  call auxeval(hL,uL,vL,mL,pL,phiL,thetaL,
  & kappa,SN,rhoR,tanpsi,D,tauR,sigbed,kperm,compress,pmL)
- tauR=0.5*tauL
- h = 0.5*hL
+ tauR=0.5d0*tauL
+ h = 0.5d0*hL
  v = vL
  mbar = mL
  else
  call auxeval(hR,uR,vR,mR,pR,phiR,thetaR,
  & kappa,SN,rhoL,tanpsi,D,tauL,sigbed,kperm,compress,pmR)
  call auxeval(hR,uR,vR,mR,pR,phiR,thetaR,
  & kappa,SN,rhoR,tanpsi,D,tauR,sigbed,kperm,compress,pmR)
- tauL=0.5*tauR
- h = 0.5*hR
+ tauL=0.5d0*tauR
+ h = 0.5d0*hR
  v = vR
  mbar = mR
  endif
@@ -262,7 +262,7 @@ subroutine riemann_dig2_aug_sswave_ez(ixy,meqn,mwaves,hL,hR,
 
  vnorm = sqrt(uR**2 + uL**2 + vR**2 + vL**2)
  !vnorm = sqrt(u**2 + v**2)
- if (vnorm>0.0) then
+ if (vnorm>0.0d0) then
  !tausource = - dx*0.5*(tauL/rhoL + tauR/rhoR)*u/vnorm
  !tausource = - dx*max(tauL/rhoL , tauR/rhoR)*u/vnorm
  !tausource = 0.0!dx*tauR*taudir/rhoR
@@ -273,37 +273,37 @@ subroutine riemann_dig2_aug_sswave_ez(ixy,meqn,mwaves,hL,hR,
  ! write(*,*) 'ixy', ixy
  !endif
 
- if ((uL**2 + vL**2)==0.0) then
+ if ((uL**2 + vL**2)==0.0d0) then
  taudirL = taudirR
  else
  taudirL = -uL/sqrt(uL**2 + vL**2)
  endif
 
- if ((uR**2 + vR**2)>0.0) then
+ if ((uR**2 + vR**2)>0.0d0) then
  taudirR = -uR/sqrt(uR**2 + vR**2)
  endif
 
- tausource = 0.0*dx*0.5*(taudirL*tauL/rhoL + tauR*taudirR/rhoR)
+ tausource=0.d0*dx*0.5d0*(taudirL*tauL/rhoL+tauR*taudirR/rhoR)
  !elseif (dx*max(abs(tauL*taudirL/rhoL),abs(tauR*taudirR/rhoR))
- elseif (dx*0.5*abs(taudirR*tauR/rhoR + tauL*taudirR/rhoL)
+ elseif (dx*0.5d0*abs(taudirR*tauR/rhoR + tauL*taudirR/rhoL)
  & .gt.abs(del(2) - source2dx)) then 
 
  !no failure
  tausource = del(2) - source2dx
- del(1) = 0.0
- del(0) = 0.0
- del(4) = 0.0
+ del(1) = 0.0d0
+ del(0) = 0.0d0
+ del(4) = 0.0d0
  else
  !failure
  !tausource = sign(abs(dx*0.5*taudir*(tauL/rhoL + tauR/rhoR))
  !tausource = sign(abs(dx*0.5*tau/rho)
  !tausource = dx*taudir*tauR/rhoR
 
- tausource = 0.5*dx*((taudirR*tauR/rhoR)+(tauL*taudirR/rhoL))
+ tausource = 0.5d0*dx*((taudirR*tauR/rhoR)+(tauL*taudirR/rhoL))
  tausource = dsign(tausource,del(2)-source2dx)
  endif
  if (wallprob) then
- tausource = 0.0
+ tausource = 0.0d0
  endif
  del(2) = del(2) - source2dx - tausource
  !del(4) = del(4) + dx*3.0*vnorm*tanpsi/(h*compress)
@@ -349,6 +349,11 @@ subroutine riemann_dig2_aug_sswave_ez(ixy,meqn,mwaves,hL,hR,
  enddo
  enddo
 
+ if (ixy.eq.2) then 
+ write(58, 699) beta(1:3), del(0:2)
+699 format("beta(1:3), del(0:2)", 6e15.7)
+ endif 
+
  !waves and fwaves for delta hum
  fw(4,1) = fw(1,1)*mL
  fw(4,3) = fw(1,3)*mR
@@ -368,10 +373,10 @@ subroutine riemann_dig2_aug_sswave_ez(ixy,meqn,mwaves,hL,hR,
 
 
  !fwaves for segregation
- seg_L = chiL*hL*uL*(1.0+(1.0-alpha_seg)*(1.0-chiL))
- seg_R = chiR*hR*uR*(1.0+(1.0-alpha_seg)*(1.0-chiR))
- fw(6,1) = fw(1,1)*chiL*(1.0+(1.0-alpha_seg)*(1.0-chiL))
- fw(6,3) = fw(1,3)*chiR*(1.0+(1.0-alpha_seg)*(1.0-chiR))
+ seg_L = chiL*hL*uL*(1.0d0+(1.0d0-alpha_seg)*(1.0d0-chiL))
+ seg_R = chiR*hR*uR*(1.0d0+(1.0d0-alpha_seg)*(1.0d0-chiR))
+ fw(6,1) = fw(1,1)*chiL*(1.0d0+(1.0d0-alpha_seg)*(1.0d0-chiL))
+ fw(6,3) = fw(1,3)*chiR*(1.0d0+(1.0d0-alpha_seg)*(1.0d0-chiR))
  fw(6,2) = seg_R - seg_L - fw(6,1) - fw(6,3)
  return
  end !subroutine riemann_dig2_aug_sswave_ez
@@ -789,8 +794,8 @@ subroutine riemanntype(hL,hR,uL,uR,hm,s1m,s2m,rare1,rare2,
  h0=(dsqrt(h0)-F0/slope)**2
  enddo
  hm=h0
- u1m=uL-(hm-hL)*dsqrt((.5d0*g)*(1/hm + 1/hL))
- u2m=uR+(hm-hR)*dsqrt((.5d0*g)*(1/hm + 1/hR))
+ u1m=uL-(hm-hL)*dsqrt((.5d0*g)*(1.d0/hm + 1.d0/hL))
+ u2m=uR+(hm-hR)*dsqrt((.5d0*g)*(1.d0/hm + 1.d0/hR))
  um=.5d0*(u1m+u2m)
 
  s1m=u1m-dsqrt(g*hm)
@@ -803,7 +808,7 @@ subroutine riemanntype(hL,hR,uL,uR,hm,s1m,s2m,rare1,rare2,
 
  do iter=1,maxiter
  F0=delu + 2.d0*(dsqrt(g*h0)-dsqrt(g*h_max))
- & + (h0-h_min)*dsqrt(.5d0*g*(1/h0+1/h_min))
+ & + (h0-h_min)*dsqrt(.5d0*g*(1.d0/h0+1.d0/h_min))
  slope=(F_max-F0)/(h_max-h_min)
  h0=h0-F0/slope
  enddo
@@ -1150,7 +1155,7 @@ SUBROUTINE gaussj(a,n,np,b,m,mp)
  if (a(icol,icol).eq.0.) write(*,*) 'singular matrix in gaussj'
 
  pivinv=1./a(icol,icol)
- a(icol,icol)=1.
+ a(icol,icol)=1.d0
  do 16 l=1,n
  a(icol,l)=a(icol,l)*pivinv
 16 continue
@@ -1160,7 +1165,7 @@ SUBROUTINE gaussj(a,n,np,b,m,mp)
  do 21 ll=1,n
  if(ll.ne.icol)then
  dum=a(ll,icol)
- a(ll,icol)=0.
+ a(ll,icol)=0.d0
  do 18 l=1,n
  a(ll,l)=a(ll,l)-a(icol,l)*dum
 18 continue