trlimit.f.html
Source file:   trlimit.f
Directory:   /Users/rjl/git/rjleveque/clawpack-4.6.3/book/chap9/acoustics/layered
Converted:   Mon Jan 21 2013 at 20:15:48   using clawcode2html
This documentation file will not reflect any later changes in the source file.

 
c
c
c ===================================================================
      subroutine trlimit(maxmx,meqn,mwaves,mbc,mx,aux,
     &                   wave,s,mthlim)
c ===================================================================
c
c     # Transmission-based limiter for acoustics equations
c     #    See Fogarty and LeVeque paper for a description.
c
c     # Use in place of limiter.f
c
c     # Assumes 
c     #    aux(i,1) = sound speed c in i'th cell
c     #    aux(i,2) = impedance Z in i'th cell
c     #    wave(i,*,1) is the left-going wave 
c     #      = alf_i^1 [Z_{i-1} ; 1]
c     #    wave(i,*,2) is the right-going wave 
c     #      = alf_i^2 [Z_i ; 1]
c     # 
c
c     --------------------------------------------------------------------
c
      implicit double precision (a-h,o-z)
      dimension  aux(1-mbc:maxmx+mbc, *)
      dimension    s(1-mbc:maxmx+mbc, mwaves)
      dimension wave(1-mbc:maxmx+mbc, meqn, mwaves)
      dimension mthlim(mwaves)
c
c
c     # transmission-based limiter:
c     # use the fact that second component of wave is strength alpha since
c     # second component of eigenvector is 1.
c
       do 100 i=0,mx+1
c          # 1-wave at this cell and neighbor:
	   alf1i = wave(i,2,1)
	   if (alf1i .ne. 0.d0) then
	      alf1ip = wave(i+1,2,1)
c             # transmitted part of neighboring 1-wave:
	      alf1ipt = (2.d0*aux(i,1)/(aux(i-1,1)+aux(i,1))) * alf1ip
	      wlimitr = philim(alf1i, alf1ipt, mthlim(1))
	      do m=1,meqn
	         wave(i,m,1) = wlimitr * wave(i,m,1)
	         enddo
	      endif

c          # 2-wave at this cell and neighbor:
	   alf2i = wave(i,2,2)
	   if (alf2i .ne. 0.d0) then
	      alf2im = wave(i-1,2,2)
c             # transmitted part of neighboring 2-wave:
  	      alf2imt = (2.d0*aux(i-1,1)/(aux(i-1,1)+aux(i,1))) * alf2im
  	      wlimitr = philim(alf2i, alf2imt, mthlim(2))
	      do m=1,meqn
	         wave(i,m,2) = wlimitr * wave(i,m,2)
	         enddo
	      endif
  100      continue
c
      return
      end