c******************************************************************************* c ... trpack.f c ... package of subroutines that access the internal routines of the c ... chemkin library to get various thermodynamic and transport property c ... values c******************************************************************************* c ... get species names c ... c ... prerequisite calls : trinit c ... c ... inputs: c ... nsp : number of species c ... c ... outputs: c ... spname(i),i=1,...,nsp : names of species, max lensym char. each c ... c------------------------------------------------------------------------------- subroutine getnames(nsp,spname) implicit double precision (a-h,o-z) include 'parameter.com' character*8 spname(nsp) logical lerr include 'comwk.com' include 'comfl.com' c------------------------------------------------------------------------------- c ... get character strings for names of the species in array spname() call cksyms (cwork(icc), lout, cwork(iks), lerr) do i=1,nsp spname(i)=cwork(iks+i-1) end do if(lerr)then write(*,*)' getnames: *** error return from cksyms' stop endif return end c******************************************************************************* c ... get multicomponent diffusion coefficients, c ... c ... prerequisites: trinit c ... c ... inputs: c ... nsp = number of species c ... temp = temperature in degrees K c ... pressure = pressure in SI units, Pascals (N/m^2) c ... x(i) = mole fractions of species i=1,...,nsp c ... c ... outputs: c ... dkj(i,j) = multicomponent diffusion coefficients, SI units, (m^2/s) c ... of species i in species j, (i,j=1,...,nsp) c ... c ... notes: c ... x() can have any dimension >= nsp in the calling program c ... if dkj() is a 2D array in the calling program (CPO) then its first c ... dimension in the CPO has to be exactly nsp c ... to avoid this restriction, use dkj() as a 1D vector in the c ... CPO and dimension it with size >= nsp**2. Then, the correspondence c ... between dkj_here(i,j) here and dkj_CPO(k) in the CPO is : c ... dkj_CPO( k ) = dkj_here(i,j) for k=(j-1)*nsp+i c ... c ... do NOT change the sizes of arrays and other parameters defined c ... in the parameter statement unless you do the same everywhere c ... else (including inside certain chemkin routines) c ... common blocks ckwblk and flflfl tie this routine to various c ... work space arrays in chemkin c------------------------------------------------------------------------------- subroutine getmcdij(temp,pressure,nsp,x,dkj) implicit double precision (a-h,o-z) include 'parameter.com' dimension x(nsp),dkj(nsp,nsp) include 'comwk.com' include 'comfl.com' c------------------------------------------------------------------------------- c ... find pressure in cgs units (dynes/cm^2) c ... 1 N = 10^5 dynes ---> 1 Pa = 10 dynes/cm^2 pcgs = pressure*10.0d0 c------------------------------------------------------------------------------- c ... call the chemkin routine for multicomponent diffusion coeff. call mcmdif(pcgs,temp,x,nsp,iwork(imcw),rwork(nmcw),dkj(1,1)) c------------------------------------------------------------------------------- c ... dkj(i,j) are outputs, in cgs units : cm^2/s, convert to SI (m^2/s) convert = 1.d-4 do i=1,nsp do j=1,nsp dkj(i,j)=convert*dkj(i,j) end do end do c------------------------------------------------------------------------------- return end c******************************************************************************* c ... get multicomponent diffusion coefficients, from Chemkin subroutine getmadc(temp,pressure,nsp,x,dkj) implicit double precision (a-h,o-z) include 'parameter.com' dimension x(nsp),dkj(nsp) include 'comwk.com' include 'comfl.com' c------------------------------------------------------------------------------- c ... find pressure in cgs units (dynes/cm^2) c ... 1 N = 10^5 dynes ---> 1 Pa = 10 dynes/cm^2 pcgs = pressure*10.0d0 c------------------------------------------------------------------------------- c ... call the chemkin routine for multicomponent diffusion coeff. call mcadif(pcgs,temp,x,rwork(nmcw),dkj(1)) c------------------------------------------------------------------------------- c ... dkj(i) are outputs, in cgs units : cm^2/s, convert to SI (m^2/s) convert = 1.d-4 do i=1,nsp dkj(i)=convert*dkj(i) end do c------------------------------------------------------------------------------- return end c******************************************************************************* c ... get multicomponent thermal diffusion coefficients c ... and the special "multicomponent thermal conductivity" c ... this is not the true thermal conductivity of the mixture c ... c ... prerequisites: trinit c ... c ... inputs: c ... nsp = number of species c ... temp = temperature in degrees K c ... pressure = pressure in SI units, Pascals (N/m^2) c ... x(i) = mole fractions of species i=1,...,nsp c ... c ... outputs: c ... tdr(i) = "multicomponent" thermal diffusion coefficients, c ... SI units, (kg/m.s) c ... of species i (i=1,...,nsp) c ... condpr = special "multicomponent thermal conductivity" c ... this is the term (lambda_prime) in HCB, used c ... in the conduction term in the multicomponent c ... diffusion coeff. formulation c ... SI units, (W/m.K) c ... c ... notes: c ... x(),tdr() can have any dimension >= nsp in the calling program c ... c ... do NOT change the sizes of arrays and other parameters defined c ... in the parameter statement unless you do the same everywhere c ... else (including inside certain chemkin routines) c ... common blocks ckwblk and flflfl tie this routine to various c ... work space arrays in chemkin c------------------------------------------------------------------------------- subroutine getmcdti(temp,pressure,nsp,x,tdr,condpr) implicit double precision (a-h,o-z) include 'parameter.com' dimension x(nsp),tdr(nsp) include 'comwk.com' include 'comfl.com' c------------------------------------------------------------------------------- c ... find pressure in cgs units (dynes/cm^2) c ... 1 N = 10^5 dynes ---> 1 Pa = 10 dynes/cm^2 pcgs = pressure*10.0d0 c------------------------------------------------------------------------------- c ... get multicomponent thermal diffusion coefficients c ... and the special "multicomponent thermal conductivity" call mcmcdt(pcgs,temp,x,iwork(imcw),rwork(nmcw),iwork(ickw), 1 rwork(nckw),tdr(1),condpr) c------------------------------------------------------------------------------- c ... tdr(i) are in cgs units: gm/cm.s, convert to SI (kg/m.s) convert = 0.1d0 do i=1,nsp tdr(i)=convert*tdr(i) end do c------------------------------------------------------------------------------- c ... condpr is in cgs units: erg/cm.K.s, convert to SI (J/m.K.s) condpr = 1.d-5*condpr c------------------------------------------------------------------------------- return end c******************************************************************************* c ... get molecular weights of species c ... c ... prerequisites: trinit c ... c ... inputs: c ... nsp = number of species c ... c ... outputs: c ... wt() = molecular weights of species i=1,...,nsp c ... c ... notes: c ... wt() can have any dimension >= nsp in the calling program (Kg/K-mol or gr/mol) c ... c ... do NOT change the sizes of arrays and other parameters defined c ... in the parameter statement unless you do the same everywhere c ... else (including inside certain chemkin routines) c ... common blocks ckwblk and flflfl tie this routine to various c ... work space arrays in chemkin c------------------------------------------------------------------------------- subroutine getwt(nsp,wref,wt) implicit double precision (a-h,o-z) include 'parameter.com' dimension wt(nsp) include 'comwk.com' include 'comfl.com' c------------------------------------------------------------------------------- call ckwt (iwork(ickw), rwork(nckw), wt) #ifdef NON_DIM do i=1,nsp wt(i) = wt(i)/wref enddo #endif c------------------------------------------------------------------------------- return end c******************************************************************************* c ... get mass fractions of species in mixture c ... c ... prerequisites: trinit c ... c ... inputs: c ... nsp = number of species c ... x() = mole fractions of species i=1,...,nsp c ... c ... outputs: c ... y() = mass fractions of species i=1,...,nsp c ... c ... notes: c ... x(),y() can have any dimension >= nsp in the calling program c ... c ... do NOT change the sizes of arrays and other parameters defined c ... in the parameter statement unless you do the same everywhere c ... else (including inside certain chemkin routines) c ... common blocks ckwblk and flflfl tie this routine to various c ... work space arrays in chemkin c------------------------------------------------------------------------------- subroutine getmsfr(nsp,x,y) implicit double precision (a-h,o-z) include 'parameter.com' dimension x(nsp),y(nsp) include 'comwk.com' include 'comfl.com' c------------------------------------------------------------------------------- call ckxty (x, iwork(ickw), rwork(nckw), y) c------------------------------------------------------------------------------- return end c******************************************************************************* c ... get mean molecular weight of mixture c ... c ... prerequisites: trinit c ... c ... inputs: c ... nsp = number of species c ... x() = mole fractions of species i=1,...,nsp c ... c ... outputs: c ... wtm = mean molecular weight of mixture, c ... SI units, kg/kg-mol c ... c ... notes: c ... x() can have any dimension >= nsp in the calling program c ... c ... do NOT change the sizes of arrays and other parameters defined c ... in the parameter statement unless you do the same everywhere c ... else (including inside certain chemkin routines) c ... common blocks ckwblk and flflfl tie this routine to various c ... work space arrays in chemkin c------------------------------------------------------------------------------- subroutine getwtm(nsp,x,wref,wtm) implicit double precision (a-h,o-z) include 'parameter.com' dimension x(nsp) include 'comwk.com' include 'comfl.com' c------------------------------------------------------------------------------- call ckmmwx (x, iwork(ickw), rwork(nckw), wtm) #ifdef NON_DIM wtm = wtm / wref #endif c------------------------------------------------------------------------------- return end c******************************************************************************* c ... get local density of mixture c ... c ... prerequisites: trinit c ... c ... inputs: c ... nsp = number of species c ... temp = temperature in degrees K c ... (or) temp = non-dimensional temperature if the NON_DIM flag is on c ... pressure = pressure in SI units, Pascals (N/m^2) c ... (or) pressure = non-dimensional pressure if the NON_DIM flag is on c ... x(i) = mole fractions of species i=1,...,nsp c ... c ... outputs: c ... rho = local mass density of mixture, c ... SI units, kg/m^3 c ... (or) rho = non-dimensional local mass density of mixture if the NON_DIM flag is on c ... c ... notes: c ... x() can have any dimension >= nsp in the calling program c ... c ... do NOT change the sizes of arrays and other parameters defined c ... in the parameter statement unless you do the same everywhere c ... else (including inside certain chemkin routines) c ... common blocks ckwblk and flflfl tie this routine to various c ... work space arrays in chemkin c------------------------------------------------------------------------------- subroutine getrho(temp,tempref,pressure,pref,nsp,x,rhoref,rho) implicit double precision (a-h,o-z) include 'parameter.com' dimension x(nsp) include 'comwk.com' include 'comfl.com' #ifdef NON_DIM templ = temp * tempref pressurel = pressure * pref #else templ = temp pressurel = pressure #endif c------------------------------------------------------------------------------- c ... find pressure in cgs units (dynes/cm^2) c ... 1 N = 10^5 dynes ---> 1 Pa = 10 dynes/cm^2 pcgs = pressurel*10.0d0 c------------------------------------------------------------------------------- c ... call chemkin routine call ckrhox (pcgs, templ, x, iwork(ickw), rwork(nckw), rho) c------------------------------------------------------------------------------- c ... convert rho from cgs to SI rho = rho*1000.0d0 ! convert to SI (kg/m^3) c ... return non-dimensional local density of the mixture #ifdef NON_DIM rho = rho / rhoref #endif c------------------------------------------------------------------------------- return end c******************************************************************************* c ... get binary diffusion coefficients, c ... c ... prerequisites: trinit c ... c ... inputs: c ... nsp = number of species c ... temp = temperature in degrees K c ... pressure = pressure in SI units, Pascals (N/m^2) c ... x(i) = mole fractions of species i=1,...,nsp c ... c ... outputs: c ... bindkj(i,j) = binary diffusion coefficients, SI units, (m^2/s) c ... of species i in species j, (i,j=1,...,nsp) c ... c ... notes: c ... x() can have any dimension >= nsp in the calling program c ... if bindkj() is a 2D array in the calling program (CPO) then its first c ... dimension in the CPO has to be exactly nsp c ... to avoid this restriction, use bindkj() as a 1D vector in the c ... CPO and dimension it with size >= nsp**2. Then, the correspondence c ... between bindkj_here(i,j) here and bindkj_CPO(k) in the CPO is : c ... bindkj_CPO( k ) = bindkj_here(i,j) for k=(j-1)*nsp+i c ... c ... do NOT change the sizes of arrays and other parameters defined c ... in the parameter statement unless you do the same everywhere c ... else (including inside certain chemkin routines) c ... common blocks ckwblk and flflfl tie this routine to various c ... work space arrays in chemkin c------------------------------------------------------------------------------- subroutine getbidij(temp,pressure,nsp,x,bindkj) implicit double precision (a-h,o-z) include 'parameter.com' dimension x(nsp),bindkj(nsp,nsp) include 'comwk.com' include 'comfl.com' c------------------------------------------------------------------------------- c ... find pressure in cgs units (dynes/cm^2) c ... 1 N = 10^5 dynes ---> 1 Pa = 10 dynes/cm^2 pcgs = pressure*10.0d0 c------------------------------------------------------------------------------- c ... call the chemkin routine for binary diffusion coeff. call mcsdif(pcgs,temp,nsp,rwork(nmcw),bindkj(1,1)) c------------------------------------------------------------------------------- c ... bindkj(i,j) are in cgs units: cm^2/s, convert to SI (m^2/s) convert = 1.d-4 do i=1,nsp do j=1,nsp bindkj(i,j)=convert*bindkj(i,j) end do end do c------------------------------------------------------------------------------- return end c******************************************************************************* c ... get cp(i)/R [dimensionless] c ... cp(i) = specific heat at constant pressure of species i, c ... units: J/(kg-mol.K) c ... R = universal gas constant, 8314.0 J/(kg-mol.K) c ... nb. kg-mol = kilogram-mole c ... c ... prerequisites: trinit c ... c ... inputs: c ... nsp = number of species c ... temp = temperature in degrees K c ... c ... outputs: c ... cpor(i) = cp(i)/R, dimensionless, i=1,...,nsp c ... ratio of specific heat of species i (J/kg-mol.K) c ... to universal gas constant (J/kg-mol.K) c ... c ... notes: c ... cpor() can have any dimension >= nsp in the calling program c ... c ... do NOT change the sizes of arrays and other parameters defined c ... in the parameter statement unless you do the same everywhere c ... else (including inside certain chemkin routines) c ... common blocks ckwblk and flflfl tie this routine to various c ... work space arrays in chemkin c------------------------------------------------------------------------------- subroutine getcpor(temp,nsp,cpor) implicit double precision (a-h,o-z) include 'parameter.com' dimension cpor(nsp) include 'comwk.com' include 'comfl.com' c------------------------------------------------------------------------------- c ... call chemkin routine call ckcpor (temp, iwork(ickw), rwork(nckw), cpor) c------------------------------------------------------------------------------- return end c******************************************************************************* c ... get h(i)/(R*T) [dimensionless] c ... h(i) = enthalpy of species i, c ... units: J/(kg-mol) c ... R = universal gas constant, 8314.0 J/(kg-mol.K) c ... T = temperaure (K) c ... nb. kg-mol = kilogram-mole c ... c ... prerequisites: trinit c ... c ... inputs: c ... nsp = number of species c ... temp = temperature in degrees K c ... c ... outputs: c ... hort(i) = h(i)/(R*T), dimensionless, i=1,...,nsp c ... ratio of enthalpy of species i (J/kg-mol) c ... to product of universal gas constant (J/kg-mol.K) c ... and the temperature (K) c ... c ... notes: c ... hort() can have any dimension >= nsp in the calling program c ... c ... do NOT change the sizes of arrays and other parameters defined c ... in the parameter statement unless you do the same everywhere c ... else (including inside certain chemkin routines) c ... common blocks ckwblk and flflfl tie this routine to various c ... work space arrays in chemkin c------------------------------------------------------------------------------- subroutine gethort(temp,nsp,hort) implicit double precision (a-h,o-z) include 'parameter.com' dimension hort(nsp) include 'comwk.com' include 'comfl.com' c------------------------------------------------------------------------------- c ... call chemkin routine call ckhort (temp, iwork(ickw), rwork(nckw), hort) c------------------------------------------------------------------------------- return end c******************************************************************************* c ... get the pure species viscosities c ... c ... prerequisites: trinit c ... c ... inputs: c ... nsp = number of species c ... temp = temperature in degrees K c ... c ... outputs: c ... spvis(i) = viscosity of the pure species i, i=1,...,nsp c ... units, SI, kg/m.s c ... notes: c ... spvis() can have any dimension >= nsp in the calling program c ... c ... do NOT change the sizes of arrays and other parameters defined c ... in the parameter statement unless you do the same everywhere c ... else (including inside certain chemkin routines) c ... common blocks ckwblk and flflfl tie this routine to various c ... work space arrays in chemkin c------------------------------------------------------------------------------- subroutine getspvis(temp,nsp,spvis) implicit double precision (a-h,o-z) include 'parameter.com' dimension spvis(nsp) include 'comwk.com' include 'comfl.com' c------------------------------------------------------------------------------- c ... call chemkin routine call mcsvis (temp,rwork(nmcw),spvis) c------------------------------------------------------------------------------- c ... spvis is in cgs units : gm/cm.s, convert to SI : kg/m.s convert=0.1d0 do i=1,nsp spvis(i)=convert*spvis(i) end do c------------------------------------------------------------------------------- return end c******************************************************************************* c ... get the viscosity of the mixture c ... c ... prerequisites: trinit c ... c ... inputs: c ... nsp = number of species c ... temp = temperature in degrees K c ... x(i) = mole fractions, species i=1,...,nsp c ... c ... outputs: c ... vismix = viscosity of the mixture c ... units, SI, kg/m.s c ... notes: c ... c ... do NOT change the sizes of arrays and other parameters defined c ... in the parameter statement unless you do the same everywhere c ... else (including inside certain chemkin routines) c ... common blocks ckwblk and flflfl tie this routine to various c ... work space arrays in chemkin c------------------------------------------------------------------------------- subroutine getmxvis(temp,nsp,x,vismix) implicit double precision (a-h,o-z) include 'parameter.com' dimension x(nsp) include 'comwk.com' include 'comfl.com' c------------------------------------------------------------------------------- c ... call chemkin routine call mcavis (temp,x,rwork(nmcw),vismix) c------------------------------------------------------------------------------- c ... vismix is in cgs units : gm/cm.s, convert to SI : kg/m.s convert=0.1d0 vismix = vismix*convert c------------------------------------------------------------------------------- return end c******************************************************************************* c ... get the pure species thermal conductivities c ... c ... prerequisites: trinit c ... c ... inputs: c ... nsp = number of species c ... temp = temperature in degrees K c ... c ... outputs: c ... spcon(i) = thermal conductivity of the pure species i, i=1,...,nsp c ... units, SI, (W/m.K) c ... notes: c ... spcon() can have any dimension >= nsp in the calling program c ... c ... do NOT change the sizes of arrays and other parameters defined c ... in the parameter statement unless you do the same everywhere c ... else (including inside certain chemkin routines) c ... common blocks ckwblk and flflfl tie this routine to various c ... work space arrays in chemkin c------------------------------------------------------------------------------- subroutine getspcon(temp,nsp,spcon) implicit double precision (a-h,o-z) include 'parameter.com' dimension spcon(nsp) include 'comwk.com' include 'comfl.com' c------------------------------------------------------------------------------- c ... call chemkin routine call mcscon (temp,rwork(nmcw),spcon) c------------------------------------------------------------------------------- c ... spcon() is in cgs units: erg/cm.K.s c ... convert to SI: J/m.K.s convert = 1.d-5 do i=1,nsp spcon(i)=convert*spcon(i) end do c------------------------------------------------------------------------------- return end c******************************************************************************* c ... get the thermal conductivity of the mixture c ... c ... prerequisites: trinit c ... c ... inputs: c ... nsp = number of species c ... temp = temperature in degrees K c ... x(i) = mole fractions, species i=1,...,nsp c ... c ... outputs: c ... conmix = thermal conductivity of the mixture c ... units, SI, kg/m.s c ... notes: c ... c ... this routine should not be used when the multicomponent diffusion c ... coefficient implementation is used, rather the condpr computed c ... from the same routine that computes the multicomponent thermal c ... diffusion coefficients should be used in the energy equation c ... c ... do NOT change the sizes of arrays and other parameters defined c ... in the parameter statement unless you do the same everywhere c ... else (including inside certain chemkin routines) c ... common blocks ckwblk and flflfl tie this routine to various c ... work space arrays in chemkin c------------------------------------------------------------------------------- subroutine getmxcon(temp,nsp,x,conmix) implicit double precision (a-h,o-z) include 'parameter.com' dimension x(nsp) include 'comwk.com' include 'comfl.com' c------------------------------------------------------------------------------- c ... call chemkin routine call mcacon (temp,x,rwork(nmcw),conmix) c------------------------------------------------------------------------------- c ... conmix is in cgs units: erg/cm.K.s c ... convert to SI: J/m.K.s convert = 1.d-5 conmix = convert*conmix c------------------------------------------------------------------------------- return end c******************************************************************************* c ... get the mass production rate wms(i) of species i=1,...,nsp c ... c ... prerequisites: trinit, and getwt (to provide the wt() input vector) c ... c ... inputs: c ... nsp = number of species c ... temp = temperature in degrees K c ... pressure = pressure in SI units, Pascals (N/m^2) c ... x(i) = mole fractions of species i=1,...,nsp c ... wt(i) = molar weights of species i=1,...,nsp, (kg/kg-mol) c ... c ... outputs: c ... wms(i) = rate of mass production of species i, i=1,...,nsp c ... SI units, (kg/m^3.s) c ... c ... notes: c ... wms() can have any dimension >= nsp in the calling program c ... c ... do NOT change the sizes of arrays and other parameters defined c ... in the parameter statement unless you do the same everywhere c ... else (including inside certain chemkin routines) c ... common blocks ckwblk and flflfl tie this routine to various c ... work space arrays in chemkin c------------------------------------------------------------------------------- subroutine getwms(temp,pressure,nsp,x,wt,wms) implicit double precision (a-h,o-z) include 'parameter.com' dimension x(nsp),wms(nsp),wt(nsp) include 'comwk.com' include 'comfl.com' c------------------------------------------------------------------------------- c ... find pressure in cgs units (dynes/cm^2) c ... 1 N = 10^5 dynes ---> 1 Pa = 10 dynes/cm^2 pcgs = pressure*10.0d0 c------------------------------------------------------------------------------- c ... call the chemkin routine call ckwxp (pcgs,temp,x,iwork(ickw),rwork(nckw),wms) c------------------------------------------------------------------------------- c ... wms() returns in gm-mol/cm^3.s c ... convert to gm/cm^3.s by multiplying by wt(i), [wt(i)]=gm/gm-mol=kg/kg-mol c ... then to kg/m^3.s by multiplying by 1.d-3/1.d-6=1.d3 convert = 1.d3 do i=1,nsp if(wt(i).eq.0)then write(*,*)' fill out wt() array before calling getwms ' stop endif wms(i)=wms(i)*wt(i)*convert end do c------------------------------------------------------------------------------- return end c******************************************************************************* c ... get mole fractions of species in mixture c ... c ... prerequisites: trinit c ... c ... inputs: c ... nsp = number of species c ... y() = mass fractions of species i=1,...,nsp c ... c ... outputs: c ... x() = mole fractions of species i=1,...,nsp c ... c ... notes: c ... x(),y() can have any dimension >= nsp in the calling program c ... c ... do NOT change the sizes of arrays and other parameters defined c ... in the parameter statement unless you do the same everywhere c ... else (including inside certain chemkin routines) c ... common blocks ckwblk and flflfl tie this routine to various c ... work space arrays in chemkin c------------------------------------------------------------------------------- subroutine getmlfr(nsp,y,x) implicit double precision (a-h,o-z) include 'parameter.com' dimension x(nsp),y(nsp) include 'comwk.com' include 'comfl.com' c------------------------------------------------------------------------------- c ... call the chemkin routine call ckytx (y, iwork(ickw), rwork(nckw), x) c------------------------------------------------------------------------------- return end c******************************************************************************* c ... get cp of mixture c ... c ... prerequisites: trinit c ... c ... inputs: c ... temp = temperature in degrees K c ... (or) temp = non-dimensional temperature if the NON_DIM flag is on c ... nsp = number of species c ... y(i) = mass fractions of species i=1,...,nsp c ... c ... outputs: c ... cpmix = cp of the mixture at given (T,Y()) c ... SI units, J/kg.K c ... (or) cpmix = cp of the mixture at given (T,Y()) non-dimensional. c ... notes: c ... y() can have any dimension >= nsp in the calling program c ... c ... do NOT change the sizes of arrays and other parameters defined c ... in the parameter statement unless you do the same everywhere c ... else (including inside certain chemkin routines) c ... common blocks ckwblk and flflfl tie this routine to various c ... work space arrays in chemkin c------------------------------------------------------------------------------- subroutine getcpmix(temp,tempref,nsp,y,cpref,cpmix) implicit double precision (a-h,o-z) include 'parameter.com' dimension y(nsp) include 'comwk.com' include 'comfl.com' #ifdef NON_DIM templ = temp * tempref #else templ = temp #endif c------------------------------------------------------------------------------- c ... call chemkin routine call ckcpbs (templ, y, iwork(ickw), rwork(nckw), cpmix) c------------------------------------------------------------------------------- c ... convert cpmix from cgs (ergs/gm.K) to SI c ... 1J = 10^7 ergs cpmix = cpmix*1.d-4 ! convert to SI (J/kg.K) c------------------------------------------------------------------------------- c ... return non-dimensional cp_mix #ifdef NON_DIM cpmix = cpmix / cpref #endif return end c******************************************************************************* c ... get cv of mixture c ... c ... prerequisites: trinit c ... c ... inputs: c ... temp = temperature in degrees K c ... (or) temp = non-dimensional temperature if the NON_DIM flag is on c ... nsp = number of species c ... y(i) = mass fractions of species i=1,...,nsp c ... c ... outputs: c ... cvmix = cv of the mixture at given (T,Y()) c ... SI units, J/kg.K c ... (or) cvmix = cv of the mixture at given (T,Y()) non-dimensional. c ... c ... notes: c ... y() can have any dimension >= nsp in the calling program c ... c ... do NOT change the sizes of arrays and other parameters defined c ... in the parameter statement unless you do the same everywhere c ... else (including inside certain chemkin routines) c ... common blocks ckwblk and flflfl tie this routine to various c ... work space arrays in chemkin c------------------------------------------------------------------------------- subroutine getcvmix(temp,tempref,nsp,y,cpref,rgref,cvmix) implicit double precision (a-h,o-z) include 'parameter.com' dimension y(nsp) include 'comwk.com' include 'comfl.com' #ifdef NON_DIM cvref = cpref - rgref templ = temp * tempref #else templ = temp #endif c------------------------------------------------------------------------------- c ... call chemkin routine call ckcvbs (templ, y, iwork(ickw), rwork(nckw), cvmix) c------------------------------------------------------------------------------- c ... convert cvmix from cgs (ergs/gm.K) to SI c ... 1J = 10^7 ergs cvmix = cvmix*1.d-4 ! convert to SI (J/kg.K) c------------------------------------------------------------------------------- c ... return non-dimensional cv_mix #ifdef NON_DIM cvmix = cvmix / cvref #endif return end c******************************************************************************* c ... get h of mixture c ... c ... prerequisites: trinit c ... c ... inputs: c ... temp = temperature in degrees K c ... (or) temp = non-dimensional temperature if the NON_DIM flag is on c ... nsp = number of species c ... y(i) = mass fractions of species i=1,...,nsp c ... c ... outputs: c ... hmix = h of the mixture at given (T,Y()) c ... SI units, J/kg c ... (or) hmix = h of the mixture at given (T,Y()) non-dimensional. c ... c ... notes: c ... y() can have any dimension >= nsp in the calling program c ... c ... do NOT change the sizes of arrays and other parameters defined c ... in the parameter statement unless you do the same everywhere c ... else (including inside certain chemkin routines) c ... common blocks ckwblk and flflfl tie this routine to various c ... work space arrays in chemkin c------------------------------------------------------------------------------- subroutine gethmix(temp,tempref,nsp,y,href,hmix) implicit double precision (a-h,o-z) include 'parameter.com' dimension y(nsp) include 'comwk.com' include 'comfl.com' #ifdef NON_DIM templ = temp * tempref #else templ = temp #endif c------------------------------------------------------------------------------- c ... call chemkin routine call ckhbms (templ, y, iwork(ickw), rwork(nckw), hmix) c------------------------------------------------------------------------------- c ... convert hmix from cgs (ergs/gm) to SI c ... 1J = 10^7 ergs hmix = hmix*1.d-4 ! convert to SI (J/kg) c------------------------------------------------------------------------------- c ... return non-dimensional h_mix #ifdef NON_DIM hmix = hmix / href #endif return end c******************************************************************************* c ... get gsp(i) [J/kg] c ... gsp(i) = hsp(i)-T*ssp(i) c ... where hsp(i) is the specific enthalpy of species i c ... and ssp(i) is the specific entropy of species i c ... NOTE: while the Chemkin manual calls this the 'standard state c ... Gibbs free energy', this is misleading, because the tabulation c ... in chemkin uses different 0-bases for hsp and ssp. c ... Namely, hsp of a naturally occuring element at 298K is zero c ... but not so for ssp. In fact, gsp=hsp-T*ssp, does not agree c ... with tabulated 'standard free energy of FORMATION', this c ... latter requires a definition with consistent hsp and ssp bases. c ... Nevertheless, the quantity gsp in Chemkin is useful, and is c ... certainly meaningful in the difference evaluations across reactions c ... since the effect of the different bases disappears. c ... c ... So, from here on we'll refer to gsp (with reservation) as the : c ... standard state Gibbs free energy, in mass units, of species i, c ... units: J/kg c ... other quants: c ... R = universal gas constant, 8314.0 J/(kg-mol.K) c ... T = temperaure (K) c ... nb. kg-mol = kilogram-mole c ... c ... prerequisites: trinit c ... c ... inputs: c ... nsp = number of species c ... temp = temperature in degrees K c ... c ... outputs: c ... gsp(i), i=1,...,nsp c ... standard state Gibbs free energy, in SI units (J/kg), of species i c ... c ... notes: c ... gsp() can have any dimension >= nsp in the calling program c ... c ... do NOT change the sizes of arrays and other parameters defined c ... in the parameter statement unless you do the same everywhere c ... else (including inside certain chemkin routines) c ... common blocks ckwblk and flflfl tie this routine to various c ... work space arrays in chemkin c------------------------------------------------------------------------------- subroutine getgibbs(temp,nsp,gsp) implicit double precision (a-h,o-z) include 'parameter.com' dimension gsp(nsp) include 'comwk.com' include 'comfl.com' c------------------------------------------------------------------------------- c ... call chemkin routine call ckgms (temp, iwork(ickw), rwork(nckw), gsp) c------------------------------------------------------------------------------- c ... convert gsp from cgs (ergs/gm) to SI (J/kg) c ... 1J = 10^7 ergs convert = 1.d-4 do i=1,nsp gsp(i) = gsp(i)*convert ! convert to SI (J/kg) end do c------------------------------------------------------------------------------- return end c******************************************************************************* c c ADDED SUBROUTINES c c******************************************************************************* c ... get the mean entropy of mixture sbms (J/Kg*K) subroutine getsbms(nsp,pressure, temp, y, sbms) implicit double precision (a-h,o-z) include 'parameter.com' dimension y(nsp) include 'comwk.com' include 'comfl.com' c------------------------------------------------------------------------------- c ... find pressure in cgs units (dynes/cm^2) c ... 1 N = 10^5 dynes ---> 1 Pa = 10 dynes/cm^2 pcgs = pressure*10.0d0 c------------------------------------------------------------------------------- c------------------------------------------------------------------------------- c ... call chemkin routine call cksbms (pcgs, temp, y, iwork(ickw), rwork(nckw), sbms_cgs) c------------------------------------------------------------------------------- c ... sbms is in cgs units (ergs/(gm*K)) c ... convert to SI (J/Kg K) sbms = sbms_cgs * (1.0d-4) return end c****************************************************************************** c ... get the standard state entropies (J/Kg*K) subroutine getsms(nsp,temp,sms) implicit double precision (a-h,o-z) include 'parameter.com' dimension sms(nsp) include 'comwk.com' include 'comfl.com' c------------------------------------------------------------------------------- c ... call chemkin routine call cksms (temp, iwork(ickw), rwork(nckw), sms) c------------------------------------------------------------------------------- c ... sms is in cgs units (ergs/(gm*K)) c ... convert to SI (J/Kg K) const = 1.0d-4 do i=1,nsp sms(i) = sms(i) * const end do return end c****************************************************************************** c ... get total element & reaction count c ... mm total elements c ... ii total reactions subroutine getindx(mm,ii) implicit double precision (a-h,o-z) include 'parameter.com' include 'comwk.com' include 'comfl.com' c------------------------------------------------------------------------------- c ... call chemkin routine call ckindx (iwork(ickw), rwork(nckw), mm, ns_indx, ii, nfit) c------------------------------------------------------------------------------- return end c******************************************************************************* c ... get matrix NCF(MM,KK) containing the elemental c ... composition of the species. c ... NCF(M,K) is the no. of atoms of element M in species K subroutine getncf(ns_indx,mm,lda,ncf) implicit double precision (a-h,o-z) include 'parameter.com' include 'comwk.com' include 'comfl.com' dimension NCF(lda,*) c------------------------------------------------------------------------------- c ... call chemkin routine call ckncf(lda,iwork(ickw),rwork(nckw),ncf) c------------------------------------------------------------------------------- return end c******************************************************************************* c ... temp = temperature in degrees K c ... (or) temp = non-dimensional temperature if the NON_DIM flag is on c ... nsp = number of species c ... OUTPUT: c ... cpms(i) = specific heat at constant pressure of species i, c ... units: J/(kg-mol.K) c ... (or) cpms(i) = non-dimensional specific heat at constant pressure of species i, c ... if the NON_DIM flag is on c ... nb. kg-mol = kilogram-mole c ... For chemkin it is: c ... cpms in (ergs/(gm*K)) subroutine getcpms(temp,tempref,nsp,cpref,cpms) implicit double precision (a-h,o-z) include 'parameter.com' dimension cpms(nsp) include 'comwk.com' include 'comfl.com' #ifdef NON_DIM templ = temp * tempref #else templ = temp #endif c------------------------------------------------------------------------------- c ... call chemkin routine call ckcpms (templ, iwork(ickw), rwork(nckw), cpms) c------------------------------------------------------------------------------- c ... spms is in cgs units (ergs/(gm*K)) c ... convert to SI (J/Kg K) do i=1,nsp cpms(i) = cpms(i) * (1.0d-4) end do #ifdef NON_DIM do i=1,nsp cpms(i) = cpms(i) / cpref end do #endif return end c****************************************************************************** c ... get the species' enthalpies (J/Kg) c ... temp = temperature in degrees K c ... (or) temp = non-dimensional temperature if the NON_DIM flag is on c ... nsp = number of species c ... OUTPUT: c ... hms(i) = species enthalpies, units: J/kg c ... (or) hms(i) = non-dimensional species enthalpies, if the NON_DIM flag is on subroutine gethms(temp,tempref,nsp,href,hms) implicit double precision (a-h,o-z) include 'parameter.com' dimension hms(nsp) include 'comwk.com' include 'comfl.com' #ifdef NON_DIM templ = temp * tempref #else templ = temp #endif c------------------------------------------------------------------------------- c ... call chemkin routine call ckhms (templ, iwork(ickw), rwork(nckw), hms) c------------------------------------------------------------------------------- c ... hms is in cgs units (ergs/gm) c ... convert to SI (J/Kg) do i=1,nsp hms(i) = hms(i) * (1.0d-4) end do #ifdef NON_DIM do i=1,nsp hms(i) = hms(i) / href end do #endif return end c****************************************************************************** c ... get the species' enthalpies (J/Kg-mole) c ... temp = temperature in degrees K c ... (or) temp = non-dimensional temperature if the NON_DIM flag is on c ... nsp = number of species c ... OUTPUT: c ... hml(i) = species enthalpies, units: J/kg-mole c ... (or) hml(i) = non-dimensional species enthalpies, if the NON_DIM flag is on subroutine gethml(temp,tempref,nsp,wref,href,hml) implicit double precision (a-h,o-z) include 'parameter.com' dimension hml(nsp) include 'comwk.com' include 'comfl.com' #ifdef NON_DIM templ = temp * tempref #else templ = temp #endif c------------------------------------------------------------------------------- c ... call chemkin routine call ckhml (templ, iwork(ickw), rwork(nckw), hml) c------------------------------------------------------------------------------- c ... hml is in cgs units (ergs/mole) c ... convert to SI (J/Kg-mole) do i=1,nsp hml(i) = hml(i) * (1.0d-4) end do #ifdef NON_DIM do i=1,nsp hml(i) = hml(i) / (href*wref) end do #endif return end c******************************************************************************* c ... Get mean Gibbs Energy in Mass units (J/Kg) c------------------------------------------------------------------------------- subroutine getgbms(pressure,temp,Y,nsp,gbms) implicit double precision (a-h,o-z) include 'parameter.com' dimension Y(nsp) include 'comwk.com' include 'comfl.com' c------------------------------------------------------------------------------- c ... find pressure in cgs units (dynes/cm^2) c ... 1 N = 10^5 dynes ---> 1 Pa = 10 dynes/cm^2 pcgs = pressure*10.0d0 c------------------------------------------------------------------------------- c------------------------------------------------------------------------------- c ... call chemkin routine call ckgbms (pcgs,temp,Y,iwork(ickw),rwork(nckw),gbmscgs) c------------------------------------------------------------------------------- c ... convert gsp from cgs (ergs/gm) to SI (J/kg) c ... 1J = 10^7 ergs convert = 1.d-4 gbms = gbmscgs*convert ! convert to SI (J/kg) c------------------------------------------------------------------------------- return end