qinit.f.html
Source file:   qinit.f
Directory:   /var/www/html/clawpack/links/shockvacuum10/shockvacuum
Converted:   Mon Aug 16 2010 at 14:13:31   using clawcode2html
This documentation file will not reflect any later changes in the source file.

 
c
c
c =========================================================
       subroutine qinit(maxmx,meqn,mbc,mx,xlower,dx,q,maux,aux)
c =========================================================
c
c     # Set initial conditions for q as a perturbation of the 
c     # steady state solution evaluated pointwise
c
c
      implicit double precision (a-h,o-z)
      dimension q(1-mbc:maxmx+mbc, meqn)
      dimension aux(1-mbc:maxmx+mbc, *)
      double precision amp, center, width, goverc0, grav
      double precision npoly, nstar
      common /ceos/ grav, gamma, gamma1, sK, npoly, nstar, Crho, ieos
c
c
c     # Parameters for perturbation:
      amp = 10.d0
      center = -11.d0
      width = 0.3d0

      write(6,*) '==========='
      write(6,*) 'ieos = ',ieos
      write(6,*) 'grav = ',grav
      write(6,*) 'npoly = ',npoly
      write(6,*) 'gamma = ',gamma
      write(6,*) 'sK = ',sK
      write(6,*) 'nstar = ',nstar
      write(6,*) 'Crho = ',Crho
      write(6,*) '==========='

      do 150 i=mx,1,-1
	     xc = xlower + (i-0.5d0)*dx

c            old version used this:
c            q(i,1) = (-grav*xc / (npoly*gamma*sK))**npoly

             q(i,1) = Crho*(-xc)**nstar

c            # add perturbations, but only where exp doesn't underflow:

             qexp = -(xc-center)**2/width**2
             if (qexp.gt.-50.d0) then
                q(i,1) = q(i,1) + amp*exp(qexp)
                endif

             q(i,2) = 0.d0
             qexp = -(xc-center-1.d0)**2/width**2
             if (qexp.gt.-50.d0) then
                q(i,2) = q(i,2) + 3.d0*amp*exp(qexp)
                endif

c            # old version:
c            q(i,2) = 3*amp*exp(-(xc-center-1.0)**2/width**2)

c            # isentropic (old version used this):
c            pressure = sK * q(i,1)**gamma 
             
c            # general stratified:
             pressure = grav/(nstar+1.d0)*Crho * (-xc)**(nstar+1.d0)
             q(i,3) = pressure / gamma1 +
     &                   0.5d0*q(i,2)**2/q(i,1)
c
  150        continue

c
      return
      end