c ... code to compute d(z)/dt for dvode
c ... yv ={ rho.Y_1, rho.Y_2, ... rho.Y_N, T, rho, P}

c ... If we are solving for (ns-1) species it is:
c ... z  ={ rho.Y_1, rho.Y_2, ... rho.Y_N-1, rho, P} ... reduced vector

c ... Otherwise if we are solving for all the species (ns) it is:
c ... z  ={ rho.Y_1, rho.Y_2, ... rho.Y_N, rho, P} ... reduced vector

      subroutine gzdot(neq,z,zdot,ns,
     + nsmx_v,idim_v,ntmx_v,maxsp_v,maxtb_v,nfar_v,npar_v,
     + xl10,cspl_diag,dtc, 
     + cgspl,hspl,cpspl,yv,rr_f,rr_r,aaf,ropl,qfac,ropl1,cthb,rklow,pr,
     + prlog,fc,xp,arg,aa,qff,bb,clog,oclog,fcent,yvdot,yy,wspl,cspl,
     + xik,rc,wbar,tc,otc,vtemp,vlntemp,rtab,ctot,vt,vot,capsp,hwt,
     + sumw,wbarbiro,cp,p0l,orcp,econ,
     + jtab,jtp,tempref,wt,rg,cgref,cpref,href,
     + ospwt,gcon,rexpo,rcape,
     + numreac,itfdep,par,cape,nspec,nprod,kflagrev,itrdep,
     + nkk,anuk,kspk,anusum,nreac,nu,nunk,nthb,ithb,
     + ntbs,nktb,aik,nfal,ifal,ifop,kfal,capel,fpar,damspwt,
     + capb,rcapb,capbl,wcheml,damkohlr,ewtl,
     + temptab,temp1tab,odtab,cgsptab,cpsptab,hsptab,
     + n_slack_species_index)

      IMPLICIT NONE

c     Global Variables

      integer neq,ns,nsmx_v,idim_v,ntmx_v,maxsp_v,maxtb_v
      integer nfar_v,npar_v,jtab,jtp,numreac,nthb,nfal
      integer n_slack_species_index

      integer itfdep(idim_v),nspec(idim_v),nprod(idim_v)
      integer kflagrev(idim_v),itrdep(idim_v),nkk(idim_v)
      integer kspk(maxsp_v,idim_v),nreac(idim_v)
      integer nu(maxsp_v,idim_v),nunk(maxsp_v,idim_v)
      integer ithb(idim_v),ntbs(idim_v),nktb(maxtb_v,idim_v)
      integer ifal(idim_v),ifop(idim_v),kfal(idim_v)


      double precision xl10,dtc,rc,wbar,tc,otc,vtemp,vlntemp
      double precision rtab,ctot,vt,vot,capsp,hwt,sumw
      double precision wbarbiro,cp,p0l,orcp,tempref,cgref
      double precision cpref,href,gcon,damkohlr,ewtl
      double precision temp1tab,odtab

c      real*8 z(1),zdot(1)
      real*8 z(nsmx_v+2),zdot(nsmx_v+2),cspl_diag(nsmx_v)
      real*8 cgspl(nsmx_v),hspl(nsmx_v),cpspl(nsmx_v)
      real*8 yv(nsmx_v+3),rr_f(idim_v),rr_r(idim_v)
      real*8 aaf(idim_v),ropl(idim_v),qfac(idim_v)
      real*8 ropl1(idim_v),cthb(idim_v),rklow(idim_v)
      real*8 pr(idim_v),prlog(idim_v),fc(idim_v),xp(idim_v)
      real*8 arg(idim_v),aa(idim_v),qff(idim_v),bb(idim_v)
      real*8 clog(idim_v),oclog(idim_v),fcent(idim_v)
      real*8 yvdot(nsmx_v+3),yy(nsmx_v),wspl(nsmx_v)
      real*8 cspl(nsmx_v),xik(idim_v),econ(5,idim_v)
      real*8 wt(nsmx_v),rg(nsmx_v),ospwt(nsmx_v),rexpo(idim_v)
      real*8 rcape(idim_v),par(npar_v,idim_v),cape(idim_v)
      real*8 anuk(maxsp_v,idim_v),anusum(idim_v)
      real*8 aik(maxtb_v,idim_v),capel(idim_v)
      real*8 fpar(nfar_v,idim_v),damspwt(nsmx_v)
      real*8 capb(idim_v),rcapb(idim_v),capbl(idim_v)
      real*8 wcheml(nsmx_v),temptab(ntmx_v)
      real*8 cgsptab(nsmx_v,ntmx_v),cpsptab(nsmx_v,ntmx_v)
      real*8 hsptab(nsmx_v,ntmx_v)	

c     Local Variables

      double precision dln10,sum
      double precision cgsv,cpv,hsv,xikv,crprod,cvp,crp1
      double precision cpprod,cpp1
      double precision fclog,xn,flog,pr1,prcm,pcor


      integer kfl,kfc,ifp,kthb,ispthb,ip3,nyp,isp
      integer i,k,nk,ksp
      integer nnn,ir

      data dln10/2.3025850929940459d0/


c-------------------------------------------------------------------------------
c	--------------------------------------------------------
c		---------------------------------------
c				start
c-------------------------------------------------------------------------------
c-------------------------------------------------------------------------------
c-------------------------------------------------------------------------------
c  ...	decompress vector z into vector yv as explained above

      sum = 0.0d0

#ifdef ALL_SPECIES
      yv(ns+2) = z(ns+1)  !density term 
      yv(ns+3) = z(ns+2)  !pressure term

      do i=1,ns
         yv(i) = z(i)
      end do

#else
      yv(ns+2) = z(ns)    !density term 
      yv(ns+3) = z(ns+1)  !pressure term

      if (n_slack_species_index.eq.1) then  
         do i=2,ns
            yv(i) = z(i-1)
            sum = sum + z(i-1)
         end do
         yv(n_slack_species_index) = z(ns) - sum

      elseif (n_slack_species_index.eq.ns) then 
         do i=1,ns-1
            yv(i) = z(i)
           sum = sum + z(i)
         end do
         yv(n_slack_species_index) = z(ns) - sum

      else
         do i=1,(n_slack_species_index-1)
            yv(i) = z(i)
            sum = sum + z(i)
         end do
 
         do i=(n_slack_species_index+1),ns
            yv(i) = z(i-1)
            sum = sum + z(i-1)
         end do
         yv(n_slack_species_index) = z(ns) - sum
 
      endif           

#endif

      rc = yv(ns+2)  !density term 

c  ...	and now for T

      do i=1,ns
         yy(i) = yv(i)/rc
      end do

      wbar = 0.
      do i=1,ns
         wbar = wbar + yy(i)*ospwt(i)
      end do
      wbar = 1./wbar

c  ...	reconstruct temperature using the equation of state
      p0l = yv(ns+3)    !pressure term
      tc  = p0l * wbar / rc
      yv(ns+1) = tc     !temperature term

      otc = 1./tc
      vlntemp = log(tc)
      vtemp   = tc*tempref

c  ...  get gibbs free energies
c  ...	get c_{p,i} for each of the substances 
c  ...	get enthalpies for each of the species
c-------------------------------------------------------------------------------

c ... linear interpolation from tables
	    
      jtab = int((tc-temp1tab)*odtab) + 1
      rtab = (tc-temptab(jtab))*odtab
      jtp  = jtab + 1

      do k=1,ns
         cgsv     = cgsptab(k,jtab)
         cgspl(k) = cgsv + rtab*(cgsptab(k,jtp)-cgsv)     
         cpv      = cpsptab(k,jtab)
         cpspl(k) = cpv + rtab*(cpsptab(k,jtp)-cpv)     
         hsv      = hsptab(k,jtab)
         hspl(k)  = hsv + rtab*(hsptab(k,jtp)-hsv)     
      end do

c-------------------------------------------------------------------------------
c------------------------------------------------------------------------------
c	PART II		PART II		PART II		PART II	
c------------------------------------------------------------------------------
c
c			species Source Terms
c
c------------------------------------------------------------------------------
c	...this should follow closely getrates and differences...
c	...should primarily be due to code inserted to compute...
c	...     derivatives which enter into the Jacobian      ...
c
c ... initialize

	do ksp=1,ns
	   wspl(ksp)  = 0.0
	   cspl(ksp)  = yv(ksp)*ospwt(ksp)	! concentrations
	end do

c-------------------------------------------------------------------------------
c ... loop over all reactions to evaluate rate-of-progress: ropl(k)

	do k = 1, numreac
	      
c ... evaluate forward reaction rate, rr_f

	   if(itfdep(k).ne.0)then
	      rr_f(k) = capb(k)*exp(par(2,k)*vlntemp - cape(k)*otc)
	   else
	      rr_f(k) = capb(k)
	   endif

c-------------------------------------------------------------------------------
c ... nspec(k) > 0 iff reaction is reversible
c ... xik = SUM_i=1,ns nu(i,k)*G_ik
c ... evaluate (1/equilibrium_reaction_rate) = orr_e
c ... evaluate reverse reaction rate, rr_r
c ... kflagrev(k) = 1 -> reverse Arrhenius rate parameters are defined
c ...                    use these rather than orr_e to find rr_r

	   if(nspec(k).gt.0)then
	      if(kflagrev(k).eq.1)then
		 if(itrdep(k).ne.0)then
		   rr_r(k) = rcapb(k)*exp(rexpo(k)*vlntemp-rcape(k)*otc)
		 else
		   rr_r(k) = rcapb(k) 
		 endif
	      else
		 xikv   = 0.0
		 do nk = 1, nkk(k)
		    xikv = xikv + anuk(nk,k)*cgspl(kspk(nk,k))
		 end do
		 xik(k)  = xikv
		 aaf(k)  = exp(xik(k)*otc+anusum(k)*vlntemp) 
	         rr_r(k) = rr_f(k)*aaf(k)
	      endif
	   else
	      rr_r(k)  = 0.0			! reaction irreversible-> rr_r=0
	   endif

c-------------------------------------------------------------------------------
c ... find reactants concentrations product PROD(c_j**nu_jk^prime)_j=1,ns
c ... it is assumed below that nyp is a positive integer (non-zero)
c ... which should be true in general

	   crprod = 1.0

	   do isp = 1,nreac(k)
	      nyp = - nu(isp,k)
	      cvp = cspl(nunk(isp,k))
	      crp1   = cvp**nyp
	      crprod = crprod*crp1
	   end do

c-------------------------------------------------------------------------------
c ... for reverse reaction,
c ... find products concentrations product PROD(c_j**nu_jk^double-prime)_j=1,ns
c ... nspec(k) > 0 if reaction is reversible
c ... it is assumed below that nu(ip3,k) is a positive integer (non-zero)
c ... which should be true in general

	   if(nspec(k).gt.0)then
	      cpprod = 1.0
	      do ip3 = 4,nprod(k)+3
		 nyp = nu(ip3,k)
		 cvp = cspl(nunk(ip3,k))
		 cpp1   = cvp**nyp
		 cpprod = cpprod * cpp1
	      end do
	   else
	      cpprod = 0.0
	   endif

c-------------------------------------------------------------------------------
c ... evaluate reaction rate-of-progress

	   ropl(k) = rr_f(k)*crprod - rr_r(k)*cpprod

c-------------------------------------------------------------------------------

	end do

c ... end of loop over all reactions for rate-of-progress
c-------------------------------------------------------------------------------
c-------------------------------------------------------------------------------
c ... initialize qfac array, product of concentrations to be used in 
c ... third body and fall off reactions

	do k=1,numreac	
	   qfac(k) = 1.0
	end do

c ... evaluate third-body-scaled reaction rates-of-progress
c ... nb. qfac is the product of concentrations, with possible 3rd
c ... body efficiencies.
c ... if the reaction is also a fall-off reaction with kfal=0, then
c ... qfac is just used to find pr below, and then reset to 1.0
c ... otherwise, it is actually used to scale the rate of progress 
c ... later on below

	if(nthb.gt.0)then
	   ctot = 0.0
	   do ksp=1,ns
	      ctot = ctot + cspl(ksp)
	   end do
	   do kthb = 1,nthb
	      k = ithb(kthb)
	      qfac(k) = ctot
	      do ispthb = 1,ntbs(kthb)
		 ksp  = nktb(ispthb,kthb)
		 qfac(k) = qfac(k)+(aik(ispthb,kthb)-1.0d0)*cspl(ksp)
	      end do
	      ropl1(k) = ropl(k)
	      qff(k)   = qfac(k)
              ropl(k)  = qfac(k) * ropl(k)
	   end do
	endif

c-------------------------------------------------------------------------------
c ... evaluate fall-off pressure corrections
c ... vt is dimensional temperature (K), vot = 1/vt 
c ... these are used in the computation of the dimensionless F

	if(nfal.gt.0)then

	   vt = tc*tempref
	   vot= 1.0/vt

	   do kfl = 1, nfal
	      k   = ifal(kfl)
              kfc = kfal(kfl)
	      ifp = ifop(kfl)
	      if(kfc.eq.0)then
		 cthb(k)    = qfac(k)
		 qfac(k) = 1.0
                 ropl(k) = ropl(k)/cthb(k)         ! reset ropl(k)
	      else
		 cthb(k)    = cspl(kfc)
	      endif

	      rklow(k) = capbl(k)*exp(fpar(2,kfl)*vlntemp-capel(k)*otc)
	      pr(k)    = rklow(k) * cthb(k) / rr_f(k)
	      prlog(k) = log10(pr(k))

	      if(ifp.eq.1)then	                    ! Lindemann form
	         fc(k) = 1.0
	      else			            ! ifop = 2, 3, or 4
		 if(ifp.eq.2)then                   ! SRI form
		    xp(k)     = 1.0/(1.0+prlog(k)*prlog(k))
		    econ(1,k) = fpar(4,kfl) * exp(-fpar(5,kfl)*vot) 
		    econ(2,k) = exp(-vt/fpar(6,kfl)) 
		    arg(k)    = econ(1,k) + econ(2,k)
		    fc(k) = (arg(k)**xp(k))*fpar(7,kfl)*vt**fpar(8,kfl)
	         else 			            ! 6 or 7-param TROE form 
		    econ(3,k) = (1.0-fpar(4,kfl))*exp(-vt/fpar(5,kfl))
                    econ(4,k) = fpar(4,kfl)*exp(-vt/fpar(6,kfl))
		    fcent(k)  = econ(3,k) + econ(4,k)
		    if(ifp.eq.4)then                ! 7-parameter TROE form
		       econ(5,k) = exp(-fpar(7,kfl)*vot)
		       fcent(k)  = fcent(k) + econ(5,k)
		    endif

		    fclog  = log10(fcent(k))
		    xn     = 0.75 - 1.27*fclog
		    aa(k)  = prlog(k) - (0.4+0.67*fclog)        ! O10 - O11
		    bb(k)  = xn-0.14*aa(k)        	        ! O10 - O11
		    clog(k)  = aa(k)/bb(k)
		    oclog(k) = 1./(1.0 + clog(k)*clog(k))       ! O9
		    flog     = fclog*oclog(k)	      		! O9
		    fc(k)    = exp(flog*dln10) ! mod fr 10.0**flog, HNN 9/23/98 (O9)

		 endif

	      endif

	      pr1 = 1.0+pr(k)
	      prcm = pr(k)/pr1
	      pcor = fc(k)*prcm

	      ropl(k) = ropl(k) * pcor

	   end do 

	endif

c-------------------------------------------------------------------------------
c ... evaluate contributions to reactants and products reaction rates
c ... initialize pdf

	do k=1,numreac   
	   do nk = 1, nkk(k)
	      ksp    = kspk(nk,k)
	      wspl(ksp) = wspl(ksp) + anuk(nk,k)*ropl(k)
	   end do
	end do

c-------------------------------------------------------------------------------
c ... final scaling

	do ksp=1,ns
       	   yvdot(ksp) = wspl(ksp)*damspwt(ksp)
	end do

c    .................................................................
c
c			Density Source Term
c
c    .................................................................

c  ...	sum_i=1^k h_i w_i

	hwt  = 0.
	sumw = 0.
	do k=1,ns
	   hwt  = hwt +  hspl(k) * yvdot(k)
	   sumw = sumw + yvdot(k)*ospwt(k)
	end do

	wbarbiro = wbar/rc

c  ...	find non-dimensional cp

	cp = 0.0
	do k=1,ns
	   cp = cp + yy(k)*cpspl(k)
	end do

	orcp = 1./(rc*cp)

c  ...	source term

	i = ns+2

	yvdot(i) = orcp*rc/tc*hwt 
     +		   + rc*capsp/p0l*( 1. - gcon/(cp*wbar) ) 
     +	           - rc*wbarbiro*sumw				! O40

c    .................................................................
c
c			Reduced State Vector
c
c    .................................................................

c  ...	zdot 

#ifdef ALL_SPECIES
      do i=1,ns
         zdot(i) = yvdot(i)
      end do
      zdot(ns+1)   = yvdot(ns+2)
      zdot(ns+2) = 0.d0 
#else
        if (n_slack_species_index.eq.1) then     
           do i=1,ns-1
              zdot(i) = yvdot(i+1)
           end do

        elseif (n_slack_species_index.eq.ns) then
           do i=1,ns-1
              zdot(i) = yvdot(i)
           end do

        else
           do i=1,(n_slack_species_index-1)
              zdot(i) = yvdot(i)
           end do
 
           do i=n_slack_species_index,ns-1
              zdot(i) = yvdot(i+1)
           end do
 
        endif          
                   
	zdot(ns)   = yvdot(ns+2)
	zdot(ns+1) = 0.d0 
#endif

	return
	end