|
bc1.f.html |
|
|
Source file: bc1.f
|
|
Directory: /Users/rjl/clawpack_src/clawpack_master/geoclaw/examples/1d_classic/bouss_wavetank_matsuyama
|
|
Converted: Mon Feb 19 2024 at 14:30:59
using clawcode2html
|
|
This documentation file will
not reflect any later changes in the source file.
|
c
c
c =================================================================
subroutine bc1(meqn,mbc,mx,xlower,dx,q,maux,aux,t,dt,mthbc)
c =================================================================
c
c # Standard boundary condition choices for claw2
c
c # Modified for 1d GeoClaw to extend topo aux(1,:) to ghost cells.
c
c # At each boundary k = 1 (left), 2 (right):
c # mthbc(k) = 0 for user-supplied BC's (must be inserted!)
c # = 1 for zero-order extrapolation
c # = 2 for periodic boundary coniditions
c # = 3 for solid walls, assuming this can be implemented
c # by reflecting the data about the boundary and then
c # negating the 2'nd component of q.
c ------------------------------------------------
c
c # Extend the data from the computational region
c # i = 1, 2, ..., mx2
c # to the virtual cells outside the region, with
c # i = 1-ibc and i = mx+ibc for ibc=1,...,mbc
c
use geoclaw_module, only: grav
implicit double precision (a-h,o-z)
dimension q(meqn,1-mbc:mx+mbc)
dimension aux(maux,1-mbc:mx+mbc)
dimension mthbc(2)
pi = 4.d0*atan(1.d0)
h0 = 4.d0
ampl_eta = 0.03
ampl_u = ampl_eta * sqrt(grav/h0)
!write(6,*) '+++ ampl_eta, ampl_u: ',ampl_eta, ampl_u
!ampl = 0.047 #0.15d0
tperiod = 20.d0
c
c
c-------------------------------------------------------
c # left boundary:
c-------------------------------------------------------
do ibc=1,mbc
aux(1,1-ibc)=aux(1,1)
enddo
go to (100,110,120,130) mthbc(1)+1
c
100 continue
c # wave maker
if (t < tperiod) then
s = ampl_u * sin(2.d0*pi*t/tperiod)
else
s = 0.d0
endif
do ibc=1,mbc
aux(1,1-ibc) = aux(1,ibc)
do m=1,meqn
q(m,1-ibc) = q(m,ibc)
enddo
u = q(2,ibc)/q(1,ibc)
q(2,1-ibc) = (2.d0*s - u) * q(1,1-ibc)
enddo
go to 199
c
110 continue
c # zero-order extrapolation:
do ibc=1,mbc
do m=1,meqn
q(m,1-ibc) = q(m,1)
end do
end do
go to 199
120 continue
c # periodic:
do ibc=1,mbc
do m=1,meqn
q(m,1-ibc) = q(m,mx+1-ibc)
end do
end do
go to 199
130 continue
c # solid wall (assumes 2'nd component is velocity or momentum in x):
do ibc=1,mbc
do m=1,meqn
q(m,1-ibc) = q(m,ibc)
end do
c # negate the normal velocity:
q(2,1-ibc) = -q(2,ibc)
end do
go to 199
199 continue
c
c-------------------------------------------------------
c # right boundary:
c-------------------------------------------------------
do ibc=1,mbc
aux(1,mx+ibc)=aux(1,mx)
enddo
go to (200,210,220,230) mthbc(2)+1
c
200 continue
c # user-specified boundary conditions go here in place of error output
write(6,*) '*** ERROR *** mthbc(2)=0 and no BCs specified in bc2'
stop
go to 299
210 continue
c # zero-order extrapolation:
do ibc=1,mbc
do m=1,meqn
q(m,mx+ibc) = q(m,mx)
end do
end do
go to 299
220 continue
c # periodic:
do ibc=1,mbc
do m=1,meqn
q(m,mx+ibc) = q(m,ibc)
end do
end do
go to 299
230 continue
c # solid wall (assumes 2'nd component is velocity or momentum in x):
do ibc=1,mbc
do m=1,meqn
q(m,mx+ibc) = q(m,mx+1-ibc)
end do
c # negate the normal velocity:
q(2,mx+ibc) = -q(2,mx+1-ibc)
end do
go to 299
299 continue
c
return
end