subroutine getrates(pressure,n_slack_species_index, 1 nspecies,ncells,n_size,data_in,data_out) IMPLICIT NONE include 'parameter.par' include 'chem.com' include 'chem_ref.com' integer n_slack_species_index,nspecies,ncells,n_size integer nsmxl,nrmx parameter (nrmx=216) !maximum number of reactions real*8 temp(ncells),rho(ncells),ys(nspecies,ncells), 1 yc(ncells),cspl(nspecies),wspl(nspecies),wcon(nspecies), 2 arop(nrmx),wcheml(nspecies,ncells),hspl(nspecies), 3 ewt(ncells),cpspl(nspecies) real*8 data_in(n_size) !Solution vector in [Temp, Y(species-1)] real*8 data_out(n_size)!Solution vector out !d[Temp, Y(species-1)]/dt real*8 pressure_in(ncells) real*8 pressure(ncells) integer n,j,k double precision sum,V,rho_k,rop,cpc n=ncells ns=nspecies c--------------------------------------------------------------------------- c assign the input data c do j=1,n c pressure(j) = pressure_in(j) c enddo do j=1,n ! j=cell id Sum = 0.0 #ifdef NON_DIM temp(j) = data_in( (j-1) * (nspecies) + 1 ) #else temp(j) = data_in( (j-1) * (nspecies) + 1 ) / tempref #endif #ifdef ALL_SPECIES do k=1,ns ys(k,j) = data_in( (j-1)*(nspecies)+ (1+k) ) Sum = Sum + ys(k,j) end do c Make sure that conservation of mass holds true if (Sum.ne.1.d0) then do k=1,ns ys(k,j) = ys(k,j) / Sum end do endif #else c Compute the species to take the blame from the rest do k=1,(ns-1) Sum = Sum + data_in ( (j-1)*(nspecies)+ (1+k) ) end do ys(n_slack_species_index,j) = 1.0 - Sum c Make sure that no negative mass fraction appears if ( ys(n_slack_species_index,j).lt.0.0 ) then ys(n_slack_species_index,j) = 0.0 end if c Fill in the mass fractions of the rest of the species if (n_slack_species_index.eq.1) then do k=2,ns ys(k,j) = data_in( (j-1)*(nspecies)+ (k) ) end do elseif (n_slack_species_index.eq.ns) then do k=1,(ns-1) ys(k,j) = data_in( (j-1)*(nspecies)+ (1+k) ) end do else do k=1,(n_slack_species_index-1) ys(k,j) = data_in( (j-1)*(nspecies)+ (1+k) ) end do do k=(n_slack_species_index+1),ns ys(k,j) = data_in( (j-1)*(nspecies)+ (k) ) end do endif #endif end do c------------------------------------------------------------------------------- c Evaluate the cell density of the mixture: c 1. rho(mix) = m(mix) / V(mix) c 2. m(mix) = 1 kgr c 3. V(mix) = V(1) + V(2) +...+ V(ns), for species 1..ns c 4. V(k) = m(k) / rho(k) c 5. m(k) = ys(k,i,j) * m(mix) = ys(k,i,j) c 6. rho(k) = pressure / ( rg(k)*temp(i,j)*tempref ) c 7. rg(k) = r0 / wt do j=1,n V = 0.0 sum = 0.0 c#ifdef NON_DIM c pressure(j) = pressure(j) * pref c#endif do k = 1,ns rg(k) = r0 / (1000.0 * wt(k)) #ifdef NON_DIM rho_k= pressure(j)*pref / (1000.0 * rg(k) * temp(j) * tempref ) #else rho_k= pressure(j) / (1000.0 * rg(k) * temp(j) * tempref ) #endif V = ys(k,j) / rho_k sum = sum + V end do c non-dimensionalize rho, since that's how it's used later. rho(j) = (1.0 / sum) / rhoref end do c------------------------------------------------------------------------------- c find the rop's do j=1,n call y2c(rho(j),ys(1,j),cspl) do k=1,numreac arop(k) = rop(k,temp(j),tempref,cspl) end do end do c find wchem's do j=1,n call wchem(temp(j),tempref,cspl, wcheml(1,j)) end do c=========================temperature source term============================= c ... hspl(k) is the enthalpy of species k at (i,j) do j=1,n c ... linear interpolation from table of h_k(T) call hsp(temp(j),hspl) c ... wspl(k) = wchem(k,j) is the rate of creation of species k at (j) sum = 0.0 do k=1,ns #ifdef NON_DIM sum = sum + damkohlr * hspl(k) * wcheml(k,j) #else sum = sum + hspl(k) * wcheml(k,j) #endif end do ewt(j) = - sum c ... linear interpolation from table of cp_k(T) call cpsp(temp(j),cpspl) !!corresponds to each species k on cell (i,j) c ... cpc is the specific heat for the cell sum = 0.0 do k=1,ns sum = sum + ys(k,j)*cpspl(k) end do cpc = sum end do #ifdef NON_DIM c============================ Create the non-dimensional output =========================== do j=1,n c ... w_T = -Da*Sum(h_k*wchem) + g(t) data_out( (j-1) * (nspecies)+1 )=ewt(j)/(rho(j)*cpc) #ifdef ALL_SPECIES do k=1,ns data_out ( (j-1)*(nspecies)+ (1+k) ) = damkohlr* 1 wcheml(k,j)/rho(j) end do #else if (n_slack_species_index.eq.1) then do k=2,ns data_out( (j-1)*(nspecies)+ (k) ) = damkohlr* 1 wcheml(k,j)/rho(j) end do elseif (n_slack_species_index.eq.ns) then do k=1,(ns-1) data_out( (j-1)*(nspecies)+ (1+k) ) = damkohlr* 1 wcheml(k,j)/rho(j) end do else do k=1,(n_slack_species_index-1) data_out( (j-1)*(nspecies)+ (1+k) ) = damkohlr* 1 wcheml(k,j)/rho(j) end do do k=(n_slack_species_index+1),ns data_out ( (j-1)*(nspecies)+ (k) ) = damkohlr* 1 wcheml(k,j)/rho(j) end do endif #endif end do #else c============================ Create the dimensional output =========================== do j=1,n c ... w_T = -Sum(h_k*wchem) + g(t) in [K/sec] data_out( (j-1) * (nspecies)+1 ) = (ewt(j)* 1 (tempref*wspref)/(rho(j)*rhoref*cpc)) c ... dY/dt in [1/sec] #ifdef ALL_SPECIES do k=1,ns data_out ( (j-1)*(nspecies)+ (1+k) ) = (wcheml(k,j)* 1 wspref)/(rho(j)*rhoref) end do #else if (n_slack_species_index.eq.1) then do k=2,ns data_out ( (j-1)*(nspecies)+ (k) ) = (wcheml(k,j)* 1 wspref)/(rho(j)*rhoref) end do elseif (n_slack_species_index.eq.ns) then do k=1,(ns-1) data_out ( (j-1)*(nspecies)+ (1+k) ) = (wcheml(k,j)* 1 wspref)/(rho(j)*rhoref) end do else do k=1,(n_slack_species_index-1) data_out ( (j-1)*(nspecies)+ (1+k) ) = (wcheml(k,j)* 1 wspref)/(rho(j)*rhoref) end do do k=(n_slack_species_index+1),ns data_out ( (j-1)*(nspecies)+ (k) ) = (wcheml(k,j)* 1 wspref)/(rho(j)*rhoref) end do endif #endif end do c========================================================================== #endif !NON_DIM flag ends here return end