Skip to content

Commit

Permalink
fix indentation in rpt2_geoclaw.f
Browse files Browse the repository at this point in the history
  • Loading branch information
@rjleveque
rjleveque committed Oct 1, 2020
1 parent cd39c82 commit d36f4bd
Showing 1 changed file with 107 additions and 108 deletions.
215 changes: 107 additions & 108 deletions src/rpt2_geoclaw.f
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@ subroutine rpt2(ixy,imp,maxm,meqn,mwaves,maux,mbc,mx,
! is split equally between bmasdq and bpasdq to improve symmetry.
!

!-----------------------last modified June 2020 ----------------------
!-----------------------last modified October 2020 ----------------------

integer ixy,maxm,meqn,maux,mwaves,mbc,mx,imp

Expand All @@ -45,11 +45,11 @@ subroutine rpt2(ixy,imp,maxm,meqn,mwaves,maux,mbc,mx,
abs_tol=tol

if (ixy.eq.1) then
mu = 2
mv = 3
mu = 2
mv = 3
else
mu = 3
mv = 2
mu = 3
mv = 2
endif

do i=2-mbc,mx+mbc
Expand All @@ -61,99 +61,102 @@ subroutine rpt2(ixy,imp,maxm,meqn,mwaves,maux,mbc,mx,
hvl=qr(mv,i-1)
hvr=ql(mv,i)

!===========determine velocity from momentum===========================
if (hl.lt.abs_tol) then
hl=0.d0
ul=0.d0
vl=0.d0
else
ul=hul/hl
vl=hvl/hl
endif

if (hr.lt.abs_tol) then
hr=0.d0
ur=0.d0
vr=0.d0
else
ur=hur/hr
vr=hvr/hr
endif

do mw=1,mwaves
s(mw)=0.d0
beta(mw)=0.d0
do m=1,meqn
r(m,mw)=0.d0
enddo
enddo
dxdcp = 1.d0
dxdcm = 1.d0

if (hl <= tol .and. hr <= tol) go to 90

* !check and see if cell that transverse waves are going in is high and dry
if (imp.eq.1) then
eta = qr(1,i-1) + aux2(1,i-1)
topo1 = aux1(1,i-1)
topo3 = aux3(1,i-1)
else
eta = ql(1,i) + aux2(1,i)
topo1 = aux1(1,i)
topo3 = aux3(1,i)
endif
if (eta.lt.max(topo1,topo3)) go to 90

if (coordinate_system.eq.2) then
if (ixy.eq.2) then
dxdcp=(earth_radius*deg2rad)
dxdcm = dxdcp
! determine velocity from momentum
if (hl.lt.abs_tol) then
hl=0.d0
ul=0.d0
vl=0.d0
else
if (imp.eq.1) then
dxdcp = earth_radius*cos(aux3(3,i-1))*deg2rad
dxdcm = earth_radius*cos(aux1(3,i-1))*deg2rad
else
dxdcp = earth_radius*cos(aux3(3,i))*deg2rad
dxdcm = earth_radius*cos(aux1(3,i))*deg2rad
endif
ul=hul/hl
vl=hvl/hl
endif
endif

if (hr.lt.abs_tol) then
hr=0.d0
ur=0.d0
vr=0.d0
else
ur=hur/hr
vr=hvr/hr
endif

do mw=1,mwaves
s(mw)=0.d0
beta(mw)=0.d0
do m=1,meqn
r(m,mw)=0.d0
enddo
enddo

c=====Determine some speeds necessary for the Jacobian=================
dxdcp = 1.d0
dxdcm = 1.d0


! In v5.7.x and prior versions,
! we used left right states to define Roe averages,
! which is consistent with those used in rpn2.
! But now we are computing upgoing, downgoing waves either in
! cell on left (if imp==1) or on right (if imp==2) so we
! should possibly use q values in cells above/below,
! but these aren't available (only aux values).
! At any rate, there is no clear justification for using cells
! on the other side of the normal-solve interface.
! Modified to use left or right state alone in defining
! Jacobian, based on imp:
if (imp == 1) then
! asdq is leftgoing, use q from cell i-1:
if (hl <= tol) go to 90
s(1) = vl-dsqrt(g*hl)
s(2) = vl
s(3) = vl+dsqrt(g*hl)
ulr = ul
hlr = hl
if (hl <= tol .and. hr <= tol) go to 90

c ! check to see if cell that transverse waves are going in is dry

if (imp.eq.1) then
eta = qr(1,i-1) + aux2(1,i-1)
topo1 = aux1(1,i-1)
topo3 = aux3(1,i-1)
else
eta = ql(1,i) + aux2(1,i)
topo1 = aux1(1,i)
topo3 = aux3(1,i)
endif

if (eta.lt.max(topo1,topo3)) go to 90

if (coordinate_system.eq.2) then
if (ixy.eq.2) then
dxdcp=(earth_radius*deg2rad)
dxdcm = dxdcp
else
! asdq is rightgoing, use q from cell i:
if (hr <= tol) go to 90
s(1) = vr-dsqrt(g*hr)
s(2) = vr
s(3) = vr+dsqrt(g*hr)
ulr = ur
hlr = hr
if (imp.eq.1) then
dxdcp = earth_radius*cos(aux3(3,i-1))*deg2rad
dxdcm = earth_radius*cos(aux1(3,i-1))*deg2rad
else
dxdcp = earth_radius*cos(aux3(3,i))*deg2rad
dxdcm = earth_radius*cos(aux1(3,i))*deg2rad
endif
endif
endif

c Determine some speeds necessary for the Jacobian

! In v5.7.x and prior versions,
! we used left right states to define Roe averages,
! which is consistent with those used in rpn2.
! But now we are computing upgoing, downgoing waves either in
! cell on left (if imp==1) or on right (if imp==2) so we
! should possibly use q values in cells above/below,
! but these aren't available (only aux values).
! At any rate, there is no clear justification for using cells
! on the other side of the normal-solve interface.
! Modified to use left or right state alone in defining
! Jacobian, based on imp:
if (imp == 1) then
! asdq is leftgoing, use q from cell i-1:
if (hl <= tol) go to 90
s(1) = vl-dsqrt(g*hl)
s(2) = vl
s(3) = vl+dsqrt(g*hl)
ulr = ul
hlr = hl
else
! asdq is rightgoing, use q from cell i:
if (hr <= tol) go to 90
s(1) = vr-dsqrt(g*hr)
s(2) = vr
s(3) = vr+dsqrt(g*hr)
ulr = ur
hlr = hr
endif
c Determine asdq decomposition (beta)
c=======================Determine asdq decomposition (beta)============
delf1=asdq(1,i)
delf2=asdq(mu,i)
delf3=asdq(mv, i)
Expand All @@ -162,9 +165,8 @@ subroutine rpt2(ixy,imp,maxm,meqn,mwaves,maux,mbc,mx,
beta(1) = (s(3)*delf1 - delf3) / (s(3) - s(1))
beta(2) = -ulr*delf1 + delf2
beta(3) = (delf3 - s(1)*delf1) / (s(3) - s(1))
c======================End =================================================
c=====================Set-up eigenvectors===================================
c Set-up eigenvectors
r(1,1) = 1.d0
r(2,1) = ulr ! fixed bug, ulr not s(2)=vlr
r(3,1) = s(1)
Expand All @@ -176,18 +178,18 @@ subroutine rpt2(ixy,imp,maxm,meqn,mwaves,maux,mbc,mx,
r(1,3) = 1.d0
r(2,3) = ulr ! fixed bug, ulr not s(2)=vlr
r(3,3) = s(3)
c============================================================================
90 continue
c============= compute fluctuations==========================================
90 continue
! compute fluctuations
! initialize all components to 0:
bmasdq(:,i) = 0.d0
bpasdq(:,i) = 0.d0
do mw=1,3
if ((abs(s(mw)) > 0.d0) .and.
& (abs(s(mw)) < 0.001d0*dsqrt(g*hlr))) then
do mw=1,3
if ((abs(s(mw)) > 0.d0) .and.
& (abs(s(mw)) < 0.001d0*dsqrt(g*hlr))) then
! split correction symmetrically if nearly zero
! Note wave drops out if s(mw)==0 exactly, so no need to split
bmasdq(1,i) =bmasdq(1,i) +
Expand All @@ -202,21 +204,18 @@ subroutine rpt2(ixy,imp,maxm,meqn,mwaves,maux,mbc,mx,
& 0.5d0*dxdcp*s(mw)*beta(mw)*r(2,mw)
bpasdq(mv,i)=bpasdq(mv,i)+
& 0.5d0*dxdcp*s(mw)*beta(mw)*r(3,mw)
elseif (s(mw).lt.0.d0) then
elseif (s(mw).lt.0.d0) then
bmasdq(1,i) =bmasdq(1,i) + dxdcm*s(mw)*beta(mw)*r(1,mw)
bmasdq(mu,i)=bmasdq(mu,i)+ dxdcm*s(mw)*beta(mw)*r(2,mw)
bmasdq(mv,i)=bmasdq(mv,i)+ dxdcm*s(mw)*beta(mw)*r(3,mw)
elseif (s(mw).gt.0.d0) then
elseif (s(mw).gt.0.d0) then
bpasdq(1,i) =bpasdq(1,i) + dxdcp*s(mw)*beta(mw)*r(1,mw)
bpasdq(mu,i)=bpasdq(mu,i)+ dxdcp*s(mw)*beta(mw)*r(2,mw)
bpasdq(mv,i)=bpasdq(mv,i)+ dxdcp*s(mw)*beta(mw)*r(3,mw)
endif
enddo
c========================================================================
enddo
c
endif
enddo ! loop on mw
c
enddo ! loop on i
return
end

0 comments on commit d36f4bd

Please sign in to comment.