INRT_AugmentDisSources.F90

Go to the documentation of this file.

! *********************************************************************
! * Rocstar Simulation Suite                                          *
! * Copyright@2015, Illinois Rocstar LLC. All rights reserved.        *
! *                                                                   *
! * Illinois Rocstar LLC                                              *
! * Champaign, IL                                                     *
! * www.illinoisrocstar.com                                           *
! * sales@illinoisrocstar.com                                         *
! *                                                                   *
! * License: See LICENSE file in top level of distribution package or *
! * http://opensource.org/licenses/NCSA                               *
! *********************************************************************
! *********************************************************************
! * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,   *
! * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   *
! * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          *
! * NONINFRINGEMENT.  IN NO EVENT SHALL THE CONTRIBUTORS OR           *
! * COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER       *
! * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,   *
! * Arising FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE    *
! * USE OR OTHER DEALINGS WITH THE SOFTWARE.                          *
! *********************************************************************
!******************************************************************************
!
! Purpose: converts the transfers of primary quantities along Edges for each
!          particle into augmentations of RHS terms for all quantities, for
!          an interaction involving Lagrangian particles
!
! Description: none.
!
! Input: iInrt = index of interaction
!
! Output: augments region%levels(iLev)%...%rhs structures
!
! Notes:
!
!   The RHS structures use opposite sign as the input source structure
!
!   For efficiency, this routine requires Nodes to be stored in this order:
!     Mixture, Lagrangian particle, Eulerian particle, Internal
!
!   The energy corresponding to a gas mass is actually taken to be the
!   enthalpy because energy added in reactions are measured as enthalpies
!
!******************************************************************************
!
! $Id: INRT_AugmentDisSources.F90,v 1.5 2008/12/06 08:44:31 mtcampbe Exp $
!
! Copyright: (c) 2003 by the University of Illinois
!
!******************************************************************************

SUBROUTINE inrt_augmentdissources( region,iInrt )

  USE moddatatypes
  USE moddatastruct, ONLY : t_region
  USE modglobal,     ONLY : t_global
  USE modinteract,   ONLY : t_inrt_input,t_inrt_interact
  USE modmixture,    ONLY : t_mixt
  USE modpartlag,    ONLY : t_plag
#ifdef RFLO
#ifdef PEUL
  USE modparteul,    ONLY : t_peul
#endif
#endif
#ifdef RFLU
  USE modspecies,    ONLY : t_spec
#endif
  USE moderror
  USE modparameters
  USE inrt_modparameters
#ifdef PLAG
  USE plag_modparameters
#endif

  IMPLICIT NONE

! ... parameters
  TYPE(t_region), INTENT(INOUT), TARGET :: region
  INTEGER,        INTENT(IN)            :: iinrt

! ... loop variables
  INTEGER :: ipcls,iplag,ipeul,iedge

! ... local variables
  INTEGER, PARAMETER :: max_nodes = 101

  CHARACTER(CHRLEN)  :: rcsidentstring

  LOGICAL :: computeaux

  INTEGER :: npcls,nplag,npeul,nintl,ninputedges,nnodes,nedges
  INTEGER :: indmixt,indplag0,indpeul0,indintl
  INTEGER :: indplagvapor,indplagn,indpeuln
  INTEGER :: indcp,gasmodel,ic,inod,off1beg,ic1beg
#ifdef RFLO
  INTEGER :: ilev
#endif
  INTEGER :: tedge,inode(2),token(2)
  INTEGER, POINTER :: pcvplagmass(:), aiv(:,:)

  REAL(RFREAL) :: temp1,temp2,spht,massdot,hcapdot,enerdot
  REAL(RFREAL) :: kinedot1,kinedot2,thrmdot1,thrmdot2,enerdot1,enerdot2
  REAL(RFREAL) :: intlmass,intlener,intltemp,intlhcap,contfac
  REAL(RFREAL) :: sphtplag(max_nodes),sphtpeul(max_nodes)
  REAL(RFREAL), DIMENSION(3) :: velo1,velo2,intlvelo,intlmome
  REAL(RFREAL), DIMENSION(3) :: momedot,momedot1,momedot2
  REAL(RFREAL)               :: src(max_nodes,5)
  REAL(RFREAL), DIMENSION(:),   POINTER :: p1begmixttemp,pplagtemp
  REAL(RFREAL), DIMENSION(:,:), POINTER :: p1begmixtvelo,pplagvelo
  REAL(RFREAL), DIMENSION(:,:), POINTER :: pmixtrhs,pplagrhs,ppeulrhs
  REAL(RFREAL), DIMENSION(:,:), POINTER :: gv,arv,primary

  TYPE(t_inrt_input),    POINTER :: input
  TYPE(t_inrt_interact), POINTER :: inrt
  TYPE(t_mixt),          POINTER :: pmixt
#ifdef RFLO
#ifdef PLAG
  TYPE(t_plag),          POINTER :: pplag
#endif
#ifdef PEUL
  TYPE(t_peul),          POINTER :: ppeul
#endif
#endif
#ifdef RFLU
  TYPE(t_plag),          POINTER :: pplag
  TYPE(t_spec),          POINTER :: ppeul
#endif
  TYPE(t_global),        POINTER :: global

!******************************************************************************

  rcsidentstring = '$RCSfile: INRT_AugmentDisSources.F90,v $ $Revision: 1.5 $'

  global => region%global

  CALL registerfunction( global,'INRT_AugmentDisSources',&
  'INRT_AugmentDisSources.F90' )

#ifdef PLAG
! begin -----------------------------------------------------------------------

! Check if there are any particles

  npcls = 0

#ifdef RFLO
  ilev  = region%currLevel
  IF (global%plagUsed) npcls = region%levels(ilev)%plag%nPcls
#endif
#ifdef RFLU
  IF (global%plagUsed) npcls = region%plag%nPcls
#endif

  IF (npcls < 1) go to 999

#ifdef RFLU
! Check that have primitive state vector --------------------------------------

  IF ( region%mixt%cvState /= cv_mixt_state_duvwp ) THEN
    CALL errorstop(global,err_cv_state_invalid,__line__)
  END IF ! region%mixt%cvState
#endif

! initialize interaction constants and pointers

  input => region%inrtInput
  inrt  => input%inrts(iinrt)

  nplag = input%nPlag
  npeul = input%nPeul

  nintl       = inrt%nIntl
  ninputedges = inrt%nInputEdges

  indmixt  = input%indMixt
  indplag0 = input%indPlag0
  indpeul0 = input%indPeul0
  indintl  = input%indIntl

  indplagvapor = input%indPlagVapor

  indplagn = indplag0 + nplag
  indpeuln = indpeul0 + npeul

  nnodes = input%nNodes
  IF (nnodes > max_nodes .OR. nplag > max_nodes .OR. npeul > max_nodes) &
    CALL errorstop( global,err_exceeds_decl_mem,__line__ )

  nedges = inrt%nEdges

  computeaux = input%computeAux

! initialize data constants and pointers

#ifdef RFLO
  pmixt => region%levels(ilev)%mixt
  pplag => region%levels(ilev)%plag
#ifdef PEUL
  ppeul => region%levels(ilev)%peul
#endif
  indcp =  pmixt%indCp
#endif
#ifdef RFLU
  pmixt => region%mixt
  pplag => region%plag
  ppeul => region%spec

  indcp =  region%mixtInput%indCp
#endif

  pcvplagmass => pplag%cvPlagMass
  aiv => pplag%aiv
  arv => pplag%arv

! Constructing pointers from sections of arrays can be dangerous
! Extensive checks are employed to ensure correct behavior

! The string "1beg" is used to identify arrays beginning at 1 instead of ibc

#ifdef RFLO
  p1begmixtvelo => pmixt%dv(dv_mixt_uvel:dv_mixt_wvel,:)
  off1beg = lbound(p1begmixtvelo,2) - lbound(pmixt%dv,2)
#endif
#ifdef RFLU
  p1begmixtvelo => pmixt%cv(cv_mixt_xvel:cv_mixt_zvel,:)
  off1beg = lbound(p1begmixtvelo,2) - lbound(pmixt%cv,2)
#endif

  IF (ubound(p1begmixtvelo,1) /= 3) THEN
    CALL errorstop( global,err_inrt_indexrange,__line__ )
  ENDIF

  p1begmixttemp => pmixt%dv(dv_mixt_temp,:)

  IF (off1beg /= lbound(p1begmixttemp,1) - lbound(pmixt%dv,2)) THEN
    CALL errorstop( global,err_inrt_indexrange,__line__ )
  ENDIF

! The string "1beg" is not needed here:  Plag arrays already begin at 1

  pplagvelo => pplag%dv(dv_plag_uvel:dv_plag_wvel,:)

  IF (ubound(pplagvelo,1) /= 3) THEN
    CALL errorstop( global,err_inrt_indexrange,__line__ )
  ENDIF

  IF (lbound(pplagvelo,2) /= lbound(pplag%dv,2)) THEN
    CALL errorstop( global,err_inrt_indexrange,__line__ )
  ENDIF

  pplagtemp => pplag%dv(dv_plag_temp,:)

  IF (lbound(pplagtemp,1) /= lbound(pplag%dv,2)) THEN
    CALL errorstop( global,err_inrt_indexrange,__line__ )
  ENDIF

  pmixtrhs => pmixt%rhs
  pplagrhs => pplag%rhs
#ifdef RFLO
#ifdef PEUL
  ppeulrhs => ppeul%rhs
#endif
#endif
#ifdef RFLU
  ppeulrhs => ppeul%rhs
#endif
  gasmodel =  region%mixtInput%gasModel
  gv       => pmixt%gv

  IF (nplag > 0) sphtplag(1:nplag) = region%plagInput%spht(1:nplag)

  DO ipeul = 1,npeul
#ifdef RFLO
#ifdef PEUL
    sphtpeul(ipeul) = region%peulInput%ptypes(ipeul)%material%spht
#endif
#endif
! use spht for now:  should be spht at constant pressure (or volume?)
#ifdef RFLU
    sphtpeul(ipeul) = region%specInput%specType(ipeul)%pMaterial%spht
#endif
  END DO ! iPeul

  primary => pplag%inrtSources

  DO ipcls = 1,npcls

    ic = aiv(aiv_plag_icells,ipcls)
    ic1beg = ic + off1beg

! ----------------------
! - Compute source terms
! ----------------------

    src(:nnodes,:) = 0._rfreal
    intlhcap       = 0._rfreal

    DO iedge = 1, nedges

      tedge = inrt%edges(iedge)%tEdge
      inode = inrt%edges(iedge)%iNode
      token = inrt%edges(iedge)%token

      SELECT CASE (tedge)

! --- Mass Edge

      CASE (inrt_edge_mass)

        massdot = primary(iedge,ipcls)
        IF (massdot == 0._rfreal) cycle ! do not waste time with null transfer

! ----- BEGIN INLINE: makeSpecificHeat(iNode,spht) ------------------------BEGI

! ----- Nodes at both ends should have the same specific heat, so use either.
! ----- Restriction: cannot use an internal Node.

        IF (inode(1) == indintl) THEN
          inod = inode(2)
        ELSE
          inod = inode(1)
        END IF ! iNode(1)

        IF (inod == indmixt) THEN

          SELECT CASE (gasmodel)

          CASE (gas_model_tcperf)
            spht = gv(gv_mixt_cp,indcp*ic) ! specific heat at constant pressure

          CASE default
            CALL errorstop( global,err_reached_default,__line__ )

          END SELECT ! gasModel

        ELSE IF (inod <= indplagn) THEN
          spht = sphtplag(inod-indplag0)

        ELSE IF (inod == indintl .OR. inod == indplagvapor) THEN
          CALL errorstop( global,err_inrt_indexrange,__line__ )

        ELSE
          spht = sphtpeul(inod-indpeul0)

        END IF ! iNod

! ----- END   INLINE: makeSpecificHeat(iNode,spht) ------------------------ENDI

        hcapdot = spht * massdot

! ----- BEGIN INLINE: makeVelocity(iNode(1),velo1) ------------------------BEGI
! ----- BEGIN INLINE: makeTemperature(iNode(1),temp1) ---------------------BEGI

        inod = inode(1)

        IF (inod == indmixt) THEN
          velo1 = p1begmixtvelo(:,ic1beg)
          temp1 = p1begmixttemp(ic1beg)
        ELSE IF (inod <= indplagn) THEN
          velo1 = pplagvelo(:,ipcls)
          temp1 = pplagtemp(ipcls)
        ELSE IF (inod == indplagvapor) THEN
          CALL errorstop( global,err_inrt_indexrange,__line__ )
        ELSE IF (inod <= indpeuln) THEN
          velo1 = p1begmixtvelo(:,ic1beg) ! Sets Smoke velocity = Fluid velocity
          temp1 = p1begmixttemp(ic1beg)   ! Sets Smoke temp     = Fluid temp
        ELSE
          velo1 = intlvelo
          temp1 = intltemp
        END IF ! iNod

! ----- END   INLINE: makeVelocity(iNode(1),velo1) ------------------------ENDI
! ----- END   INLINE: makeTemperature(iNode(1),temp1) ---------------------ENDI

! ----- BEGIN INLINE: makeVelocity(iNode(2),velo2) ------------------------BEGI
! ----- BEGIN INLINE: makeTemperature(iNode(2),temp2) ---------------------BEGI

        inod = inode(2)

        IF (inod == indmixt) THEN
          velo2 = p1begmixtvelo(:,ic1beg)
          temp2 = p1begmixttemp(ic1beg)
        ELSE IF (inod <= indplagn) THEN
          velo2 = pplagvelo(:,ipcls)
          temp2 = pplagtemp(ipcls)
        ELSE IF (inod == indplagvapor) THEN
          CALL errorstop( global,err_inrt_indexrange,__line__ )
        ELSE IF (inod <= indpeuln) THEN
          velo2 = p1begmixtvelo(:,ic1beg) ! Sets Smoke velocity = Fluid velocity
          temp2 = p1begmixttemp(ic1beg)   ! Sets Smoke temp     = Fluid temp
        ELSE
          velo2 = 0._rfreal ! intlVelo not defined yet
          temp2 = 0._rfreal ! intlTemp not defined yet
! ------- If downwind Node is Internal, then augment its heat capacity
          intlhcap = intlhcap + hcapdot
        END IF ! iNod

! ----- END   INLINE: makeVelocity(iNode(2),velo2) ------------------------ENDI
! ----- END   INLINE: makeTemperature(iNode(2),temp2) ---------------------BEGI

        IF (computeaux) THEN

          momedot1 = massdot * velo1
          momedot2 = massdot * velo2

          kinedot1 = 0.5_rfreal*dot_product(velo1,momedot1)
          kinedot2 = 0.5_rfreal*dot_product(velo2,momedot2)

          thrmdot1 = hcapdot * temp1
          thrmdot2 = hcapdot * temp2

        ELSE

          momedot1 = 0._rfreal
          momedot2 = 0._rfreal

          kinedot1 = 0._rfreal
          kinedot2 = 0._rfreal

          thrmdot1 = 0._rfreal
          thrmdot2 = 0._rfreal

        ENDIF

! ----- Compute sources for upwind Node

        SELECT CASE(token(1))

        CASE(inrt_perm_pmass)
          src(inode(1),1)   = src(inode(1),1)   - massdot
          src(inode(1),2:4) = src(inode(1),2:4) - momedot1
          src(inode(1),5  ) = src(inode(1),5  ) - kinedot1 - thrmdot1

        CASE(inrt_perm_pmome)
          src(inode(1),1)   = src(inode(1),1)   - massdot
          src(inode(1),2:4) = src(inode(1),2:4) - momedot2
          src(inode(1),5  ) = src(inode(1),5  ) - kinedot2 - thrmdot1

        CASE(inrt_perm_pall)
          src(inode(1),1)   = src(inode(1),1)   - massdot
          src(inode(1),2:4) = src(inode(1),2:4) - momedot2
          src(inode(1),5  ) = src(inode(1),5  ) - kinedot2 - thrmdot2

        END SELECT ! token(1)

! ----- Compute sources for downwind Node

        SELECT CASE(token(2))

        CASE(inrt_perm_pmass)
          src(inode(2),1)   = src(inode(2),1)   + massdot
          src(inode(2),2:4) = src(inode(2),2:4) + momedot2
          src(inode(2),5  ) = src(inode(2),5  ) + kinedot2 + thrmdot2

        CASE(inrt_perm_pmome)
          src(inode(2),1)   = src(inode(2),1)   + massdot
          src(inode(2),2:4) = src(inode(2),2:4) + momedot1
          src(inode(2),5  ) = src(inode(2),5  ) + kinedot1 + thrmdot2

        CASE(inrt_perm_pall)
          src(inode(2),1)   = src(inode(2),1)   + massdot
          src(inode(2),2:4) = src(inode(2),2:4) + momedot1
          src(inode(2),5  ) = src(inode(2),5  ) + kinedot1 + thrmdot1

        END SELECT ! token(2)

! --- Momentum Edge

      CASE (inrt_edge_mome)

        momedot = primary(iedge:iedge+2,ipcls)

! ----- BEGIN INLINE: makeVelocity(iNode(1),velo1) ------------------------BEGI

        inod = inode(1)

        IF (inod == indmixt) THEN
          velo1 = p1begmixtvelo(:,ic1beg)
        ELSE IF (inod <= indplagn) THEN
          velo1 = pplagvelo(:,ipcls)
        ELSE IF (inod == indplagvapor) THEN
          CALL errorstop( global,err_inrt_indexrange,__line__ )
        ELSE IF (inod <= indpeuln) THEN
          velo1 = p1begmixtvelo(:,ic1beg) ! Sets Smoke velocity = Fluid velocity
        ELSE
          velo1 = intlvelo
        END IF ! iNod

! ----- END   INLINE: makeVelocity(iNode(1),velo1) ------------------------ENDI

! ----- BEGIN INLINE: makeVelocity(iNode(2),velo2) ------------------------BEGI

        inod = inode(2)

        IF (inod == indmixt) THEN
          velo2 = p1begmixtvelo(:,ic1beg)
        ELSE IF (inod <= indplagn) THEN
          velo2 = pplagvelo(:,ipcls)
        ELSE IF (inod == indplagvapor) THEN
          CALL errorstop( global,err_inrt_indexrange,__line__ )
        ELSE IF (inod <= indpeuln) THEN
          velo2 = p1begmixtvelo(:,ic1beg) ! Sets Smoke velocity = Fluid velocity
        ELSE
          velo2 = 0._rfreal ! intlVelo not defined yet
        END IF ! iNod

! ----- END   INLINE: makeVelocity(iNode(2),velo2) ------------------------ENDI

        IF (computeaux) THEN

          enerdot1 = dot_product(momedot,velo1)
          enerdot2 = dot_product(momedot,velo2)

        ELSE

          enerdot1 = 0._rfreal
          enerdot2 = 0._rfreal

        ENDIF

! ----- Compute sources for upwind Node

        SELECT CASE(token(1))

        CASE(inrt_perm_pmome)
          src(inode(1),2:4) = src(inode(1),2:4) - momedot
          src(inode(1),5  ) = src(inode(1),5  ) - enerdot1

        CASE(inrt_perm_pall)
          src(inode(1),2:4) = src(inode(1),2:4) - momedot
          src(inode(1),5  ) = src(inode(1),5  ) - enerdot2

        END SELECT ! token(1)

! ----- Compute sources for downwind Node

        SELECT CASE(token(2))

        CASE(inrt_perm_pmome)
          src(inode(2),2:4) = src(inode(2),2:4) + momedot
          src(inode(2),5)   = src(inode(2),5)   + enerdot2

        CASE(inrt_perm_pall)
          src(inode(2),2:4) = src(inode(2),2:4) + momedot
          src(inode(2),5)   = src(inode(2),5)   + enerdot1

        END SELECT ! token(2)

! --- Dummy momentum Edge

      CASE (inrt_edge_mome_dum)
        cycle

! --- Energy Edge

      CASE (inrt_edge_ener)

        enerdot = primary(iedge,ipcls)

        IF (token(1) == inrt_perm_pall) &
          src(inode(1),5) = src(inode(1),5) - enerdot

        IF (token(2) == inrt_perm_pall) &
          src(inode(2),5) = src(inode(2),5) + enerdot

! --- Ghost mass Edge

      CASE (inrt_edge_mass_gho)

! ----- Downwind node must have Block Token.  Check value of upwind node

        IF (token(1) >= inrt_perm_pmass) THEN

          massdot = primary(iedge,ipcls)
          src(inode(1),1) = src(inode(1),1) - massdot

        END IF ! token(1)

      CASE default
        CALL errorstop( global,err_reached_default,__line__ )

      END SELECT ! tEdge

! --- Compute velocity and temperature at internal Node if necessary

      IF (nintl > 0) THEN
        IF (iedge == ninputedges) THEN

          intlmass = src(indintl,1)
          intlmome = src(indintl,2:4)
          intlener = src(indintl,5)

          IF (intlmass > 0._rfreal) THEN
            intlvelo = intlmome / intlmass
          ELSE
            intlvelo = 0._rfreal
          END IF ! intlMass

          IF (intlhcap > 0._rfreal) THEN
            intltemp = (intlener - 0.5_rfreal*dot_product(intlmome, &
                        intlvelo)) / intlhcap
          ELSE
            intltemp = 0._rfreal
          END IF ! intlHcap

        END IF ! iEdge
      END IF ! nIntl

    END DO ! iEdge

! ------------------
! - Augment RHS data
! ------------------

! - Continuum conversion factor: converts from single particle value to
! - superparticle value

! - Note: does not incorporate volume because rhs values are not per volume
! - for either particles or continuua, in contrast to cv values, which are
! - not per volume for particles, but are per volume for continuua.

    contfac = arv(arv_plag_spload,ipcls)

! - Augment Gas Sources

    pmixtrhs(  cv_mixt_dens             ,ic) = &
      pmixtrhs(cv_mixt_dens             ,ic) - contfac*src(indmixt,1  )

    pmixtrhs(  cv_mixt_xmom:cv_mixt_zmom,ic) = &
      pmixtrhs(cv_mixt_xmom:cv_mixt_zmom,ic) - contfac*src(indmixt,2:4)

    pmixtrhs(  cv_mixt_ener             ,ic) = &
      pmixtrhs(cv_mixt_ener             ,ic) - contfac*src(indmixt,5  )

! - Augment Lagrangian Particle Sources

    DO iplag = 1, nplag

      pplagrhs(  pcvplagmass(iplag)       ,ipcls) = &
        pplagrhs(pcvplagmass(iplag)       ,ipcls) - src(indplag0+iplag,1)

      pplagrhs(  cv_plag_xmom:cv_plag_zmom,ipcls) = &
        pplagrhs(cv_plag_xmom:cv_plag_zmom,ipcls) - src(indplag0+iplag,2:4)

      pplagrhs(  cv_plag_ener             ,ipcls) = &
        pplagrhs(cv_plag_ener             ,ipcls) - src(indplag0+iplag,5)

    END DO ! iPlag

! - Augment Lagrangian Particle Vapor Energy

    pplagrhs(  cv_plag_enervapor,ipcls) = &
      pplagrhs(cv_plag_enervapor,ipcls) - src(indplagvapor,5)

! - Augment Smoke Sources
#ifdef RFLO
#ifdef PEUL
    DO ipeul = 1, npeul
      ppeulrhs(ipeul,ic) = ppeulrhs(ipeul,ic) - contfac * src(indpeul0+ipeul,1)
    END DO ! iPeul
#endif
#endif
#ifdef RFLU
    DO ipeul = 1, npeul
      ppeulrhs(ipeul,ic) = ppeulrhs(ipeul,ic) - contfac * src(indpeul0+ipeul,1)
    END DO ! iPeul
#endif
  END DO ! iPcls

! finalize --------------------------------------------------------------------

999  CONTINUE
#endif
  CALL deregisterfunction( global )

END SUBROUTINE inrt_augmentdissources

!******************************************************************************
!
! RCS Revision history:
!
! $Log: INRT_AugmentDisSources.F90,v $
! Revision 1.5  2008/12/06 08:44:31  mtcampbe
! Updated license.
!
! Revision 1.4  2008/11/19 22:17:44  mtcampbe
! Added Illinois Open Source License/Copyright
!
! Revision 1.3  2006/02/15 20:18:30  wasistho
! put peul within ifdef
!
! Revision 1.2  2005/10/31 21:09:37  haselbac
! Changed specModel and SPEC_MODEL_NONE
!
! Revision 1.1  2004/12/01 21:56:10  fnajjar
! Initial revision after changing case
!
! Revision 1.15  2004/07/28 15:42:12  jferry
! deleted defunct constructs: useDetangle, useSmokeDrag, useSmokeHeatTransfer
!
! Revision 1.14  2004/07/23 22:43:16  jferry
! Integrated rocspecies into rocinteract
!
! Revision 1.13  2004/03/25 21:14:53  jferry
! fixed pointer offset bug
!
! Revision 1.12  2004/03/05 22:09:03  jferry
! created global variables for peul, plag, and inrt use
!
! Revision 1.11  2004/03/02 21:47:29  jferry
! Added After Update interactions
!
! Revision 1.10  2004/01/31 03:59:22  haselbac
! Initial integration for Rocflu and Rocpart
!
! Revision 1.9  2003/09/19 20:35:26  jferry
! Implemented oxidizer species for burning interaction
!
! Revision 1.8  2003/05/08 17:17:14  jferry
! changed energy associated with mass to enthalpy
!
! Revision 1.7  2003/05/07 15:13:10  jferry
! Rearranged for efficiency
!
! Revision 1.6  2003/04/09 15:02:39  jferry
! removed erroneous volume normalization for continuum rhs
!
! Revision 1.5  2003/04/03 21:10:17  jferry
! implemented additional safety checks for rocinteract
!
! Revision 1.4  2003/04/02 22:32:03  jferry
! codified Activeness and Permission structures for rocinteract
!
! Revision 1.3  2003/03/24 23:30:52  jferry
! overhauled rocinteract to allow interaction design to use user input
!
! Revision 1.2  2003/03/11 16:05:54  jferry
! Created data type for material properties
!
! Revision 1.1  2003/03/04 22:12:35  jferry
! Initial import of Rocinteract
!
!******************************************************************************