module MO_EXP_PROD_LOSS

      CONTAINS

      subroutine EXP_PROD_LOSS( prod, loss, y, rxt, het_rates )

      use CHEM_MODS, only : clscnt1, rxntot, hetcnt
      use MO_GRID,   only : pcnstm1, plnplv

      implicit none

!--------------------------------------------------------------------
!     ... Dummy args                                                                      
!--------------------------------------------------------------------
      real, dimension(plnplv,clscnt1), intent(out) :: &
            prod, &
            loss
      real, intent(in)    ::  y(plnplv,pcnstm1)
      real, intent(in)    ::  rxt(plnplv,rxntot)
      real, intent(in)    ::  het_rates(plnplv,hetcnt)


!--------------------------------------------------------------------
!       ... Loss and production for Explicit method
!--------------------------------------------------------------------

      loss(:,1) = (rxt(:,56)* y(:,1) +rxt(:,55)* y(:,15))* y(:,10)
      prod(:,1) = 0.
      loss(:,2) = ((rxt(:,39) +rxt(:,40))* y(:,1) + rxt(:,4))* y(:,2)
      prod(:,2) = 0.
      loss(:,3) = (rxt(:,64)* y(:,15))* y(:,14)
      prod(:,3) = 0.
      loss(:,4) = ( + rxt(:,168))* y(:,54)
      prod(:,4) = 0.
      loss(:,5) = ( + het_rates(:,15))* y(:,55)
      prod(:,5) =rxt(:,168)*y(:,54)
      loss(:,6) = (rxt(:,65)* y(:,1) +rxt(:,74)* y(:,15))* y(:,61)
      prod(:,6) =.050*rxt(:,56)*y(:,10)*y(:,1)
      loss(:,7) = 0.
      prod(:,7) = 0.
      loss(:,8) = 0.
      prod(:,8) = 0.

      end subroutine EXP_PROD_LOSS

      end module MO_EXP_PROD_LOSS

      module MO_IMP_PROD_LOSS

      CONTAINS

      subroutine IMP_PROD_LOSS( prod, loss, y, rxt, het_rates )

      use CHEM_MODS, only : clscnt4, rxntot, hetcnt, clsze
      use MO_GRID,   only : pcnstm1, plnplv

      implicit none

!--------------------------------------------------------------------
!     ... Dummy args                                                                      
!--------------------------------------------------------------------
      real, dimension(clsze,clscnt4), intent(out) :: &
            prod, &
            loss
      real, intent(in)    ::  y(clsze,pcnstm1)
      real, intent(in)    ::  rxt(clsze,rxntot)
      real, intent(in)    ::  het_rates(clsze,hetcnt)


!--------------------------------------------------------------------
!     ... Local variables                                                                 
!--------------------------------------------------------------------
      integer :: k

!--------------------------------------------------------------------
!       ... Loss and production for Implicit method
!--------------------------------------------------------------------

      do k = 1,clsze
         loss(k,54) = (2.*rxt(k,35)* y(k,1) + (rxt(k,39) +rxt(k,40))* y(k,2) &
                  +rxt(k,42)* y(k,4) + (rxt(k,43) +rxt(k,44))* y(k,5) +rxt(k,56) &
                 * y(k,10) + (rxt(k,66) +rxt(k,68))* y(k,15) + (rxt(k,67) +rxt(k,69)) &
                 * y(k,16) +rxt(k,76)* y(k,18) +rxt(k,90)* y(k,19) +rxt(k,101) &
                 * y(k,28) +rxt(k,110)* y(k,32) +rxt(k,112)* y(k,33) +rxt(k,129) &
                 * y(k,39) +rxt(k,65)* y(k,61) + rxt(k,38))* y(k,1)
         prod(k,54) = (.100*rxt(k,90)*y(k,19) +.200*rxt(k,110)*y(k,32) + &
                 .200*rxt(k,112)*y(k,33) +.765*rxt(k,129)*y(k,39))*y(k,1) &
                  + (rxt(k,133) +rxt(k,134)*y(k,4))*y(k,3) + (.300*rxt(k,85)*y(k,23) + &
                 .300*rxt(k,122)*y(k,36))*y(k,16) +rxt(k,5)*y(k,4) +rxt(k,6)*y(k,5) &
                  +.890*rxt(k,9)*y(k,6) +rxt(k,73)*y(k,15)*y(k,15)
         loss(k,1) = (rxt(k,134)* y(k,4) + rxt(k,133))* y(k,3)
         prod(k,1) =rxt(k,5)*y(k,4)
         loss(k,55) = (rxt(k,42)* y(k,1) +rxt(k,134)* y(k,3) +rxt(k,51)* y(k,6) &
                  +rxt(k,57)* y(k,11) +rxt(k,41)* y(k,16) +rxt(k,78)* y(k,20) &
                  +rxt(k,83)* y(k,23) +rxt(k,104)* y(k,31) + (rxt(k,113) +rxt(k,114)) &
                 * y(k,34) +rxt(k,120)* y(k,36) +rxt(k,92)* y(k,37) +rxt(k,136) &
                 * y(k,41) +rxt(k,98)* y(k,49) +rxt(k,141)* y(k,52) +rxt(k,146) &
                 * y(k,56) +rxt(k,154)* y(k,58) + rxt(k,5))* y(k,4)
         prod(k,55) = (2.000*rxt(k,39)*y(k,2) +rxt(k,43)*y(k,5))*y(k,1) +rxt(k,133) &
                 *y(k,3) +rxt(k,6)*y(k,5) +.110*rxt(k,9)*y(k,6)
         loss(k,48) = ((rxt(k,43) +rxt(k,44))* y(k,1) +rxt(k,46)* y(k,6) +rxt(k,49) &
                 * y(k,15) +rxt(k,52)* y(k,16) +rxt(k,84)* y(k,23) +rxt(k,126) &
                 * y(k,36) + rxt(k,6))* y(k,5)
         prod(k,48) = (rxt(k,41)*y(k,16) +rxt(k,42)*y(k,1) +2.000*rxt(k,51)*y(k,6) + &
                 rxt(k,57)*y(k,11) +rxt(k,78)*y(k,20) +rxt(k,83)*y(k,23) + &
                 rxt(k,92)*y(k,37) +rxt(k,98)*y(k,49) +.920*rxt(k,104)*y(k,31) + &
                 rxt(k,113)*y(k,34) +rxt(k,120)*y(k,36) +rxt(k,136)*y(k,41) + &
                 rxt(k,141)*y(k,52) +1.206*rxt(k,146)*y(k,56) +rxt(k,154)*y(k,58)) &
                 *y(k,4) + (.890*rxt(k,9) +rxt(k,45)*y(k,16) +rxt(k,105)*y(k,31) + &
                 rxt(k,115)*y(k,34) +rxt(k,121)*y(k,36) +rxt(k,130)*y(k,39) + &
                 1.206*rxt(k,147)*y(k,56) +rxt(k,152)*y(k,57) +rxt(k,155)*y(k,58)) &
                 *y(k,6) + (rxt(k,10) +rxt(k,54) +rxt(k,53)*y(k,15))*y(k,8) &
                  + (rxt(k,7) +rxt(k,47))*y(k,9) + (rxt(k,144)*y(k,26) + &
                 rxt(k,151)*y(k,57))*y(k,15) + (.600*rxt(k,19) +rxt(k,88))*y(k,25) &
                  + (rxt(k,20) +rxt(k,127))*y(k,30) +rxt(k,8)*y(k,7) &
                  +.206*rxt(k,148)*y(k,56)*y(k,16) +rxt(k,30)*y(k,57)
         loss(k,52) = (rxt(k,51)* y(k,4) +rxt(k,46)* y(k,5) +rxt(k,62)* y(k,13) &
                  +rxt(k,45)* y(k,16) +rxt(k,77)* y(k,18) +rxt(k,145)* y(k,19) &
                  +rxt(k,82)* y(k,21) +rxt(k,105)* y(k,31) +rxt(k,115)* y(k,34) &
                  +rxt(k,121)* y(k,36) +rxt(k,130)* y(k,39) +rxt(k,150)* y(k,53) &
                  +rxt(k,147)* y(k,56) +rxt(k,152)* y(k,57) +rxt(k,155)* y(k,58) &
                  + rxt(k,9) + rxt(k,132))* y(k,6)
         prod(k,52) = (rxt(k,50)*y(k,7) +.500*rxt(k,164)*y(k,30) +rxt(k,165)*y(k,25)) &
                 *y(k,15) + (rxt(k,7) +rxt(k,47))*y(k,9) +rxt(k,44)*y(k,5)*y(k,1) &
                  +.400*rxt(k,19)*y(k,25)
         loss(k,24) = (rxt(k,50)* y(k,15) + rxt(k,8) + het_rates(k,2))* y(k,7)
         prod(k,24) = (rxt(k,132) +rxt(k,62)*y(k,13) +rxt(k,82)*y(k,21) + &
                 rxt(k,150)*y(k,53))*y(k,6) + (2.000*rxt(k,48) +2.000*rxt(k,131)) &
                 *y(k,9) +rxt(k,49)*y(k,15)*y(k,5)
         loss(k,11) = (rxt(k,53)* y(k,15) + rxt(k,10) + rxt(k,54) + het_rates(k,7)) &
                 * y(k,8)
         prod(k,11) =rxt(k,52)*y(k,16)*y(k,5)
         loss(k,6) = ( + rxt(k,7) + rxt(k,47) + rxt(k,48) + rxt(k,131))* y(k,9)
         prod(k,6) =rxt(k,46)*y(k,6)*y(k,5)
         loss(k,50) = (rxt(k,57)* y(k,4) + 2.*(rxt(k,58) +rxt(k,59))* y(k,11) &
                  +rxt(k,60)* y(k,16) +rxt(k,86)* y(k,23) +rxt(k,107)* y(k,31) &
                  +rxt(k,117)* y(k,34) +rxt(k,123)* y(k,36) +rxt(k,94)* y(k,37) &
                  +rxt(k,138)* y(k,41) +rxt(k,157)* y(k,58))* y(k,11)
         prod(k,50) = (rxt(k,83)*y(k,4) +.900*rxt(k,86)*y(k,11) + &
                 2.000*rxt(k,89)*y(k,23) +rxt(k,108)*y(k,31) +rxt(k,118)*y(k,34) + &
                 rxt(k,124)*y(k,36) +rxt(k,158)*y(k,58))*y(k,23) &
                  + (.750*rxt(k,56)*y(k,10) +.310*rxt(k,76)*y(k,18))*y(k,1) &
                  + (rxt(k,55)*y(k,10) +.700*rxt(k,61)*y(k,12))*y(k,15) +rxt(k,16) &
                 *y(k,21) +rxt(k,18)*y(k,24) +.400*rxt(k,19)*y(k,25) +.300*rxt(k,23) &
                 *y(k,32) +rxt(k,27)*y(k,43)
         loss(k,17) = (rxt(k,61)* y(k,15) + rxt(k,11) + het_rates(k,4))* y(k,12)
         prod(k,17) =rxt(k,60)*y(k,16)*y(k,11)
         loss(k,40) = (rxt(k,62)* y(k,6) +rxt(k,63)* y(k,15) + rxt(k,12) + rxt(k,13) &
                  + het_rates(k,3))* y(k,13)
         prod(k,40) = (rxt(k,57)*y(k,4) +2.000*rxt(k,58)*y(k,11) +rxt(k,59)*y(k,11) + &
                 rxt(k,86)*y(k,23) +.700*rxt(k,94)*y(k,37) +1.200*rxt(k,107)*y(k,31) + &
                 .880*rxt(k,117)*y(k,34) +2.000*rxt(k,123)*y(k,36) + &
                 rxt(k,138)*y(k,41) +.700*rxt(k,157)*y(k,58))*y(k,11) &
                  + (.250*rxt(k,56)*y(k,10) +.540*rxt(k,76)*y(k,18) + &
                 .600*rxt(k,90)*y(k,19) +rxt(k,101)*y(k,28) +.800*rxt(k,110)*y(k,32) + &
                 .700*rxt(k,112)*y(k,33) +1.326*rxt(k,129)*y(k,39))*y(k,1) &
                  + (.300*rxt(k,61)*y(k,12) +.500*rxt(k,87)*y(k,24) + &
                 .500*rxt(k,97)*y(k,28) +.500*rxt(k,151)*y(k,57) +rxt(k,162)*y(k,45) + &
                 .500*rxt(k,164)*y(k,30) +rxt(k,165)*y(k,25))*y(k,15) &
                  + (rxt(k,78)*y(k,20) +.510*rxt(k,104)*y(k,31) + &
                 .250*rxt(k,113)*y(k,34) +rxt(k,120)*y(k,36) +rxt(k,141)*y(k,52) + &
                 .072*rxt(k,146)*y(k,56))*y(k,4) + (.600*rxt(k,105)*y(k,31) + &
                 .250*rxt(k,115)*y(k,34) +rxt(k,121)*y(k,36) +.072*rxt(k,147)*y(k,56)) &
                 *y(k,6) + (.600*rxt(k,108)*y(k,31) +.250*rxt(k,118)*y(k,34) + &
                 rxt(k,124)*y(k,36))*y(k,23) +rxt(k,11)*y(k,12) &
                  +.008*rxt(k,148)*y(k,56)*y(k,16) +rxt(k,17)*y(k,22) +1.340*rxt(k,21) &
                 *y(k,33) +2.000*rxt(k,125)*y(k,36)*y(k,36) +rxt(k,26)*y(k,44) &
                  +rxt(k,33)*y(k,47) +rxt(k,32)*y(k,48) +2.000*rxt(k,100)*y(k,50) &
                  +rxt(k,30)*y(k,57) +.690*rxt(k,31)*y(k,60)
         loss(k,47) = ((rxt(k,66) +rxt(k,68))* y(k,1) +rxt(k,49)* y(k,5) +rxt(k,50) &
                 * y(k,7) +rxt(k,53)* y(k,8) +rxt(k,55)* y(k,10) +rxt(k,61)* y(k,12) &
                  +rxt(k,63)* y(k,13) +rxt(k,64)* y(k,14) + 2.*rxt(k,73)* y(k,15) &
                  +rxt(k,72)* y(k,16) +rxt(k,71)* y(k,17) +rxt(k,75)* y(k,18) &
                  +rxt(k,102)* y(k,19) +rxt(k,81)* y(k,21) +rxt(k,80)* y(k,22) &
                  +rxt(k,87)* y(k,24) +rxt(k,165)* y(k,25) +rxt(k,144)* y(k,26) &
                  +rxt(k,91)* y(k,27) +rxt(k,97)* y(k,28) +rxt(k,103)* y(k,29) &
                  +rxt(k,164)* y(k,30) +rxt(k,109)* y(k,32) +rxt(k,111)* y(k,33) &
                  +rxt(k,119)* y(k,35) +rxt(k,96)* y(k,38) +rxt(k,128)* y(k,39) &
                  +rxt(k,135)* y(k,40) +rxt(k,139)* y(k,42) +rxt(k,140)* y(k,43) &
                  +rxt(k,143)* y(k,44) +rxt(k,162)* y(k,45) +rxt(k,163)* y(k,46) &
                  +rxt(k,167)* y(k,47) +rxt(k,166)* y(k,48) +rxt(k,153)* y(k,51) &
                  +rxt(k,149)* y(k,53) +rxt(k,151)* y(k,57) +rxt(k,159)* y(k,59) &
                  +rxt(k,161)* y(k,60) +rxt(k,74)* y(k,61))* y(k,15)
         prod(k,47) = (2.000*rxt(k,38) +.750*rxt(k,56)*y(k,10) +rxt(k,65)*y(k,61) + &
                 rxt(k,67)*y(k,16) +rxt(k,69)*y(k,16) +.330*rxt(k,76)*y(k,18) + &
                 .270*rxt(k,90)*y(k,19) +.120*rxt(k,101)*y(k,28) + &
                 .080*rxt(k,110)*y(k,32) +.215*rxt(k,112)*y(k,33) + &
                 1.156*rxt(k,129)*y(k,39))*y(k,1) + (.300*rxt(k,61)*y(k,12) + &
                 .500*rxt(k,80)*y(k,22) +.500*rxt(k,96)*y(k,38) + &
                 .100*rxt(k,119)*y(k,35))*y(k,15) + (rxt(k,41)*y(k,4) + &
                 rxt(k,45)*y(k,6))*y(k,16) +rxt(k,8)*y(k,7) +rxt(k,11)*y(k,12) &
                  +2.000*rxt(k,15)*y(k,17) +rxt(k,17)*y(k,22) +rxt(k,18)*y(k,24) &
                  +.660*rxt(k,22)*y(k,33) +rxt(k,24)*y(k,38) +rxt(k,25)*y(k,42) &
                  +rxt(k,26)*y(k,44) +rxt(k,29)*y(k,59)
         loss(k,51) = ((rxt(k,67) +rxt(k,69))* y(k,1) +rxt(k,41)* y(k,4) +rxt(k,52) &
                 * y(k,5) +rxt(k,45)* y(k,6) +rxt(k,60)* y(k,11) +rxt(k,72)* y(k,15) &
                  + 2.*rxt(k,70)* y(k,16) +rxt(k,79)* y(k,20) +rxt(k,85)* y(k,23) &
                  +rxt(k,106)* y(k,31) +rxt(k,116)* y(k,34) +rxt(k,122)* y(k,36) &
                  +rxt(k,93)* y(k,37) +rxt(k,137)* y(k,41) +rxt(k,142)* y(k,52) &
                  +rxt(k,148)* y(k,56) +rxt(k,156)* y(k,58))* y(k,16)
         prod(k,51) = (rxt(k,64)*y(k,14) +rxt(k,74)*y(k,61) +rxt(k,63)*y(k,13) + &
                 rxt(k,66)*y(k,1) +rxt(k,68)*y(k,1) +rxt(k,71)*y(k,17) + &
                 .250*rxt(k,97)*y(k,28) +.200*rxt(k,119)*y(k,35) +rxt(k,151)*y(k,57) + &
                 rxt(k,162)*y(k,45) +rxt(k,163)*y(k,46) +.500*rxt(k,164)*y(k,30) + &
                 rxt(k,166)*y(k,48) +.600*rxt(k,167)*y(k,47))*y(k,15) &
                  + (rxt(k,57)*y(k,4) +2.000*rxt(k,58)*y(k,11) + &
                 .900*rxt(k,86)*y(k,23) +rxt(k,94)*y(k,37) +rxt(k,107)*y(k,31) + &
                 .730*rxt(k,117)*y(k,34) +rxt(k,123)*y(k,36) +rxt(k,138)*y(k,41) + &
                 rxt(k,157)*y(k,58))*y(k,11) + (.400*rxt(k,56)*y(k,10) + &
                 rxt(k,65)*y(k,61) +.190*rxt(k,76)*y(k,18) +.060*rxt(k,90)*y(k,19) + &
                 .120*rxt(k,101)*y(k,28) +.060*rxt(k,110)*y(k,32) + &
                 .275*rxt(k,112)*y(k,33) +.102*rxt(k,129)*y(k,39))*y(k,1) &
                  + (rxt(k,78)*y(k,20) +rxt(k,92)*y(k,37) +rxt(k,104)*y(k,31) + &
                 .470*rxt(k,113)*y(k,34) +rxt(k,136)*y(k,41) + &
                 .794*rxt(k,146)*y(k,56) +1.500*rxt(k,154)*y(k,58))*y(k,4) &
                  + (rxt(k,62)*y(k,13) +rxt(k,105)*y(k,31) +.470*rxt(k,115)*y(k,34) + &
                 .794*rxt(k,147)*y(k,56) +rxt(k,152)*y(k,57) + &
                 1.500*rxt(k,155)*y(k,58))*y(k,6) + (rxt(k,108)*y(k,31) + &
                 .470*rxt(k,118)*y(k,34) +1.500*rxt(k,158)*y(k,58))*y(k,23) &
                  + (rxt(k,10) +rxt(k,54))*y(k,8) + (rxt(k,99) +rxt(k,100))*y(k,50) &
                  +rxt(k,11)*y(k,12) +2.000*rxt(k,12)*y(k,13) +.794*rxt(k,148)*y(k,56) &
                 *y(k,16) +rxt(k,16)*y(k,21) +rxt(k,17)*y(k,22) +1.340*rxt(k,21) &
                 *y(k,33) +1.200*rxt(k,95)*y(k,37)*y(k,37) +rxt(k,24)*y(k,38) &
                  +rxt(k,25)*y(k,42) +2.000*rxt(k,33)*y(k,47) +rxt(k,32)*y(k,48) &
                  +rxt(k,28)*y(k,53) +rxt(k,30)*y(k,57) +rxt(k,31)*y(k,60)
         loss(k,5) = (rxt(k,71)* y(k,15) + rxt(k,15) + het_rates(k,1))* y(k,17)
         prod(k,5) =rxt(k,70)*y(k,16)*y(k,16)
         loss(k,36) = (rxt(k,76)* y(k,1) +rxt(k,77)* y(k,6) +rxt(k,75)* y(k,15)) &
                 * y(k,18)
         prod(k,36) = (.070*rxt(k,90)*y(k,19) +.119*rxt(k,129)*y(k,39))*y(k,1) &
                  +.700*rxt(k,23)*y(k,32)
         loss(k,31) = (rxt(k,90)* y(k,1) +rxt(k,145)* y(k,6) +rxt(k,102)* y(k,15)) &
                 * y(k,19)
         prod(k,31) = 0.
         loss(k,29) = (rxt(k,78)* y(k,4) +rxt(k,79)* y(k,16))* y(k,20)
         prod(k,29) = (rxt(k,75)*y(k,18) +.500*rxt(k,80)*y(k,22))*y(k,15)
         loss(k,37) = (rxt(k,82)* y(k,6) +rxt(k,81)* y(k,15) + rxt(k,16) &
                  + het_rates(k,25))* y(k,21)
         prod(k,37) = (rxt(k,78)*y(k,20) +rxt(k,92)*y(k,37) +.270*rxt(k,136)*y(k,41)) &
                 *y(k,4) + (.500*rxt(k,76)*y(k,18) +.040*rxt(k,110)*y(k,32))*y(k,1) &
                  + (.500*rxt(k,96)*y(k,38) +rxt(k,163)*y(k,46))*y(k,15) &
                  + (.800*rxt(k,94)*y(k,11) +1.600*rxt(k,95)*y(k,37))*y(k,37) &
                  +rxt(k,17)*y(k,22) +rxt(k,24)*y(k,38)
         loss(k,23) = (rxt(k,80)* y(k,15) + rxt(k,17) + het_rates(k,5))* y(k,22)
         prod(k,23) =rxt(k,79)*y(k,20)*y(k,16)
         loss(k,53) = (rxt(k,83)* y(k,4) +rxt(k,84)* y(k,5) +rxt(k,86)* y(k,11) &
                  +rxt(k,85)* y(k,16) + 2.*rxt(k,89)* y(k,23) +rxt(k,108)* y(k,31) &
                  +rxt(k,118)* y(k,34) +rxt(k,158)* y(k,58))* y(k,23)
         prod(k,53) = (rxt(k,82)*y(k,21) +.530*rxt(k,115)*y(k,34) + &
                 rxt(k,121)*y(k,36) +rxt(k,150)*y(k,53))*y(k,6) &
                  + (.530*rxt(k,113)*y(k,34) +rxt(k,120)*y(k,36) +rxt(k,141)*y(k,52)) &
                 *y(k,4) + (rxt(k,81)*y(k,21) +.500*rxt(k,87)*y(k,24) + &
                 rxt(k,149)*y(k,53))*y(k,15) + (.260*rxt(k,117)*y(k,34) + &
                 rxt(k,123)*y(k,36))*y(k,11) + (.600*rxt(k,19) +rxt(k,88))*y(k,25) &
                  +.530*rxt(k,118)*y(k,34)*y(k,23) +.300*rxt(k,23)*y(k,32) &
                  +1.340*rxt(k,21)*y(k,33) +2.000*rxt(k,125)*y(k,36)*y(k,36) &
                  +rxt(k,27)*y(k,43) +rxt(k,26)*y(k,44) +rxt(k,32)*y(k,48) +rxt(k,28) &
                 *y(k,53)
         loss(k,20) = (rxt(k,87)* y(k,15) + rxt(k,18) + het_rates(k,6))* y(k,24)
         prod(k,20) = (.700*rxt(k,85)*y(k,23) +.700*rxt(k,122)*y(k,36))*y(k,16)
         loss(k,21) = (rxt(k,165)* y(k,15) + rxt(k,19) + rxt(k,88))* y(k,25)
         prod(k,21) =rxt(k,84)*y(k,23)*y(k,5)
         loss(k,9) = (rxt(k,144)* y(k,15) + het_rates(k,8))* y(k,26)
         prod(k,9) =rxt(k,77)*y(k,18)*y(k,6)
         loss(k,2) = (rxt(k,91)* y(k,15))* y(k,27)
         prod(k,2) = 0.
         loss(k,12) = (rxt(k,101)* y(k,1) +rxt(k,97)* y(k,15))* y(k,28)
         prod(k,12) = 0.
         loss(k,3) = (rxt(k,103)* y(k,15))* y(k,29)
         prod(k,3) = 0.
         loss(k,25) = (rxt(k,164)* y(k,15) + rxt(k,20) + rxt(k,127))* y(k,30)
         prod(k,25) =rxt(k,126)*y(k,36)*y(k,5)
         loss(k,46) = (rxt(k,104)* y(k,4) +rxt(k,105)* y(k,6) +rxt(k,107)* y(k,11) &
                  +rxt(k,106)* y(k,16) +rxt(k,108)* y(k,23))* y(k,31)
         prod(k,46) = (rxt(k,102)*y(k,19) +1.640*rxt(k,128)*y(k,39) + &
                 .500*rxt(k,161)*y(k,60))*y(k,15) +1.700*rxt(k,130)*y(k,39)*y(k,6)
         loss(k,43) = (rxt(k,110)* y(k,1) +rxt(k,109)* y(k,15) + rxt(k,23) &
                  + het_rates(k,9))* y(k,32)
         prod(k,43) = (.320*rxt(k,104)*y(k,4) +.350*rxt(k,105)*y(k,6) + &
                 .260*rxt(k,107)*y(k,11) +.350*rxt(k,108)*y(k,23))*y(k,31) &
                  + (.039*rxt(k,146)*y(k,4) +.039*rxt(k,147)*y(k,6) + &
                 .039*rxt(k,148)*y(k,16))*y(k,56) + (.200*rxt(k,90)*y(k,19) + &
                 .442*rxt(k,129)*y(k,39))*y(k,1) +.402*rxt(k,31)*y(k,60)
         loss(k,42) = (rxt(k,112)* y(k,1) +rxt(k,111)* y(k,15) + rxt(k,21) + rxt(k,22) &
                  + het_rates(k,10))* y(k,33)
         prod(k,42) = (.230*rxt(k,104)*y(k,4) +.250*rxt(k,105)*y(k,6) + &
                 .190*rxt(k,107)*y(k,11) +.250*rxt(k,108)*y(k,23))*y(k,31) &
                  + (.167*rxt(k,146)*y(k,4) +.167*rxt(k,147)*y(k,6) + &
                 .167*rxt(k,148)*y(k,16))*y(k,56) + (.400*rxt(k,90)*y(k,19) + &
                 1.122*rxt(k,129)*y(k,39))*y(k,1) +.288*rxt(k,31)*y(k,60)
         loss(k,44) = ((rxt(k,113) +rxt(k,114))* y(k,4) +rxt(k,115)* y(k,6) &
                  +rxt(k,117)* y(k,11) +rxt(k,116)* y(k,16) +rxt(k,118)* y(k,23)) &
                 * y(k,34)
         prod(k,44) = (rxt(k,109)*y(k,32) +.500*rxt(k,111)*y(k,33) + &
                 .200*rxt(k,119)*y(k,35))*y(k,15)
         loss(k,13) = (rxt(k,119)* y(k,15) + het_rates(k,16))* y(k,35)
         prod(k,13) =rxt(k,116)*y(k,34)*y(k,16)
         loss(k,49) = (rxt(k,120)* y(k,4) +rxt(k,126)* y(k,5) +rxt(k,121)* y(k,6) &
                  +rxt(k,123)* y(k,11) +rxt(k,122)* y(k,16) +rxt(k,124)* y(k,23) &
                  + 2.*rxt(k,125)* y(k,36))* y(k,36)
         prod(k,49) = (.500*rxt(k,111)*y(k,33) +.500*rxt(k,119)*y(k,35) + &
                 .800*rxt(k,167)*y(k,47))*y(k,15) + (rxt(k,20) +rxt(k,127))*y(k,30) &
                  +.200*rxt(k,90)*y(k,19)*y(k,1) +.660*rxt(k,21)*y(k,33)
         loss(k,28) = (rxt(k,92)* y(k,4) +rxt(k,94)* y(k,11) +rxt(k,93)* y(k,16) &
                  + 2.*rxt(k,95)* y(k,37))* y(k,37)
         prod(k,28) = (rxt(k,91)*y(k,27) +.500*rxt(k,96)*y(k,38))*y(k,15)
         loss(k,14) = (rxt(k,96)* y(k,15) + rxt(k,24) + het_rates(k,11))* y(k,38)
         prod(k,14) =rxt(k,93)*y(k,37)*y(k,16)
         loss(k,34) = (rxt(k,129)* y(k,1) +rxt(k,130)* y(k,6) +rxt(k,128)* y(k,15)) &
                 * y(k,39)
         prod(k,34) = 0.
         loss(k,4) = (rxt(k,135)* y(k,15))* y(k,40)
         prod(k,4) = 0.
         loss(k,32) = (rxt(k,136)* y(k,4) +rxt(k,138)* y(k,11) +rxt(k,137)* y(k,16)) &
                 * y(k,41)
         prod(k,32) = (1.330*rxt(k,103)*y(k,29) +rxt(k,135)*y(k,40) + &
                 rxt(k,139)*y(k,42))*y(k,15)
         loss(k,15) = (rxt(k,139)* y(k,15) + rxt(k,25) + het_rates(k,12))* y(k,42)
         prod(k,15) =rxt(k,137)*y(k,41)*y(k,16)
         loss(k,27) = (rxt(k,140)* y(k,15) + rxt(k,27))* y(k,43)
         prod(k,27) = (.820*rxt(k,136)*y(k,4) +.820*rxt(k,138)*y(k,11))*y(k,41) &
                  +.100*rxt(k,128)*y(k,39)*y(k,15) +.820*rxt(k,25)*y(k,42)
         loss(k,16) = (rxt(k,143)* y(k,15) + rxt(k,26) + het_rates(k,13))* y(k,44)
         prod(k,16) =rxt(k,142)*y(k,52)*y(k,16)
         loss(k,22) = (rxt(k,162)* y(k,15) + het_rates(k,20))* y(k,45)
         prod(k,22) = (rxt(k,59)*y(k,11) +.300*rxt(k,94)*y(k,37) + &
                 .250*rxt(k,107)*y(k,31) +.250*rxt(k,117)*y(k,34) + &
                 .300*rxt(k,157)*y(k,58))*y(k,11)
         loss(k,10) = (rxt(k,163)* y(k,15) + het_rates(k,21))* y(k,46)
         prod(k,10) = (.200*rxt(k,94)*y(k,11) +.400*rxt(k,95)*y(k,37))*y(k,37)
         loss(k,30) = (rxt(k,167)* y(k,15) + rxt(k,33) + het_rates(k,22))* y(k,47)
         prod(k,30) = (.530*rxt(k,113)*y(k,4) +.530*rxt(k,115)*y(k,6) + &
                 .260*rxt(k,117)*y(k,11) +.530*rxt(k,118)*y(k,23))*y(k,34) &
                  + (.250*rxt(k,154)*y(k,4) +.250*rxt(k,155)*y(k,6) + &
                 .100*rxt(k,157)*y(k,11) +.250*rxt(k,158)*y(k,23))*y(k,58) +rxt(k,99) &
                 *y(k,50)
         loss(k,41) = (rxt(k,166)* y(k,15) + rxt(k,32) + het_rates(k,23))* y(k,48)
         prod(k,41) = (.220*rxt(k,113)*y(k,4) +.220*rxt(k,115)*y(k,6) + &
                 .230*rxt(k,117)*y(k,11) +.220*rxt(k,118)*y(k,23))*y(k,34) &
                  + (.250*rxt(k,154)*y(k,4) +.250*rxt(k,155)*y(k,6) + &
                 .100*rxt(k,157)*y(k,11) +.250*rxt(k,158)*y(k,23))*y(k,58) &
                  + (.500*rxt(k,80)*y(k,22) +.500*rxt(k,164)*y(k,30))*y(k,15)
         loss(k,19) = (rxt(k,98)* y(k,4))* y(k,49)
         prod(k,19) =.750*rxt(k,97)*y(k,28)*y(k,15)
         loss(k,7) = ( + rxt(k,99) + rxt(k,100))* y(k,50)
         prod(k,7) =rxt(k,98)*y(k,49)*y(k,4)
         loss(k,18) = (rxt(k,153)* y(k,15) + het_rates(k,24))* y(k,51)
         prod(k,18) = (.370*rxt(k,104)*y(k,4) +.400*rxt(k,105)*y(k,6) + &
                 .300*rxt(k,107)*y(k,11) +.400*rxt(k,108)*y(k,23))*y(k,31) &
                  + (rxt(k,151)*y(k,15) +rxt(k,152)*y(k,6))*y(k,57)
         loss(k,33) = (rxt(k,141)* y(k,4) +rxt(k,142)* y(k,16))* y(k,52)
         prod(k,33) = (rxt(k,140)*y(k,43) +rxt(k,143)*y(k,44))*y(k,15)
         loss(k,39) = (rxt(k,150)* y(k,6) +rxt(k,149)* y(k,15) + rxt(k,28) &
                  + het_rates(k,14))* y(k,53)
         prod(k,39) = (.250*rxt(k,113)*y(k,4) +.250*rxt(k,115)*y(k,6) + &
                 .240*rxt(k,117)*y(k,11) +.250*rxt(k,118)*y(k,23))*y(k,34) &
                  + (.250*rxt(k,154)*y(k,4) +.250*rxt(k,155)*y(k,6) + &
                 .100*rxt(k,157)*y(k,11) +.250*rxt(k,158)*y(k,23))*y(k,58) &
                  + (.950*rxt(k,110)*y(k,32) +.800*rxt(k,112)*y(k,33))*y(k,1) &
                  + (rxt(k,144)*y(k,26) +rxt(k,166)*y(k,48))*y(k,15)
         loss(k,38) = (rxt(k,146)* y(k,4) +rxt(k,147)* y(k,6) +rxt(k,148)* y(k,16) &
                  + het_rates(k,26))* y(k,56)
         prod(k,38) =rxt(k,145)*y(k,19)*y(k,6)
         loss(k,35) = (rxt(k,152)* y(k,6) +rxt(k,151)* y(k,15) + rxt(k,30) &
                  + het_rates(k,18))* y(k,57)
         prod(k,35) = (.080*rxt(k,104)*y(k,31) +rxt(k,114)*y(k,34) + &
                 .794*rxt(k,146)*y(k,56))*y(k,4) + (.794*rxt(k,147)*y(k,6) + &
                 .794*rxt(k,148)*y(k,16))*y(k,56)
         loss(k,45) = (rxt(k,154)* y(k,4) +rxt(k,155)* y(k,6) +rxt(k,157)* y(k,11) &
                  +rxt(k,156)* y(k,16) +rxt(k,158)* y(k,23))* y(k,58)
         prod(k,45) = (rxt(k,153)*y(k,51) +rxt(k,159)*y(k,59) + &
                 .500*rxt(k,161)*y(k,60))*y(k,15)
         loss(k,8) = ((rxt(k,159) +rxt(k,160))* y(k,15) + rxt(k,29) + het_rates(k,17)) &
                 * y(k,59)
         prod(k,8) = (rxt(k,148)*y(k,56) +rxt(k,156)*y(k,58))*y(k,16)
         loss(k,26) = (rxt(k,161)* y(k,15) + rxt(k,31) + het_rates(k,19))* y(k,60)
         prod(k,26) =rxt(k,106)*y(k,31)*y(k,16)
      end do

      end subroutine IMP_PROD_LOSS

      end module MO_IMP_PROD_LOSS

      module MO_RODAS_PROD_LOSS

      CONTAINS

      subroutine RODAS_PROD_LOSS( prod, loss, y, rxt, het_rates )

      use CHEM_MODS, only : clscnt5, rxntot, hetcnt, clsze
      use MO_GRID,   only : pcnstm1, plnplv

      implicit none

!--------------------------------------------------------------------
!     ... Dummy args                                                                      
!--------------------------------------------------------------------
      real, dimension(clsze,clscnt5), intent(out) :: &
            prod, &
            loss
      real, intent(in)    ::  y(clsze,pcnstm1)
      real, intent(in)    ::  rxt(clsze,rxntot)
      real, intent(in)    ::  het_rates(clsze,hetcnt)


      end subroutine RODAS_PROD_LOSS

      end module MO_RODAS_PROD_LOSS

      module MO_INDPRD

      CONTAINS

      subroutine INDPRD( class, prod, y, extfrc, rxt )

      use CHEM_MODS, only : rxntot, extcnt
      use MO_GRID,   only : pcnstm1, plnplv

      implicit none

!--------------------------------------------------------------------
!       ... Dummy arguments
!--------------------------------------------------------------------
      integer, intent(in) :: class
      real, intent(in)    :: y(plnplv,pcnstm1)
      real, intent(in)    :: rxt(plnplv,rxntot)
      real, intent(in)    :: extfrc(plnplv,extcnt)
      real, intent(inout) :: prod(plnplv,*)

!--------------------------------------------------------------------
!       ... "Independent" production for Explicit species
!--------------------------------------------------------------------
      if( class == 1 ) then
         prod(:,1) =.080*rxt(:,76)*y(:,18)*y(:,1) + extfrc(:,3)
                                                                                          
         prod(:,2) = 0.
                                                                                          
         prod(:,3) = (.560*rxt(:,76)*y(:,18) +.300*rxt(:,90)*y(:,19) + &
                 .500*rxt(:,101)*y(:,28) +.050*rxt(:,110)*y(:,32) + &
                 .200*rxt(:,112)*y(:,33) +.323*rxt(:,129)*y(:,39))*y(:,1) &
                  + (rxt(:,62)*y(:,13) +.220*rxt(:,115)*y(:,34) +rxt(:,150)*y(:,53) + &
                 rxt(:,155)*y(:,58))*y(:,6) + (rxt(:,63)*y(:,13) +rxt(:,149)*y(:,53) + &
                 .500*rxt(:,151)*y(:,57) +.400*rxt(:,167)*y(:,47))*y(:,15) &
                  + (.220*rxt(:,113)*y(:,4) +.110*rxt(:,117)*y(:,11) + &
                 .220*rxt(:,118)*y(:,23))*y(:,34) + (rxt(:,154)*y(:,4) + &
                 .400*rxt(:,157)*y(:,11) +rxt(:,158)*y(:,23))*y(:,58) + (rxt(:,12) + &
                 rxt(:,13))*y(:,13) +rxt(:,16)*y(:,21) +.700*rxt(:,23)*y(:,32) &
                  +1.340*rxt(:,22)*y(:,33) +rxt(:,33)*y(:,47) +rxt(:,28)*y(:,53) &
                  +rxt(:,30)*y(:,57) + extfrc(:,2)
                                                                                          
         prod(:,4) = 0.
                                                                                          
         prod(:,5) = 0.
                                                                                          
         prod(:,6) =rxt(:,13)*y(:,13)
                                                                                          
         prod(:,7) = 0.
                                                                                          
         prod(:,8) = 0.
                                                                                          
!--------------------------------------------------------------------
!       ... "Independent" production for Implicit species
!--------------------------------------------------------------------
      else if( class == 4 ) then
         prod(:,54) =rxt(:,4)*y(:,2) +2.000*rxt(:,1)
                                                                                          
         prod(:,1) = 0.
                                                                                          
         prod(:,55) = + extfrc(:,1)
                                                                                          
         prod(:,48) = 0.
                                                                                          
         prod(:,52) = 0.
                                                                                          
         prod(:,24) = 0.
                                                                                          
         prod(:,11) = 0.
                                                                                          
         prod(:,6) = 0.
                                                                                          
         prod(:,50) = 0.
                                                                                          
         prod(:,17) = 0.
                                                                                          
         prod(:,40) = 0.
                                                                                          
         prod(:,47) =rxt(:,14)
                                                                                          
         prod(:,51) =rxt(:,14)
                                                                                          
         prod(:,5) = 0.
                                                                                          
         prod(:,36) = 0.
                                                                                          
         prod(:,31) = 0.
                                                                                          
         prod(:,29) = 0.
                                                                                          
         prod(:,37) = 0.
                                                                                          
         prod(:,23) = 0.
                                                                                          
         prod(:,53) = 0.
                                                                                          
         prod(:,20) = 0.
                                                                                          
         prod(:,21) = 0.
                                                                                          
         prod(:,9) = 0.
                                                                                          
         prod(:,2) = 0.
                                                                                          
         prod(:,12) = 0.
                                                                                          
         prod(:,3) = 0.
                                                                                          
         prod(:,25) = 0.
                                                                                          
         prod(:,46) = 0.
                                                                                          
         prod(:,43) = 0.
                                                                                          
         prod(:,42) = 0.
                                                                                          
         prod(:,44) = 0.
                                                                                          
         prod(:,13) = 0.
                                                                                          
         prod(:,49) = 0.
                                                                                          
         prod(:,28) = 0.
                                                                                          
         prod(:,14) = 0.
                                                                                          
         prod(:,34) = 0.
                                                                                          
         prod(:,4) = 0.
                                                                                          
         prod(:,32) = 0.
                                                                                          
         prod(:,15) = 0.
                                                                                          
         prod(:,27) = 0.
                                                                                          
         prod(:,16) = 0.
                                                                                          
         prod(:,22) = 0.
                                                                                          
         prod(:,10) = 0.
                                                                                          
         prod(:,30) = 0.
                                                                                          
         prod(:,41) = 0.
                                                                                          
         prod(:,19) = 0.
                                                                                          
         prod(:,7) = 0.
                                                                                          
         prod(:,18) = 0.
                                                                                          
         prod(:,33) = 0.
                                                                                          
         prod(:,39) = 0.
                                                                                          
         prod(:,38) = 0.
                                                                                          
         prod(:,35) = 0.
                                                                                          
         prod(:,45) = 0.
                                                                                          
         prod(:,8) = 0.
                                                                                          
         prod(:,26) = 0.
                                                                                          
      end if                                                                              
                                                                                          
      end subroutine INDPRD                                                               
                                                                                          
      end module MO_INDPRD                                                                

      module MO_IMP_LIN_MATRIX

      CONTAINS

      subroutine IMP_LINMAT01( mat, y, rxt, het_rates )
!----------------------------------------------
!       ... Linear Matrix entries for Implicit species
!----------------------------------------------

      use MO_GRID,   only : pcnstm1
      use CHEM_MODS, only : rxntot, hetcnt, imp_nzcnt, clsze

      implicit none

!----------------------------------------------
!       ... Dummy args
!----------------------------------------------
      real, intent(in)    ::  y(clsze,pcnstm1)
      real, intent(in)    ::  rxt(clsze,rxntot)
      real, intent(in)    ::  het_rates(clsze,hetcnt)
      real, intent(inout) ::  mat(clsze,imp_nzcnt)
!----------------------------------------------
!       ... Local variables
!----------------------------------------------
      integer :: k


      do k = 1,clsze

         mat(k,556) = -( rxt(k,38) + rxt(k,39)*y(k,2) + rxt(k,40)*y(k,2) &
                      + rxt(k,56)*y(k,10) + rxt(k,65)*y(k,61) )
         mat(k,586) = rxt(k,5)
         mat(k,407) = rxt(k,6)
         mat(k,510) = .890*rxt(k,9)
         mat(k,2) = rxt(k,133)

         mat(k,1) = -( rxt(k,133) )
         mat(k,558) = rxt(k,5)

         mat(k,587) = -( rxt(k,5) )
         mat(k,408) = rxt(k,6)
         mat(k,511) = .110*rxt(k,9)
         mat(k,3) = rxt(k,133)
         mat(k,557) = 2.000*rxt(k,39)*y(k,2)

         mat(k,401) = -( rxt(k,6) )
         mat(k,18) = rxt(k,7) + rxt(k,47)
         mat(k,105) = rxt(k,8)
         mat(k,504) = .890*rxt(k,9)
         mat(k,37) = rxt(k,10) + rxt(k,54)
         mat(k,88) = .600*rxt(k,19) + rxt(k,88)
         mat(k,111) = rxt(k,20) + rxt(k,127)
         mat(k,206) = rxt(k,30)

         mat(k,508) = -( rxt(k,9) + rxt(k,132) )
         mat(k,19) = rxt(k,7) + rxt(k,47)
         mat(k,90) = .400*rxt(k,19)

         mat(k,103) = -( rxt(k,8) + het_rates(k,2) )
         mat(k,17) = 2.000*rxt(k,48) + 2.000*rxt(k,131)
         mat(k,487) = rxt(k,132)

         mat(k,35) = -( rxt(k,10) + rxt(k,54) + het_rates(k,7) )

         mat(k,16) = -( rxt(k,7) + rxt(k,47) + rxt(k,48) + rxt(k,131) )

         mat(k,227) = rxt(k,16)
         mat(k,83) = rxt(k,18)
         mat(k,89) = .400*rxt(k,19)
         mat(k,284) = .300*rxt(k,23)
         mat(k,126) = rxt(k,27)
         mat(k,387) = rxt(k,55)*y(k,10)
         mat(k,552) = .750*rxt(k,56)*y(k,10)

         mat(k,65) = -( rxt(k,11) + het_rates(k,4) )

         mat(k,251) = -( rxt(k,12) + rxt(k,13) + het_rates(k,3) )
         mat(k,66) = rxt(k,11)
         mat(k,99) = rxt(k,17)
         mat(k,265) = 1.340*rxt(k,21)
         mat(k,62) = rxt(k,26)
         mat(k,203) = rxt(k,30)
         mat(k,116) = .690*rxt(k,31)
         mat(k,257) = rxt(k,32)
         mat(k,149) = rxt(k,33)
         mat(k,22) = 2.000*rxt(k,100)
         mat(k,542) = .250*rxt(k,56)*y(k,10)

         mat(k,384) = -( rxt(k,55)*y(k,10) + rxt(k,64)*y(k,14) + rxt(k,74)*y(k,61) )
         mat(k,104) = rxt(k,8)
         mat(k,67) = rxt(k,11)
         mat(k,14) = 2.000*rxt(k,15)
         mat(k,101) = rxt(k,17)
         mat(k,82) = rxt(k,18)
         mat(k,268) = .660*rxt(k,22)
         mat(k,53) = rxt(k,24)
         mat(k,58) = rxt(k,25)
         mat(k,63) = rxt(k,26)
         mat(k,26) = rxt(k,29)
         mat(k,549) = 2.000*rxt(k,38) + .750*rxt(k,56)*y(k,10) + rxt(k,65)*y(k,61)

         mat(k,38) = rxt(k,10) + rxt(k,54)
         mat(k,69) = rxt(k,11)
         mat(k,254) = 2.000*rxt(k,12)
         mat(k,228) = rxt(k,16)
         mat(k,102) = rxt(k,17)
         mat(k,271) = 1.340*rxt(k,21)
         mat(k,54) = rxt(k,24)
         mat(k,59) = rxt(k,25)
         mat(k,247) = rxt(k,28)
         mat(k,207) = rxt(k,30)
         mat(k,122) = rxt(k,31)
         mat(k,261) = rxt(k,32)
         mat(k,152) = 2.000*rxt(k,33)
         mat(k,23) = rxt(k,99) + rxt(k,100)
         mat(k,553) = .400*rxt(k,56)*y(k,10) + rxt(k,65)*y(k,61)
         mat(k,388) = rxt(k,64)*y(k,14) + rxt(k,74)*y(k,61)

         mat(k,13) = -( rxt(k,15) + het_rates(k,1) )

         mat(k,275) = .700*rxt(k,23)



         mat(k,224) = -( rxt(k,16) + het_rates(k,25) )
         mat(k,98) = rxt(k,17)
         mat(k,52) = rxt(k,24)

         mat(k,96) = -( rxt(k,17) + het_rates(k,5) )

         mat(k,91) = .600*rxt(k,19) + rxt(k,88)
         mat(k,273) = 1.340*rxt(k,21)
         mat(k,287) = .300*rxt(k,23)
         mat(k,64) = rxt(k,26)
         mat(k,127) = rxt(k,27)
         mat(k,249) = rxt(k,28)
         mat(k,263) = rxt(k,32)

         mat(k,80) = -( rxt(k,18) + het_rates(k,6) )

         mat(k,85) = -( rxt(k,19) + rxt(k,88) )

         mat(k,27) = -( het_rates(k,8) )




         mat(k,107) = -( rxt(k,20) + rxt(k,127) )


         mat(k,280) = -( rxt(k,23) + het_rates(k,9) )
         mat(k,118) = .402*rxt(k,31)

         mat(k,266) = -( rxt(k,21) + rxt(k,22) + het_rates(k,10) )
         mat(k,117) = .288*rxt(k,31)


         mat(k,45) = -( het_rates(k,16) )

         mat(k,112) = rxt(k,20) + rxt(k,127)
         mat(k,270) = .660*rxt(k,21)


         mat(k,50) = -( rxt(k,24) + het_rates(k,11) )




         mat(k,55) = -( rxt(k,25) + het_rates(k,12) )

         mat(k,123) = -( rxt(k,27) )
         mat(k,56) = .820*rxt(k,25)

         mat(k,60) = -( rxt(k,26) + het_rates(k,13) )

         mat(k,92) = -( het_rates(k,20) )

         mat(k,31) = -( het_rates(k,21) )

         mat(k,148) = -( rxt(k,33) + het_rates(k,22) )
         mat(k,21) = rxt(k,99)

         mat(k,258) = -( rxt(k,32) + het_rates(k,23) )


         mat(k,20) = -( rxt(k,99) + rxt(k,100) )

         mat(k,70) = -( het_rates(k,24) )


         mat(k,244) = -( rxt(k,28) + het_rates(k,14) )

         mat(k,233) = -( het_rates(k,26) )

         mat(k,202) = -( rxt(k,30) + het_rates(k,18) )


         mat(k,24) = -( rxt(k,29) + het_rates(k,17) )

         mat(k,115) = -( rxt(k,31) + het_rates(k,19) )

      end do

      end subroutine IMP_LINMAT01

      subroutine IMP_LINMAT( mat, y, rxt, het_rates )
!----------------------------------------------
!       ... Linear Matrix entries for Implicit species
!----------------------------------------------

      use MO_GRID,   only : pcnstm1
      use CHEM_MODS, only : rxntot, hetcnt, imp_nzcnt, clsze

      implicit none

!----------------------------------------------
!       ... Dummy args
!----------------------------------------------
      real, intent(in)    ::  y(clsze,pcnstm1)
      real, intent(in)    ::  rxt(clsze,rxntot)
      real, intent(in)    ::  het_rates(clsze,hetcnt)
      real, intent(inout) ::  mat(clsze,imp_nzcnt)

      call IMP_LINMAT01( mat, y, rxt, het_rates )

      end subroutine IMP_LINMAT

      end module MO_IMP_LIN_MATRIX

      module MO_ROD_LIN_MATRIX

      CONTAINS

      subroutine ROD_LINMAT( mat, y, rxt, het_rates )
!----------------------------------------------
!       ... Linear Matrix entries for Implicit species
!----------------------------------------------

      use MO_GRID,   only : pcnstm1
      use CHEM_MODS, only : rxntot, hetcnt, rod_nzcnt, clsze

      implicit none

!----------------------------------------------
!       ... Dummy args
!----------------------------------------------
      real, intent(in)    ::  y(clsze,pcnstm1)
      real, intent(in)    ::  rxt(clsze,rxntot)
      real, intent(in)    ::  het_rates(clsze,hetcnt)
      real, intent(inout) ::  mat(clsze,rod_nzcnt)


      end subroutine ROD_LINMAT

      end module MO_ROD_LIN_MATRIX

      module MO_IMP_NLN_MATRIX

      CONTAINS

      subroutine IMP_NLNMAT01( mat, y, rxt )

      use MO_GRID,   only : pcnstm1
      use CHEM_MODS, only : rxntot, imp_nzcnt, clsze

      implicit none

!----------------------------------------------
!       ... Dummy args
!----------------------------------------------
      real, intent(in)    ::  y(clsze,pcnstm1)
      real, intent(in)    ::  rxt(clsze,rxntot)
      real, intent(inout) ::  mat(clsze,imp_nzcnt)


!----------------------------------------------
!       ... Local variables
!----------------------------------------------
      integer :: k

!----------------------------------------------
!       ... Nonlinear Matrix entries for Implicit species
!----------------------------------------------

      do k = 1,clsze
         mat(k,556) = -( 4.*rxt(k,35)*y(k,1) + rxt(k,42)*y(k,4) + (rxt(k,43) + rxt(k,44) &
                      ) * y(k,5) + (rxt(k,66) + rxt(k,68)) * y(k,15) + (rxt(k,67) &
                      + rxt(k,69)) * y(k,16) + rxt(k,76)*y(k,18) + rxt(k,90)*y(k,19) &
                      + rxt(k,101)*y(k,28) + rxt(k,110)*y(k,32) + rxt(k,112)*y(k,33) &
                      + rxt(k,129)*y(k,39) )
         mat(k,586) = - rxt(k,42)*y(k,1 )
         mat(k,407) = - (rxt(k,43) + rxt(k,44)) * y(k,1)
         mat(k,391) = - (rxt(k,66) + rxt(k,68)) * y(k,1)
         mat(k,482) = - (rxt(k,67) + rxt(k,69)) * y(k,1)
         mat(k,221) = - rxt(k,76)*y(k,1 )
         mat(k,164) = - rxt(k,90)*y(k,1 )
         mat(k,44) = - rxt(k,101)*y(k,1  )
         mat(k,288) = - rxt(k,110)*y(k,1  )
         mat(k,274) = - rxt(k,112)*y(k,1  )
         mat(k,199) = - rxt(k,129)*y(k,1  )

         mat(k,391) = mat(k,391) + 2.000*rxt(k,73)*y(k,15)
         mat(k,482) = mat(k,482) + .300*rxt(k,85)*y(k,23)
         mat(k,531) = mat(k,531) + .300*rxt(k,85)*y(k,16)
         mat(k,556) = mat(k,556) + .100*rxt(k,90)*y(k,19)
         mat(k,164) = mat(k,164) + .100*rxt(k,90)*y(k,1)
         mat(k,556) = mat(k,556) + .200*rxt(k,110)*y(k,32)
         mat(k,288) = mat(k,288) + .200*rxt(k,110)*y(k,1)
         mat(k,556) = mat(k,556) + .200*rxt(k,112)*y(k,33)
         mat(k,274) = mat(k,274) + .200*rxt(k,112)*y(k,1)
         mat(k,482) = mat(k,482) + .300*rxt(k,122)*y(k,36)
         mat(k,420) = mat(k,420) + .300*rxt(k,122)*y(k,16)
         mat(k,556) = mat(k,556) + .765*rxt(k,129)*y(k,39)
         mat(k,199) = mat(k,199) + .765*rxt(k,129)*y(k,1)
         mat(k,2) = mat(k,2) + rxt(k,134)*y(k,4)
         mat(k,586) = mat(k,586) + rxt(k,134)*y(k,3)

         mat(k,1) = -( rxt(k,134)*y(k,4) )
         mat(k,558) = - rxt(k,134)*y(k,3  )


         mat(k,587) = -( rxt(k,41)*y(k,16) + rxt(k,42)*y(k,1) + rxt(k,51)*y(k,6) &
                      + rxt(k,57)*y(k,11) + rxt(k,78)*y(k,20) + rxt(k,83)*y(k,23) &
                      + rxt(k,92)*y(k,37) + rxt(k,98)*y(k,49) + rxt(k,104)*y(k,31) &
                      + (rxt(k,113) + rxt(k,114)) * y(k,34) + rxt(k,120)*y(k,36) &
                      + rxt(k,134)*y(k,3) + rxt(k,136)*y(k,41) + rxt(k,141)*y(k,52) &
                      + rxt(k,146)*y(k,56) + rxt(k,154)*y(k,58) )
         mat(k,483) = - rxt(k,41)*y(k,4 )
         mat(k,557) = - rxt(k,42)*y(k,4 )
         mat(k,511) = - rxt(k,51)*y(k,4 )
         mat(k,448) = - rxt(k,57)*y(k,4 )
         mat(k,147) = - rxt(k,78)*y(k,4 )
         mat(k,532) = - rxt(k,83)*y(k,4 )
         mat(k,138) = - rxt(k,92)*y(k,4 )
         mat(k,79) = - rxt(k,98)*y(k,4 )
         mat(k,340) = - rxt(k,104)*y(k,4  )
         mat(k,306) = - (rxt(k,113) + rxt(k,114)) * y(k,4)
         mat(k,421) = - rxt(k,120)*y(k,4  )
         mat(k,3) = - rxt(k,134)*y(k,4  )
         mat(k,176) = - rxt(k,136)*y(k,4  )
         mat(k,184) = - rxt(k,141)*y(k,4  )
         mat(k,242) = - rxt(k,146)*y(k,4  )
         mat(k,321) = - rxt(k,154)*y(k,4  )

         mat(k,557) = mat(k,557) + rxt(k,43)*y(k,5)
         mat(k,408) = mat(k,408) + rxt(k,43)*y(k,1)

         mat(k,401) = -( (rxt(k,43) + rxt(k,44)) * y(k,1) + rxt(k,46)*y(k,6) + rxt(k,49) &
                      *y(k,15) + rxt(k,52)*y(k,16) + rxt(k,84)*y(k,23) + rxt(k,126) &
                      *y(k,36) )
         mat(k,550) = - (rxt(k,43) + rxt(k,44)) * y(k,5)
         mat(k,504) = - rxt(k,46)*y(k,5 )
         mat(k,385) = - rxt(k,49)*y(k,5 )
         mat(k,476) = - rxt(k,52)*y(k,5 )
         mat(k,525) = - rxt(k,84)*y(k,5 )
         mat(k,414) = - rxt(k,126)*y(k,5  )

         mat(k,580) = mat(k,580) + rxt(k,41)*y(k,16)
         mat(k,476) = mat(k,476) + rxt(k,41)*y(k,4)
         mat(k,550) = mat(k,550) + rxt(k,42)*y(k,4)
         mat(k,580) = mat(k,580) + rxt(k,42)*y(k,1)
         mat(k,504) = mat(k,504) + rxt(k,45)*y(k,16)
         mat(k,476) = mat(k,476) + rxt(k,45)*y(k,6)
         mat(k,580) = mat(k,580) + 2.000*rxt(k,51)*y(k,6)
         mat(k,504) = mat(k,504) + 2.000*rxt(k,51)*y(k,4)
         mat(k,37) = mat(k,37) + rxt(k,53)*y(k,15)
         mat(k,385) = mat(k,385) + rxt(k,53)*y(k,8)
         mat(k,580) = mat(k,580) + rxt(k,57)*y(k,11)
         mat(k,441) = mat(k,441) + rxt(k,57)*y(k,4)
         mat(k,580) = mat(k,580) + rxt(k,78)*y(k,20)
         mat(k,145) = mat(k,145) + rxt(k,78)*y(k,4)
         mat(k,580) = mat(k,580) + rxt(k,83)*y(k,23)
         mat(k,525) = mat(k,525) + rxt(k,83)*y(k,4)
         mat(k,580) = mat(k,580) + rxt(k,92)*y(k,37)
         mat(k,135) = mat(k,135) + rxt(k,92)*y(k,4)
         mat(k,580) = mat(k,580) + rxt(k,98)*y(k,49)
         mat(k,77) = mat(k,77) + rxt(k,98)*y(k,4)
         mat(k,580) = mat(k,580) + .920*rxt(k,104)*y(k,31)
         mat(k,333) = mat(k,333) + .920*rxt(k,104)*y(k,4)
         mat(k,504) = mat(k,504) + rxt(k,105)*y(k,31)
         mat(k,333) = mat(k,333) + rxt(k,105)*y(k,6)
         mat(k,580) = mat(k,580) + rxt(k,113)*y(k,34)
         mat(k,300) = mat(k,300) + rxt(k,113)*y(k,4)
         mat(k,504) = mat(k,504) + rxt(k,115)*y(k,34)
         mat(k,300) = mat(k,300) + rxt(k,115)*y(k,6)
         mat(k,580) = mat(k,580) + rxt(k,120)*y(k,36)
         mat(k,414) = mat(k,414) + rxt(k,120)*y(k,4)
         mat(k,504) = mat(k,504) + rxt(k,121)*y(k,36)
         mat(k,414) = mat(k,414) + rxt(k,121)*y(k,6)
         mat(k,504) = mat(k,504) + rxt(k,130)*y(k,39)
         mat(k,194) = mat(k,194) + rxt(k,130)*y(k,6)
         mat(k,580) = mat(k,580) + rxt(k,136)*y(k,41)
         mat(k,172) = mat(k,172) + rxt(k,136)*y(k,4)
         mat(k,580) = mat(k,580) + rxt(k,141)*y(k,52)
         mat(k,181) = mat(k,181) + rxt(k,141)*y(k,4)
         mat(k,385) = mat(k,385) + rxt(k,144)*y(k,26)
         mat(k,30) = mat(k,30) + rxt(k,144)*y(k,15)
         mat(k,580) = mat(k,580) + 1.206*rxt(k,146)*y(k,56)
         mat(k,239) = mat(k,239) + 1.206*rxt(k,146)*y(k,4)
         mat(k,504) = mat(k,504) + 1.206*rxt(k,147)*y(k,56)
         mat(k,239) = mat(k,239) + 1.206*rxt(k,147)*y(k,6)
         mat(k,476) = mat(k,476) + .206*rxt(k,148)*y(k,56)
         mat(k,239) = mat(k,239) + .206*rxt(k,148)*y(k,16)
         mat(k,385) = mat(k,385) + rxt(k,151)*y(k,57)
         mat(k,206) = mat(k,206) + rxt(k,151)*y(k,15)
         mat(k,504) = mat(k,504) + rxt(k,152)*y(k,57)
         mat(k,206) = mat(k,206) + rxt(k,152)*y(k,6)
         mat(k,580) = mat(k,580) + rxt(k,154)*y(k,58)
         mat(k,315) = mat(k,315) + rxt(k,154)*y(k,4)
         mat(k,504) = mat(k,504) + rxt(k,155)*y(k,58)
         mat(k,315) = mat(k,315) + rxt(k,155)*y(k,6)

         mat(k,508) = -( rxt(k,45)*y(k,16) + rxt(k,46)*y(k,5) + rxt(k,51)*y(k,4) &
                      + rxt(k,62)*y(k,13) + rxt(k,77)*y(k,18) + rxt(k,82)*y(k,21) &
                      + rxt(k,105)*y(k,31) + rxt(k,115)*y(k,34) + rxt(k,121)*y(k,36) &
                      + rxt(k,130)*y(k,39) + rxt(k,145)*y(k,19) + rxt(k,147)*y(k,56) &
                      + rxt(k,150)*y(k,53) + rxt(k,152)*y(k,57) + rxt(k,155)*y(k,58) )
         mat(k,480) = - rxt(k,45)*y(k,6 )
         mat(k,405) = - rxt(k,46)*y(k,6 )
         mat(k,584) = - rxt(k,51)*y(k,6 )
         mat(k,255) = - rxt(k,62)*y(k,6 )
         mat(k,220) = - rxt(k,77)*y(k,6 )
         mat(k,229) = - rxt(k,82)*y(k,6 )
         mat(k,337) = - rxt(k,105)*y(k,6  )
         mat(k,304) = - rxt(k,115)*y(k,6  )
         mat(k,418) = - rxt(k,121)*y(k,6  )
         mat(k,197) = - rxt(k,130)*y(k,6  )
         mat(k,163) = - rxt(k,145)*y(k,6  )
         mat(k,241) = - rxt(k,147)*y(k,6  )
         mat(k,248) = - rxt(k,150)*y(k,6  )
         mat(k,208) = - rxt(k,152)*y(k,6  )
         mat(k,319) = - rxt(k,155)*y(k,6  )

         mat(k,554) = mat(k,554) + rxt(k,44)*y(k,5)
         mat(k,405) = mat(k,405) + rxt(k,44)*y(k,1)
         mat(k,106) = mat(k,106) + rxt(k,50)*y(k,15)
         mat(k,389) = mat(k,389) + rxt(k,50)*y(k,7)
         mat(k,389) = mat(k,389) + .500*rxt(k,164)*y(k,30)
         mat(k,114) = mat(k,114) + .500*rxt(k,164)*y(k,15)
         mat(k,389) = mat(k,389) + rxt(k,165)*y(k,25)
         mat(k,90) = mat(k,90) + rxt(k,165)*y(k,15)

         mat(k,103) = -( rxt(k,50)*y(k,15) )
         mat(k,361) = - rxt(k,50)*y(k,7 )

         mat(k,396) = mat(k,396) + rxt(k,49)*y(k,15)
         mat(k,361) = mat(k,361) + rxt(k,49)*y(k,5)
         mat(k,487) = mat(k,487) + rxt(k,62)*y(k,13)
         mat(k,250) = mat(k,250) + rxt(k,62)*y(k,6)
         mat(k,487) = mat(k,487) + rxt(k,82)*y(k,21)
         mat(k,223) = mat(k,223) + rxt(k,82)*y(k,6)
         mat(k,487) = mat(k,487) + rxt(k,150)*y(k,53)
         mat(k,243) = mat(k,243) + rxt(k,150)*y(k,6)

         mat(k,35) = -( rxt(k,53)*y(k,15) )
         mat(k,348) = - rxt(k,53)*y(k,8 )

         mat(k,394) = mat(k,394) + rxt(k,52)*y(k,16)
         mat(k,451) = mat(k,451) + rxt(k,52)*y(k,5)


         mat(k,393) = mat(k,393) + rxt(k,46)*y(k,6)
         mat(k,484) = mat(k,484) + rxt(k,46)*y(k,5)

         mat(k,443) = -( rxt(k,57)*y(k,4) + (4.*rxt(k,58) + 4.*rxt(k,59)) * y(k,11) &
                      + rxt(k,60)*y(k,16) + rxt(k,86)*y(k,23) + rxt(k,94)*y(k,37) &
                      + rxt(k,107)*y(k,31) + rxt(k,117)*y(k,34) + rxt(k,123)*y(k,36) &
                      + rxt(k,138)*y(k,41) + rxt(k,157)*y(k,58) )
         mat(k,582) = - rxt(k,57)*y(k,11)
         mat(k,478) = - rxt(k,60)*y(k,11)
         mat(k,527) = - rxt(k,86)*y(k,11)
         mat(k,136) = - rxt(k,94)*y(k,11)
         mat(k,335) = - rxt(k,107)*y(k,11 )
         mat(k,302) = - rxt(k,117)*y(k,11 )
         mat(k,416) = - rxt(k,123)*y(k,11 )
         mat(k,173) = - rxt(k,138)*y(k,11 )
         mat(k,317) = - rxt(k,157)*y(k,11 )

         mat(k,68) = mat(k,68) + .700*rxt(k,61)*y(k,15)
         mat(k,387) = mat(k,387) + .700*rxt(k,61)*y(k,12)
         mat(k,552) = mat(k,552) + .310*rxt(k,76)*y(k,18)
         mat(k,218) = mat(k,218) + .310*rxt(k,76)*y(k,1)
         mat(k,582) = mat(k,582) + rxt(k,83)*y(k,23)
         mat(k,527) = mat(k,527) + rxt(k,83)*y(k,4)
         mat(k,443) = mat(k,443) + .900*rxt(k,86)*y(k,23)
         mat(k,527) = mat(k,527) + .900*rxt(k,86)*y(k,11)
         mat(k,527) = mat(k,527) + 4.000*rxt(k,89)*y(k,23)
         mat(k,527) = mat(k,527) + rxt(k,108)*y(k,31)
         mat(k,335) = mat(k,335) + rxt(k,108)*y(k,23)
         mat(k,527) = mat(k,527) + rxt(k,118)*y(k,34)
         mat(k,302) = mat(k,302) + rxt(k,118)*y(k,23)
         mat(k,527) = mat(k,527) + rxt(k,124)*y(k,36)
         mat(k,416) = mat(k,416) + rxt(k,124)*y(k,23)
         mat(k,527) = mat(k,527) + rxt(k,158)*y(k,58)
         mat(k,317) = mat(k,317) + rxt(k,158)*y(k,23)

         mat(k,65) = -( rxt(k,61)*y(k,15) )
         mat(k,354) = - rxt(k,61)*y(k,12)

         mat(k,423) = mat(k,423) + rxt(k,60)*y(k,16)
         mat(k,456) = mat(k,456) + rxt(k,60)*y(k,11)

         mat(k,251) = -( rxt(k,62)*y(k,6) + rxt(k,63)*y(k,15) )
         mat(k,496) = - rxt(k,62)*y(k,13)
         mat(k,377) = - rxt(k,63)*y(k,13)

         mat(k,572) = mat(k,572) + rxt(k,57)*y(k,11)
         mat(k,433) = mat(k,433) + rxt(k,57)*y(k,4)
         mat(k,433) = mat(k,433) + ( 4.000*rxt(k,58) + 2.000*rxt(k,59)  )*y(k,11)
         mat(k,66) = mat(k,66) + .300*rxt(k,61)*y(k,15)
         mat(k,377) = mat(k,377) + .300*rxt(k,61)*y(k,12)
         mat(k,542) = mat(k,542) + .540*rxt(k,76)*y(k,18)
         mat(k,214) = mat(k,214) + .540*rxt(k,76)*y(k,1)
         mat(k,572) = mat(k,572) + rxt(k,78)*y(k,20)
         mat(k,142) = mat(k,142) + rxt(k,78)*y(k,4)
         mat(k,433) = mat(k,433) + rxt(k,86)*y(k,23)
         mat(k,517) = mat(k,517) + rxt(k,86)*y(k,11)
         mat(k,377) = mat(k,377) + .500*rxt(k,87)*y(k,24)
         mat(k,81) = mat(k,81) + .500*rxt(k,87)*y(k,15)
         mat(k,542) = mat(k,542) + .600*rxt(k,90)*y(k,19)
         mat(k,156) = mat(k,156) + .600*rxt(k,90)*y(k,1)
         mat(k,433) = mat(k,433) + .700*rxt(k,94)*y(k,37)
         mat(k,133) = mat(k,133) + .700*rxt(k,94)*y(k,11)
         mat(k,377) = mat(k,377) + .500*rxt(k,97)*y(k,28)
         mat(k,41) = mat(k,41) + .500*rxt(k,97)*y(k,15)
         mat(k,542) = mat(k,542) + rxt(k,101)*y(k,28)
         mat(k,41) = mat(k,41) + rxt(k,101)*y(k,1)
         mat(k,572) = mat(k,572) + .510*rxt(k,104)*y(k,31)
         mat(k,326) = mat(k,326) + .510*rxt(k,104)*y(k,4)
         mat(k,496) = mat(k,496) + .600*rxt(k,105)*y(k,31)
         mat(k,326) = mat(k,326) + .600*rxt(k,105)*y(k,6)
         mat(k,433) = mat(k,433) + 1.200*rxt(k,107)*y(k,31)
         mat(k,326) = mat(k,326) + 1.200*rxt(k,107)*y(k,11)
         mat(k,517) = mat(k,517) + .600*rxt(k,108)*y(k,31)
         mat(k,326) = mat(k,326) + .600*rxt(k,108)*y(k,23)
         mat(k,542) = mat(k,542) + .800*rxt(k,110)*y(k,32)
         mat(k,278) = mat(k,278) + .800*rxt(k,110)*y(k,1)
         mat(k,542) = mat(k,542) + .700*rxt(k,112)*y(k,33)
         mat(k,265) = mat(k,265) + .700*rxt(k,112)*y(k,1)
         mat(k,572) = mat(k,572) + .250*rxt(k,113)*y(k,34)
         mat(k,295) = mat(k,295) + .250*rxt(k,113)*y(k,4)
         mat(k,496) = mat(k,496) + .250*rxt(k,115)*y(k,34)
         mat(k,295) = mat(k,295) + .250*rxt(k,115)*y(k,6)
         mat(k,433) = mat(k,433) + .880*rxt(k,117)*y(k,34)
         mat(k,295) = mat(k,295) + .880*rxt(k,117)*y(k,11)
         mat(k,517) = mat(k,517) + .250*rxt(k,118)*y(k,34)
         mat(k,295) = mat(k,295) + .250*rxt(k,118)*y(k,23)
         mat(k,572) = mat(k,572) + rxt(k,120)*y(k,36)
         mat(k,411) = mat(k,411) + rxt(k,120)*y(k,4)
         mat(k,496) = mat(k,496) + rxt(k,121)*y(k,36)
         mat(k,411) = mat(k,411) + rxt(k,121)*y(k,6)
         mat(k,433) = mat(k,433) + 2.000*rxt(k,123)*y(k,36)
         mat(k,411) = mat(k,411) + 2.000*rxt(k,123)*y(k,11)
         mat(k,517) = mat(k,517) + rxt(k,124)*y(k,36)
         mat(k,411) = mat(k,411) + rxt(k,124)*y(k,23)
         mat(k,411) = mat(k,411) + 4.000*rxt(k,125)*y(k,36)
         mat(k,542) = mat(k,542) + 1.326*rxt(k,129)*y(k,39)
         mat(k,189) = mat(k,189) + 1.326*rxt(k,129)*y(k,1)
         mat(k,433) = mat(k,433) + rxt(k,138)*y(k,41)
         mat(k,170) = mat(k,170) + rxt(k,138)*y(k,11)
         mat(k,572) = mat(k,572) + rxt(k,141)*y(k,52)
         mat(k,179) = mat(k,179) + rxt(k,141)*y(k,4)
         mat(k,572) = mat(k,572) + .072*rxt(k,146)*y(k,56)
         mat(k,234) = mat(k,234) + .072*rxt(k,146)*y(k,4)
         mat(k,496) = mat(k,496) + .072*rxt(k,147)*y(k,56)
         mat(k,234) = mat(k,234) + .072*rxt(k,147)*y(k,6)
         mat(k,468) = mat(k,468) + .008*rxt(k,148)*y(k,56)
         mat(k,234) = mat(k,234) + .008*rxt(k,148)*y(k,16)
         mat(k,377) = mat(k,377) + .500*rxt(k,151)*y(k,57)
         mat(k,203) = mat(k,203) + .500*rxt(k,151)*y(k,15)
         mat(k,433) = mat(k,433) + .700*rxt(k,157)*y(k,58)
         mat(k,311) = mat(k,311) + .700*rxt(k,157)*y(k,11)
         mat(k,377) = mat(k,377) + rxt(k,162)*y(k,45)
         mat(k,93) = mat(k,93) + rxt(k,162)*y(k,15)
         mat(k,377) = mat(k,377) + .500*rxt(k,164)*y(k,30)
         mat(k,108) = mat(k,108) + .500*rxt(k,164)*y(k,15)
         mat(k,377) = mat(k,377) + rxt(k,165)*y(k,25)
         mat(k,86) = mat(k,86) + rxt(k,165)*y(k,15)

      end do

      end subroutine IMP_NLNMAT01

      subroutine IMP_NLNMAT02( mat, y, rxt )

      use MO_GRID,   only : pcnstm1
      use CHEM_MODS, only : rxntot, imp_nzcnt, clsze

      implicit none

!----------------------------------------------
!       ... Dummy args
!----------------------------------------------
      real, intent(in)    ::  y(clsze,pcnstm1)
      real, intent(in)    ::  rxt(clsze,rxntot)
      real, intent(inout) ::  mat(clsze,imp_nzcnt)


!----------------------------------------------
!       ... Local variables
!----------------------------------------------
      integer :: k

!----------------------------------------------
!       ... Nonlinear Matrix entries for Implicit species
!----------------------------------------------

      do k = 1,clsze
         mat(k,384) = -( rxt(k,49)*y(k,5) + rxt(k,50)*y(k,7) + rxt(k,53)*y(k,8) &
                      + rxt(k,61)*y(k,12) + rxt(k,63)*y(k,13) + (rxt(k,66) + rxt(k,68) &
                      ) * y(k,1) + rxt(k,71)*y(k,17) + rxt(k,72)*y(k,16) + 4.*rxt(k,73) &
                      *y(k,15) + rxt(k,75)*y(k,18) + rxt(k,80)*y(k,22) + rxt(k,81) &
                      *y(k,21) + rxt(k,87)*y(k,24) + rxt(k,91)*y(k,27) + rxt(k,96) &
                      *y(k,38) + rxt(k,97)*y(k,28) + rxt(k,102)*y(k,19) + rxt(k,103) &
                      *y(k,29) + rxt(k,109)*y(k,32) + rxt(k,111)*y(k,33) + rxt(k,119) &
                      *y(k,35) + rxt(k,128)*y(k,39) + rxt(k,135)*y(k,40) + rxt(k,139) &
                      *y(k,42) + rxt(k,140)*y(k,43) + rxt(k,143)*y(k,44) + rxt(k,144) &
                      *y(k,26) + rxt(k,149)*y(k,53) + rxt(k,151)*y(k,57) + rxt(k,153) &
                      *y(k,51) + rxt(k,159)*y(k,59) + rxt(k,161)*y(k,60) + rxt(k,162) &
                      *y(k,45) + rxt(k,163)*y(k,46) + rxt(k,164)*y(k,30) + rxt(k,165) &
                      *y(k,25) + rxt(k,166)*y(k,48) + rxt(k,167)*y(k,47) )
         mat(k,400) = - rxt(k,49)*y(k,15)
         mat(k,104) = - rxt(k,50)*y(k,15)
         mat(k,36) = - rxt(k,53)*y(k,15)
         mat(k,67) = - rxt(k,61)*y(k,15)
         mat(k,252) = - rxt(k,63)*y(k,15)
         mat(k,549) = - (rxt(k,66) + rxt(k,68)) * y(k,15)
         mat(k,14) = - rxt(k,71)*y(k,15)
         mat(k,475) = - rxt(k,72)*y(k,15)
         mat(k,216) = - rxt(k,75)*y(k,15)
         mat(k,101) = - rxt(k,80)*y(k,15)
         mat(k,225) = - rxt(k,81)*y(k,15)
         mat(k,82) = - rxt(k,87)*y(k,15)
         mat(k,6) = - rxt(k,91)*y(k,15)
         mat(k,53) = - rxt(k,96)*y(k,15)
         mat(k,42) = - rxt(k,97)*y(k,15)
         mat(k,160) = - rxt(k,102)*y(k,15 )
         mat(k,9) = - rxt(k,103)*y(k,15 )
         mat(k,282) = - rxt(k,109)*y(k,15 )
         mat(k,268) = - rxt(k,111)*y(k,15 )
         mat(k,47) = - rxt(k,119)*y(k,15 )
         mat(k,193) = - rxt(k,128)*y(k,15 )
         mat(k,12) = - rxt(k,135)*y(k,15 )
         mat(k,58) = - rxt(k,139)*y(k,15 )
         mat(k,125) = - rxt(k,140)*y(k,15 )
         mat(k,63) = - rxt(k,143)*y(k,15 )
         mat(k,29) = - rxt(k,144)*y(k,15 )
         mat(k,245) = - rxt(k,149)*y(k,15 )
         mat(k,205) = - rxt(k,151)*y(k,15 )
         mat(k,72) = - rxt(k,153)*y(k,15 )
         mat(k,26) = - rxt(k,159)*y(k,15 )
         mat(k,121) = - rxt(k,161)*y(k,15 )
         mat(k,94) = - rxt(k,162)*y(k,15 )
         mat(k,33) = - rxt(k,163)*y(k,15 )
         mat(k,110) = - rxt(k,164)*y(k,15 )
         mat(k,87) = - rxt(k,165)*y(k,15 )
         mat(k,259) = - rxt(k,166)*y(k,15 )
         mat(k,150) = - rxt(k,167)*y(k,15 )

         mat(k,579) = mat(k,579) + rxt(k,41)*y(k,16)
         mat(k,475) = mat(k,475) + rxt(k,41)*y(k,4)
         mat(k,503) = mat(k,503) + rxt(k,45)*y(k,16)
         mat(k,475) = mat(k,475) + rxt(k,45)*y(k,6)
         mat(k,67) = mat(k,67) + .300*rxt(k,61)*y(k,15)
         mat(k,384) = mat(k,384) + .300*rxt(k,61)*y(k,12)
         mat(k,549) = mat(k,549) + ( rxt(k,67) + rxt(k,69)  )*y(k,16)
         mat(k,475) = mat(k,475) + ( rxt(k,67) + rxt(k,69)  )*y(k,1)
         mat(k,549) = mat(k,549) + .330*rxt(k,76)*y(k,18)
         mat(k,216) = mat(k,216) + .330*rxt(k,76)*y(k,1)
         mat(k,384) = mat(k,384) + .500*rxt(k,80)*y(k,22)
         mat(k,101) = mat(k,101) + .500*rxt(k,80)*y(k,15)
         mat(k,549) = mat(k,549) + .270*rxt(k,90)*y(k,19)
         mat(k,160) = mat(k,160) + .270*rxt(k,90)*y(k,1)
         mat(k,384) = mat(k,384) + .500*rxt(k,96)*y(k,38)
         mat(k,53) = mat(k,53) + .500*rxt(k,96)*y(k,15)
         mat(k,549) = mat(k,549) + .120*rxt(k,101)*y(k,28)
         mat(k,42) = mat(k,42) + .120*rxt(k,101)*y(k,1)
         mat(k,549) = mat(k,549) + .080*rxt(k,110)*y(k,32)
         mat(k,282) = mat(k,282) + .080*rxt(k,110)*y(k,1)
         mat(k,549) = mat(k,549) + .215*rxt(k,112)*y(k,33)
         mat(k,268) = mat(k,268) + .215*rxt(k,112)*y(k,1)
         mat(k,384) = mat(k,384) + .100*rxt(k,119)*y(k,35)
         mat(k,47) = mat(k,47) + .100*rxt(k,119)*y(k,15)
         mat(k,549) = mat(k,549) + 1.156*rxt(k,129)*y(k,39)
         mat(k,193) = mat(k,193) + 1.156*rxt(k,129)*y(k,1)

         mat(k,479) = -( rxt(k,41)*y(k,4) + rxt(k,45)*y(k,6) + rxt(k,52)*y(k,5) &
                      + rxt(k,60)*y(k,11) + (rxt(k,67) + rxt(k,69)) * y(k,1) &
                      + 4.*rxt(k,70)*y(k,16) + rxt(k,72)*y(k,15) + rxt(k,79)*y(k,20) &
                      + rxt(k,85)*y(k,23) + rxt(k,93)*y(k,37) + rxt(k,106)*y(k,31) &
                      + rxt(k,116)*y(k,34) + rxt(k,122)*y(k,36) + rxt(k,137)*y(k,41) &
                      + rxt(k,142)*y(k,52) + rxt(k,148)*y(k,56) + rxt(k,156)*y(k,58) )
         mat(k,583) = - rxt(k,41)*y(k,16)
         mat(k,507) = - rxt(k,45)*y(k,16)
         mat(k,404) = - rxt(k,52)*y(k,16)
         mat(k,444) = - rxt(k,60)*y(k,16)
         mat(k,553) = - (rxt(k,67) + rxt(k,69)) * y(k,16)
         mat(k,388) = - rxt(k,72)*y(k,16)
         mat(k,146) = - rxt(k,79)*y(k,16)
         mat(k,528) = - rxt(k,85)*y(k,16)
         mat(k,137) = - rxt(k,93)*y(k,16)
         mat(k,336) = - rxt(k,106)*y(k,16 )
         mat(k,303) = - rxt(k,116)*y(k,16 )
         mat(k,417) = - rxt(k,122)*y(k,16 )
         mat(k,174) = - rxt(k,137)*y(k,16 )
         mat(k,182) = - rxt(k,142)*y(k,16 )
         mat(k,240) = - rxt(k,148)*y(k,16 )
         mat(k,318) = - rxt(k,156)*y(k,16 )

         mat(k,583) = mat(k,583) + rxt(k,57)*y(k,11)
         mat(k,444) = mat(k,444) + rxt(k,57)*y(k,4)
         mat(k,444) = mat(k,444) + 4.000*rxt(k,58)*y(k,11)
         mat(k,507) = mat(k,507) + rxt(k,62)*y(k,13)
         mat(k,254) = mat(k,254) + rxt(k,62)*y(k,6)
         mat(k,254) = mat(k,254) + rxt(k,63)*y(k,15)
         mat(k,388) = mat(k,388) + rxt(k,63)*y(k,13)
         mat(k,553) = mat(k,553) + ( rxt(k,66) + rxt(k,68)  )*y(k,15)
         mat(k,388) = mat(k,388) + ( rxt(k,66) + rxt(k,68)  )*y(k,1)
         mat(k,388) = mat(k,388) + rxt(k,71)*y(k,17)
         mat(k,15) = mat(k,15) + rxt(k,71)*y(k,15)
         mat(k,553) = mat(k,553) + .190*rxt(k,76)*y(k,18)
         mat(k,219) = mat(k,219) + .190*rxt(k,76)*y(k,1)
         mat(k,583) = mat(k,583) + rxt(k,78)*y(k,20)
         mat(k,146) = mat(k,146) + rxt(k,78)*y(k,4)
         mat(k,444) = mat(k,444) + .900*rxt(k,86)*y(k,23)
         mat(k,528) = mat(k,528) + .900*rxt(k,86)*y(k,11)
         mat(k,553) = mat(k,553) + .060*rxt(k,90)*y(k,19)
         mat(k,162) = mat(k,162) + .060*rxt(k,90)*y(k,1)
         mat(k,583) = mat(k,583) + rxt(k,92)*y(k,37)
         mat(k,137) = mat(k,137) + rxt(k,92)*y(k,4)
         mat(k,444) = mat(k,444) + rxt(k,94)*y(k,37)
         mat(k,137) = mat(k,137) + rxt(k,94)*y(k,11)
         mat(k,137) = mat(k,137) + 2.400*rxt(k,95)*y(k,37)
         mat(k,388) = mat(k,388) + .250*rxt(k,97)*y(k,28)
         mat(k,43) = mat(k,43) + .250*rxt(k,97)*y(k,15)
         mat(k,553) = mat(k,553) + .120*rxt(k,101)*y(k,28)
         mat(k,43) = mat(k,43) + .120*rxt(k,101)*y(k,1)
         mat(k,583) = mat(k,583) + rxt(k,104)*y(k,31)
         mat(k,336) = mat(k,336) + rxt(k,104)*y(k,4)
         mat(k,507) = mat(k,507) + rxt(k,105)*y(k,31)
         mat(k,336) = mat(k,336) + rxt(k,105)*y(k,6)
         mat(k,444) = mat(k,444) + rxt(k,107)*y(k,31)
         mat(k,336) = mat(k,336) + rxt(k,107)*y(k,11)
         mat(k,528) = mat(k,528) + rxt(k,108)*y(k,31)
         mat(k,336) = mat(k,336) + rxt(k,108)*y(k,23)
         mat(k,553) = mat(k,553) + .060*rxt(k,110)*y(k,32)
         mat(k,285) = mat(k,285) + .060*rxt(k,110)*y(k,1)
         mat(k,553) = mat(k,553) + .275*rxt(k,112)*y(k,33)
         mat(k,271) = mat(k,271) + .275*rxt(k,112)*y(k,1)
         mat(k,583) = mat(k,583) + .470*rxt(k,113)*y(k,34)
         mat(k,303) = mat(k,303) + .470*rxt(k,113)*y(k,4)
         mat(k,507) = mat(k,507) + .470*rxt(k,115)*y(k,34)
         mat(k,303) = mat(k,303) + .470*rxt(k,115)*y(k,6)
         mat(k,444) = mat(k,444) + .730*rxt(k,117)*y(k,34)
         mat(k,303) = mat(k,303) + .730*rxt(k,117)*y(k,11)
         mat(k,528) = mat(k,528) + .470*rxt(k,118)*y(k,34)
         mat(k,303) = mat(k,303) + .470*rxt(k,118)*y(k,23)
         mat(k,388) = mat(k,388) + .200*rxt(k,119)*y(k,35)
         mat(k,49) = mat(k,49) + .200*rxt(k,119)*y(k,15)
         mat(k,444) = mat(k,444) + rxt(k,123)*y(k,36)
         mat(k,417) = mat(k,417) + rxt(k,123)*y(k,11)
         mat(k,553) = mat(k,553) + .102*rxt(k,129)*y(k,39)
         mat(k,196) = mat(k,196) + .102*rxt(k,129)*y(k,1)
         mat(k,583) = mat(k,583) + rxt(k,136)*y(k,41)
         mat(k,174) = mat(k,174) + rxt(k,136)*y(k,4)
         mat(k,444) = mat(k,444) + rxt(k,138)*y(k,41)
         mat(k,174) = mat(k,174) + rxt(k,138)*y(k,11)
         mat(k,583) = mat(k,583) + .794*rxt(k,146)*y(k,56)
         mat(k,240) = mat(k,240) + .794*rxt(k,146)*y(k,4)
         mat(k,507) = mat(k,507) + .794*rxt(k,147)*y(k,56)
         mat(k,240) = mat(k,240) + .794*rxt(k,147)*y(k,6)
         mat(k,479) = mat(k,479) + .794*rxt(k,148)*y(k,56)
         mat(k,240) = mat(k,240) + .794*rxt(k,148)*y(k,16)
         mat(k,388) = mat(k,388) + rxt(k,151)*y(k,57)
         mat(k,207) = mat(k,207) + rxt(k,151)*y(k,15)
         mat(k,507) = mat(k,507) + rxt(k,152)*y(k,57)
         mat(k,207) = mat(k,207) + rxt(k,152)*y(k,6)
         mat(k,583) = mat(k,583) + 1.500*rxt(k,154)*y(k,58)
         mat(k,318) = mat(k,318) + 1.500*rxt(k,154)*y(k,4)
         mat(k,507) = mat(k,507) + 1.500*rxt(k,155)*y(k,58)
         mat(k,318) = mat(k,318) + 1.500*rxt(k,155)*y(k,6)
         mat(k,444) = mat(k,444) + rxt(k,157)*y(k,58)
         mat(k,318) = mat(k,318) + rxt(k,157)*y(k,11)
         mat(k,528) = mat(k,528) + 1.500*rxt(k,158)*y(k,58)
         mat(k,318) = mat(k,318) + 1.500*rxt(k,158)*y(k,23)
         mat(k,388) = mat(k,388) + rxt(k,162)*y(k,45)
         mat(k,95) = mat(k,95) + rxt(k,162)*y(k,15)
         mat(k,388) = mat(k,388) + rxt(k,163)*y(k,46)
         mat(k,34) = mat(k,34) + rxt(k,163)*y(k,15)
         mat(k,388) = mat(k,388) + .500*rxt(k,164)*y(k,30)
         mat(k,113) = mat(k,113) + .500*rxt(k,164)*y(k,15)
         mat(k,388) = mat(k,388) + rxt(k,166)*y(k,48)
         mat(k,261) = mat(k,261) + rxt(k,166)*y(k,15)
         mat(k,388) = mat(k,388) + .600*rxt(k,167)*y(k,47)
         mat(k,152) = mat(k,152) + .600*rxt(k,167)*y(k,15)

         mat(k,13) = -( rxt(k,71)*y(k,15) )
         mat(k,344) = - rxt(k,71)*y(k,17)

         mat(k,449) = mat(k,449) + 2.000*rxt(k,70)*y(k,16)

         mat(k,211) = -( rxt(k,75)*y(k,15) + rxt(k,76)*y(k,1) + rxt(k,77)*y(k,6) )
         mat(k,373) = - rxt(k,75)*y(k,18)
         mat(k,538) = - rxt(k,76)*y(k,18)
         mat(k,492) = - rxt(k,77)*y(k,18)

         mat(k,538) = mat(k,538) + .070*rxt(k,90)*y(k,19)
         mat(k,154) = mat(k,154) + .070*rxt(k,90)*y(k,1)
         mat(k,538) = mat(k,538) + .119*rxt(k,129)*y(k,39)
         mat(k,188) = mat(k,188) + .119*rxt(k,129)*y(k,1)

         mat(k,153) = -( rxt(k,90)*y(k,1) + rxt(k,102)*y(k,15) + rxt(k,145)*y(k,6) )
         mat(k,536) = - rxt(k,90)*y(k,19)
         mat(k,368) = - rxt(k,102)*y(k,19 )
         mat(k,489) = - rxt(k,145)*y(k,19 )


         mat(k,140) = -( rxt(k,78)*y(k,4) + rxt(k,79)*y(k,16) )
         mat(k,564) = - rxt(k,78)*y(k,20)
         mat(k,462) = - rxt(k,79)*y(k,20)

         mat(k,366) = mat(k,366) + rxt(k,75)*y(k,18)
         mat(k,210) = mat(k,210) + rxt(k,75)*y(k,15)
         mat(k,366) = mat(k,366) + .500*rxt(k,80)*y(k,22)
         mat(k,97) = mat(k,97) + .500*rxt(k,80)*y(k,15)

      end do

      end subroutine IMP_NLNMAT02

      subroutine IMP_NLNMAT03( mat, y, rxt )

      use MO_GRID,   only : pcnstm1
      use CHEM_MODS, only : rxntot, imp_nzcnt, clsze

      implicit none

!----------------------------------------------
!       ... Dummy args
!----------------------------------------------
      real, intent(in)    ::  y(clsze,pcnstm1)
      real, intent(in)    ::  rxt(clsze,rxntot)
      real, intent(inout) ::  mat(clsze,imp_nzcnt)


!----------------------------------------------
!       ... Local variables
!----------------------------------------------
      integer :: k

!----------------------------------------------
!       ... Nonlinear Matrix entries for Implicit species
!----------------------------------------------

      do k = 1,clsze
         mat(k,224) = -( rxt(k,81)*y(k,15) + rxt(k,82)*y(k,6) )
         mat(k,374) = - rxt(k,81)*y(k,21)
         mat(k,493) = - rxt(k,82)*y(k,21)

         mat(k,539) = mat(k,539) + .500*rxt(k,76)*y(k,18)
         mat(k,212) = mat(k,212) + .500*rxt(k,76)*y(k,1)
         mat(k,569) = mat(k,569) + rxt(k,78)*y(k,20)
         mat(k,141) = mat(k,141) + rxt(k,78)*y(k,4)
         mat(k,569) = mat(k,569) + rxt(k,92)*y(k,37)
         mat(k,132) = mat(k,132) + rxt(k,92)*y(k,4)
         mat(k,431) = mat(k,431) + .800*rxt(k,94)*y(k,37)
         mat(k,132) = mat(k,132) + .800*rxt(k,94)*y(k,11)
         mat(k,132) = mat(k,132) + 3.200*rxt(k,95)*y(k,37)
         mat(k,374) = mat(k,374) + .500*rxt(k,96)*y(k,38)
         mat(k,52) = mat(k,52) + .500*rxt(k,96)*y(k,15)
         mat(k,539) = mat(k,539) + .040*rxt(k,110)*y(k,32)
         mat(k,276) = mat(k,276) + .040*rxt(k,110)*y(k,1)
         mat(k,569) = mat(k,569) + .270*rxt(k,136)*y(k,41)
         mat(k,169) = mat(k,169) + .270*rxt(k,136)*y(k,4)
         mat(k,374) = mat(k,374) + rxt(k,163)*y(k,46)
         mat(k,32) = mat(k,32) + rxt(k,163)*y(k,15)

         mat(k,96) = -( rxt(k,80)*y(k,15) )
         mat(k,360) = - rxt(k,80)*y(k,22)

         mat(k,458) = mat(k,458) + rxt(k,79)*y(k,20)
         mat(k,139) = mat(k,139) + rxt(k,79)*y(k,16)

         mat(k,530) = -( rxt(k,83)*y(k,4) + rxt(k,84)*y(k,5) + rxt(k,85)*y(k,16) &
                      + rxt(k,86)*y(k,11) + 4.*rxt(k,89)*y(k,23) + rxt(k,108)*y(k,31) &
                      + rxt(k,118)*y(k,34) + rxt(k,158)*y(k,58) )
         mat(k,585) = - rxt(k,83)*y(k,23)
         mat(k,406) = - rxt(k,84)*y(k,23)
         mat(k,481) = - rxt(k,85)*y(k,23)
         mat(k,446) = - rxt(k,86)*y(k,23)
         mat(k,338) = - rxt(k,108)*y(k,23 )
         mat(k,305) = - rxt(k,118)*y(k,23 )
         mat(k,320) = - rxt(k,158)*y(k,23 )

         mat(k,390) = mat(k,390) + rxt(k,81)*y(k,21)
         mat(k,230) = mat(k,230) + rxt(k,81)*y(k,15)
         mat(k,509) = mat(k,509) + rxt(k,82)*y(k,21)
         mat(k,230) = mat(k,230) + rxt(k,82)*y(k,6)
         mat(k,390) = mat(k,390) + .500*rxt(k,87)*y(k,24)
         mat(k,84) = mat(k,84) + .500*rxt(k,87)*y(k,15)
         mat(k,585) = mat(k,585) + .530*rxt(k,113)*y(k,34)
         mat(k,305) = mat(k,305) + .530*rxt(k,113)*y(k,4)
         mat(k,509) = mat(k,509) + .530*rxt(k,115)*y(k,34)
         mat(k,305) = mat(k,305) + .530*rxt(k,115)*y(k,6)
         mat(k,446) = mat(k,446) + .260*rxt(k,117)*y(k,34)
         mat(k,305) = mat(k,305) + .260*rxt(k,117)*y(k,11)
         mat(k,530) = mat(k,530) + .530*rxt(k,118)*y(k,34)
         mat(k,305) = mat(k,305) + .530*rxt(k,118)*y(k,23)
         mat(k,585) = mat(k,585) + rxt(k,120)*y(k,36)
         mat(k,419) = mat(k,419) + rxt(k,120)*y(k,4)
         mat(k,509) = mat(k,509) + rxt(k,121)*y(k,36)
         mat(k,419) = mat(k,419) + rxt(k,121)*y(k,6)
         mat(k,446) = mat(k,446) + rxt(k,123)*y(k,36)
         mat(k,419) = mat(k,419) + rxt(k,123)*y(k,11)
         mat(k,419) = mat(k,419) + 4.000*rxt(k,125)*y(k,36)
         mat(k,585) = mat(k,585) + rxt(k,141)*y(k,52)
         mat(k,183) = mat(k,183) + rxt(k,141)*y(k,4)
         mat(k,390) = mat(k,390) + rxt(k,149)*y(k,53)
         mat(k,249) = mat(k,249) + rxt(k,149)*y(k,15)
         mat(k,509) = mat(k,509) + rxt(k,150)*y(k,53)
         mat(k,249) = mat(k,249) + rxt(k,150)*y(k,6)

         mat(k,80) = -( rxt(k,87)*y(k,15) )
         mat(k,357) = - rxt(k,87)*y(k,24)

         mat(k,457) = mat(k,457) + .700*rxt(k,85)*y(k,23)
         mat(k,513) = mat(k,513) + .700*rxt(k,85)*y(k,16)
         mat(k,457) = mat(k,457) + .700*rxt(k,122)*y(k,36)
         mat(k,409) = mat(k,409) + .700*rxt(k,122)*y(k,16)

         mat(k,85) = -( rxt(k,165)*y(k,15) )
         mat(k,358) = - rxt(k,165)*y(k,25 )

         mat(k,395) = mat(k,395) + rxt(k,84)*y(k,23)
         mat(k,514) = mat(k,514) + rxt(k,84)*y(k,5)

         mat(k,27) = -( rxt(k,144)*y(k,15) )
         mat(k,346) = - rxt(k,144)*y(k,26 )

         mat(k,485) = mat(k,485) + rxt(k,77)*y(k,18)
         mat(k,209) = mat(k,209) + rxt(k,77)*y(k,6)

         mat(k,4) = -( rxt(k,91)*y(k,15) )
         mat(k,341) = - rxt(k,91)*y(k,27)


         mat(k,39) = -( rxt(k,97)*y(k,15) + rxt(k,101)*y(k,1) )
         mat(k,349) = - rxt(k,97)*y(k,28)
         mat(k,533) = - rxt(k,101)*y(k,28 )


         mat(k,7) = -( rxt(k,103)*y(k,15) )
         mat(k,342) = - rxt(k,103)*y(k,29 )


         mat(k,107) = -( rxt(k,164)*y(k,15) )
         mat(k,362) = - rxt(k,164)*y(k,30 )

         mat(k,397) = mat(k,397) + rxt(k,126)*y(k,36)
         mat(k,410) = mat(k,410) + rxt(k,126)*y(k,5)

         mat(k,331) = -( rxt(k,104)*y(k,4) + rxt(k,105)*y(k,6) + rxt(k,106)*y(k,16) &
                      + rxt(k,107)*y(k,11) + rxt(k,108)*y(k,23) )
         mat(k,578) = - rxt(k,104)*y(k,31 )
         mat(k,502) = - rxt(k,105)*y(k,31 )
         mat(k,474) = - rxt(k,106)*y(k,31 )
         mat(k,439) = - rxt(k,107)*y(k,31 )
         mat(k,523) = - rxt(k,108)*y(k,31 )

         mat(k,383) = mat(k,383) + rxt(k,102)*y(k,19)
         mat(k,159) = mat(k,159) + rxt(k,102)*y(k,15)
         mat(k,383) = mat(k,383) + 1.640*rxt(k,128)*y(k,39)
         mat(k,192) = mat(k,192) + 1.640*rxt(k,128)*y(k,15)
         mat(k,502) = mat(k,502) + 1.700*rxt(k,130)*y(k,39)
         mat(k,192) = mat(k,192) + 1.700*rxt(k,130)*y(k,6)
         mat(k,383) = mat(k,383) + .500*rxt(k,161)*y(k,60)
         mat(k,120) = mat(k,120) + .500*rxt(k,161)*y(k,15)

         mat(k,280) = -( rxt(k,109)*y(k,15) + rxt(k,110)*y(k,1) )
         mat(k,380) = - rxt(k,109)*y(k,32 )
         mat(k,545) = - rxt(k,110)*y(k,32 )

         mat(k,545) = mat(k,545) + .200*rxt(k,90)*y(k,19)
         mat(k,158) = mat(k,158) + .200*rxt(k,90)*y(k,1)
         mat(k,575) = mat(k,575) + .320*rxt(k,104)*y(k,31)
         mat(k,328) = mat(k,328) + .320*rxt(k,104)*y(k,4)
         mat(k,499) = mat(k,499) + .350*rxt(k,105)*y(k,31)
         mat(k,328) = mat(k,328) + .350*rxt(k,105)*y(k,6)
         mat(k,436) = mat(k,436) + .260*rxt(k,107)*y(k,31)
         mat(k,328) = mat(k,328) + .260*rxt(k,107)*y(k,11)
         mat(k,520) = mat(k,520) + .350*rxt(k,108)*y(k,31)
         mat(k,328) = mat(k,328) + .350*rxt(k,108)*y(k,23)
         mat(k,545) = mat(k,545) + .442*rxt(k,129)*y(k,39)
         mat(k,191) = mat(k,191) + .442*rxt(k,129)*y(k,1)
         mat(k,575) = mat(k,575) + .039*rxt(k,146)*y(k,56)
         mat(k,236) = mat(k,236) + .039*rxt(k,146)*y(k,4)
         mat(k,499) = mat(k,499) + .039*rxt(k,147)*y(k,56)
         mat(k,236) = mat(k,236) + .039*rxt(k,147)*y(k,6)
         mat(k,471) = mat(k,471) + .039*rxt(k,148)*y(k,56)
         mat(k,236) = mat(k,236) + .039*rxt(k,148)*y(k,16)

         mat(k,266) = -( rxt(k,111)*y(k,15) + rxt(k,112)*y(k,1) )
         mat(k,379) = - rxt(k,111)*y(k,33 )
         mat(k,544) = - rxt(k,112)*y(k,33 )

         mat(k,544) = mat(k,544) + .400*rxt(k,90)*y(k,19)
         mat(k,157) = mat(k,157) + .400*rxt(k,90)*y(k,1)
         mat(k,574) = mat(k,574) + .230*rxt(k,104)*y(k,31)
         mat(k,327) = mat(k,327) + .230*rxt(k,104)*y(k,4)
         mat(k,498) = mat(k,498) + .250*rxt(k,105)*y(k,31)
         mat(k,327) = mat(k,327) + .250*rxt(k,105)*y(k,6)
         mat(k,435) = mat(k,435) + .190*rxt(k,107)*y(k,31)
         mat(k,327) = mat(k,327) + .190*rxt(k,107)*y(k,11)
         mat(k,519) = mat(k,519) + .250*rxt(k,108)*y(k,31)
         mat(k,327) = mat(k,327) + .250*rxt(k,108)*y(k,23)
         mat(k,544) = mat(k,544) + 1.122*rxt(k,129)*y(k,39)
         mat(k,190) = mat(k,190) + 1.122*rxt(k,129)*y(k,1)
         mat(k,574) = mat(k,574) + .167*rxt(k,146)*y(k,56)
         mat(k,235) = mat(k,235) + .167*rxt(k,146)*y(k,4)
         mat(k,498) = mat(k,498) + .167*rxt(k,147)*y(k,56)
         mat(k,235) = mat(k,235) + .167*rxt(k,147)*y(k,6)
         mat(k,470) = mat(k,470) + .167*rxt(k,148)*y(k,56)
         mat(k,235) = mat(k,235) + .167*rxt(k,148)*y(k,16)

         mat(k,297) = -( (rxt(k,113) + rxt(k,114)) * y(k,4) + rxt(k,115)*y(k,6) &
                      + rxt(k,116)*y(k,16) + rxt(k,117)*y(k,11) + rxt(k,118)*y(k,23) )
         mat(k,576) = - (rxt(k,113) + rxt(k,114)) * y(k,34)
         mat(k,500) = - rxt(k,115)*y(k,34 )
         mat(k,472) = - rxt(k,116)*y(k,34 )
         mat(k,437) = - rxt(k,117)*y(k,34 )
         mat(k,521) = - rxt(k,118)*y(k,34 )

         mat(k,381) = mat(k,381) + rxt(k,109)*y(k,32)
         mat(k,281) = mat(k,281) + rxt(k,109)*y(k,15)
         mat(k,381) = mat(k,381) + .500*rxt(k,111)*y(k,33)
         mat(k,267) = mat(k,267) + .500*rxt(k,111)*y(k,15)
         mat(k,381) = mat(k,381) + .200*rxt(k,119)*y(k,35)
         mat(k,46) = mat(k,46) + .200*rxt(k,119)*y(k,15)

         mat(k,45) = -( rxt(k,119)*y(k,15) )
         mat(k,350) = - rxt(k,119)*y(k,35 )

         mat(k,452) = mat(k,452) + rxt(k,116)*y(k,34)
         mat(k,290) = mat(k,290) + rxt(k,116)*y(k,16)

         mat(k,415) = -( rxt(k,120)*y(k,4) + rxt(k,121)*y(k,6) + rxt(k,122)*y(k,16) &
                      + rxt(k,123)*y(k,11) + rxt(k,124)*y(k,23) + 4.*rxt(k,125) &
                      *y(k,36) + rxt(k,126)*y(k,5) )
         mat(k,581) = - rxt(k,120)*y(k,36 )
         mat(k,505) = - rxt(k,121)*y(k,36 )
         mat(k,477) = - rxt(k,122)*y(k,36 )
         mat(k,442) = - rxt(k,123)*y(k,36 )
         mat(k,526) = - rxt(k,124)*y(k,36 )
         mat(k,402) = - rxt(k,126)*y(k,36 )

         mat(k,551) = mat(k,551) + .200*rxt(k,90)*y(k,19)
         mat(k,161) = mat(k,161) + .200*rxt(k,90)*y(k,1)
         mat(k,386) = mat(k,386) + .500*rxt(k,111)*y(k,33)
         mat(k,270) = mat(k,270) + .500*rxt(k,111)*y(k,15)
         mat(k,386) = mat(k,386) + .500*rxt(k,119)*y(k,35)
         mat(k,48) = mat(k,48) + .500*rxt(k,119)*y(k,15)
         mat(k,386) = mat(k,386) + .800*rxt(k,167)*y(k,47)
         mat(k,151) = mat(k,151) + .800*rxt(k,167)*y(k,15)

         mat(k,131) = -( rxt(k,92)*y(k,4) + rxt(k,93)*y(k,16) + rxt(k,94)*y(k,11) &
                      + 4.*rxt(k,95)*y(k,37) )
         mat(k,563) = - rxt(k,92)*y(k,37)
         mat(k,461) = - rxt(k,93)*y(k,37)
         mat(k,427) = - rxt(k,94)*y(k,37)

         mat(k,365) = mat(k,365) + rxt(k,91)*y(k,27)
         mat(k,5) = mat(k,5) + rxt(k,91)*y(k,15)
         mat(k,365) = mat(k,365) + .500*rxt(k,96)*y(k,38)
         mat(k,51) = mat(k,51) + .500*rxt(k,96)*y(k,15)

         mat(k,50) = -( rxt(k,96)*y(k,15) )
         mat(k,351) = - rxt(k,96)*y(k,38)

         mat(k,453) = mat(k,453) + rxt(k,93)*y(k,37)
         mat(k,129) = mat(k,129) + rxt(k,93)*y(k,16)

         mat(k,187) = -( rxt(k,128)*y(k,15) + rxt(k,129)*y(k,1) + rxt(k,130)*y(k,6) )
         mat(k,371) = - rxt(k,128)*y(k,39 )
         mat(k,537) = - rxt(k,129)*y(k,39 )
         mat(k,490) = - rxt(k,130)*y(k,39 )


         mat(k,10) = -( rxt(k,135)*y(k,15) )
         mat(k,343) = - rxt(k,135)*y(k,40 )


         mat(k,167) = -( rxt(k,136)*y(k,4) + rxt(k,137)*y(k,16) + rxt(k,138)*y(k,11) )
         mat(k,566) = - rxt(k,136)*y(k,41 )
         mat(k,463) = - rxt(k,137)*y(k,41 )
         mat(k,429) = - rxt(k,138)*y(k,41 )

         mat(k,369) = mat(k,369) + 1.330*rxt(k,103)*y(k,29)
         mat(k,8) = mat(k,8) + 1.330*rxt(k,103)*y(k,15)
         mat(k,369) = mat(k,369) + rxt(k,135)*y(k,40)
         mat(k,11) = mat(k,11) + rxt(k,135)*y(k,15)
         mat(k,369) = mat(k,369) + rxt(k,139)*y(k,42)
         mat(k,57) = mat(k,57) + rxt(k,139)*y(k,15)

      end do

      end subroutine IMP_NLNMAT03

      subroutine IMP_NLNMAT04( mat, y, rxt )

      use MO_GRID,   only : pcnstm1
      use CHEM_MODS, only : rxntot, imp_nzcnt, clsze

      implicit none

!----------------------------------------------
!       ... Dummy args
!----------------------------------------------
      real, intent(in)    ::  y(clsze,pcnstm1)
      real, intent(in)    ::  rxt(clsze,rxntot)
      real, intent(inout) ::  mat(clsze,imp_nzcnt)


!----------------------------------------------
!       ... Local variables
!----------------------------------------------
      integer :: k

!----------------------------------------------
!       ... Nonlinear Matrix entries for Implicit species
!----------------------------------------------

      do k = 1,clsze
         mat(k,55) = -( rxt(k,139)*y(k,15) )
         mat(k,352) = - rxt(k,139)*y(k,42 )

         mat(k,454) = mat(k,454) + rxt(k,137)*y(k,41)
         mat(k,165) = mat(k,165) + rxt(k,137)*y(k,16)

         mat(k,123) = -( rxt(k,140)*y(k,15) )
         mat(k,364) = - rxt(k,140)*y(k,43 )

         mat(k,364) = mat(k,364) + .100*rxt(k,128)*y(k,39)
         mat(k,185) = mat(k,185) + .100*rxt(k,128)*y(k,15)
         mat(k,562) = mat(k,562) + .820*rxt(k,136)*y(k,41)
         mat(k,166) = mat(k,166) + .820*rxt(k,136)*y(k,4)
         mat(k,426) = mat(k,426) + .820*rxt(k,138)*y(k,41)
         mat(k,166) = mat(k,166) + .820*rxt(k,138)*y(k,11)

         mat(k,60) = -( rxt(k,143)*y(k,15) )
         mat(k,353) = - rxt(k,143)*y(k,44 )

         mat(k,455) = mat(k,455) + rxt(k,142)*y(k,52)
         mat(k,177) = mat(k,177) + rxt(k,142)*y(k,16)

         mat(k,92) = -( rxt(k,162)*y(k,15) )
         mat(k,359) = - rxt(k,162)*y(k,45 )

         mat(k,425) = mat(k,425) + 2.000*rxt(k,59)*y(k,11)
         mat(k,425) = mat(k,425) + .300*rxt(k,94)*y(k,37)
         mat(k,130) = mat(k,130) + .300*rxt(k,94)*y(k,11)
         mat(k,425) = mat(k,425) + .250*rxt(k,107)*y(k,31)
         mat(k,323) = mat(k,323) + .250*rxt(k,107)*y(k,11)
         mat(k,425) = mat(k,425) + .250*rxt(k,117)*y(k,34)
         mat(k,291) = mat(k,291) + .250*rxt(k,117)*y(k,11)
         mat(k,425) = mat(k,425) + .300*rxt(k,157)*y(k,58)
         mat(k,308) = mat(k,308) + .300*rxt(k,157)*y(k,11)

         mat(k,31) = -( rxt(k,163)*y(k,15) )
         mat(k,347) = - rxt(k,163)*y(k,46 )

         mat(k,422) = mat(k,422) + .200*rxt(k,94)*y(k,37)
         mat(k,128) = mat(k,128) + .200*rxt(k,94)*y(k,11)
         mat(k,128) = mat(k,128) + .800*rxt(k,95)*y(k,37)

         mat(k,148) = -( rxt(k,167)*y(k,15) )
         mat(k,367) = - rxt(k,167)*y(k,47 )

         mat(k,565) = mat(k,565) + .530*rxt(k,113)*y(k,34)
         mat(k,292) = mat(k,292) + .530*rxt(k,113)*y(k,4)
         mat(k,488) = mat(k,488) + .530*rxt(k,115)*y(k,34)
         mat(k,292) = mat(k,292) + .530*rxt(k,115)*y(k,6)
         mat(k,428) = mat(k,428) + .260*rxt(k,117)*y(k,34)
         mat(k,292) = mat(k,292) + .260*rxt(k,117)*y(k,11)
         mat(k,515) = mat(k,515) + .530*rxt(k,118)*y(k,34)
         mat(k,292) = mat(k,292) + .530*rxt(k,118)*y(k,23)
         mat(k,565) = mat(k,565) + .250*rxt(k,154)*y(k,58)
         mat(k,309) = mat(k,309) + .250*rxt(k,154)*y(k,4)
         mat(k,488) = mat(k,488) + .250*rxt(k,155)*y(k,58)
         mat(k,309) = mat(k,309) + .250*rxt(k,155)*y(k,6)
         mat(k,428) = mat(k,428) + .100*rxt(k,157)*y(k,58)
         mat(k,309) = mat(k,309) + .100*rxt(k,157)*y(k,11)
         mat(k,515) = mat(k,515) + .250*rxt(k,158)*y(k,58)
         mat(k,309) = mat(k,309) + .250*rxt(k,158)*y(k,23)

         mat(k,258) = -( rxt(k,166)*y(k,15) )
         mat(k,378) = - rxt(k,166)*y(k,48 )

         mat(k,378) = mat(k,378) + .500*rxt(k,80)*y(k,22)
         mat(k,100) = mat(k,100) + .500*rxt(k,80)*y(k,15)
         mat(k,573) = mat(k,573) + .220*rxt(k,113)*y(k,34)
         mat(k,296) = mat(k,296) + .220*rxt(k,113)*y(k,4)
         mat(k,497) = mat(k,497) + .220*rxt(k,115)*y(k,34)
         mat(k,296) = mat(k,296) + .220*rxt(k,115)*y(k,6)
         mat(k,434) = mat(k,434) + .230*rxt(k,117)*y(k,34)
         mat(k,296) = mat(k,296) + .230*rxt(k,117)*y(k,11)
         mat(k,518) = mat(k,518) + .220*rxt(k,118)*y(k,34)
         mat(k,296) = mat(k,296) + .220*rxt(k,118)*y(k,23)
         mat(k,573) = mat(k,573) + .250*rxt(k,154)*y(k,58)
         mat(k,312) = mat(k,312) + .250*rxt(k,154)*y(k,4)
         mat(k,497) = mat(k,497) + .250*rxt(k,155)*y(k,58)
         mat(k,312) = mat(k,312) + .250*rxt(k,155)*y(k,6)
         mat(k,434) = mat(k,434) + .100*rxt(k,157)*y(k,58)
         mat(k,312) = mat(k,312) + .100*rxt(k,157)*y(k,11)
         mat(k,518) = mat(k,518) + .250*rxt(k,158)*y(k,58)
         mat(k,312) = mat(k,312) + .250*rxt(k,158)*y(k,23)
         mat(k,378) = mat(k,378) + .500*rxt(k,164)*y(k,30)
         mat(k,109) = mat(k,109) + .500*rxt(k,164)*y(k,15)

         mat(k,74) = -( rxt(k,98)*y(k,4) )
         mat(k,561) = - rxt(k,98)*y(k,49)

         mat(k,356) = mat(k,356) + .750*rxt(k,97)*y(k,28)
         mat(k,40) = mat(k,40) + .750*rxt(k,97)*y(k,15)


         mat(k,559) = mat(k,559) + rxt(k,98)*y(k,49)
         mat(k,73) = mat(k,73) + rxt(k,98)*y(k,4)

         mat(k,70) = -( rxt(k,153)*y(k,15) )
         mat(k,355) = - rxt(k,153)*y(k,51 )

         mat(k,560) = mat(k,560) + .370*rxt(k,104)*y(k,31)
         mat(k,322) = mat(k,322) + .370*rxt(k,104)*y(k,4)
         mat(k,486) = mat(k,486) + .400*rxt(k,105)*y(k,31)
         mat(k,322) = mat(k,322) + .400*rxt(k,105)*y(k,6)
         mat(k,424) = mat(k,424) + .300*rxt(k,107)*y(k,31)
         mat(k,322) = mat(k,322) + .300*rxt(k,107)*y(k,11)
         mat(k,512) = mat(k,512) + .400*rxt(k,108)*y(k,31)
         mat(k,322) = mat(k,322) + .400*rxt(k,108)*y(k,23)
         mat(k,355) = mat(k,355) + rxt(k,151)*y(k,57)
         mat(k,201) = mat(k,201) + rxt(k,151)*y(k,15)
         mat(k,486) = mat(k,486) + rxt(k,152)*y(k,57)
         mat(k,201) = mat(k,201) + rxt(k,152)*y(k,6)

         mat(k,178) = -( rxt(k,141)*y(k,4) + rxt(k,142)*y(k,16) )
         mat(k,567) = - rxt(k,141)*y(k,52 )
         mat(k,464) = - rxt(k,142)*y(k,52 )

         mat(k,370) = mat(k,370) + rxt(k,140)*y(k,43)
         mat(k,124) = mat(k,124) + rxt(k,140)*y(k,15)
         mat(k,370) = mat(k,370) + rxt(k,143)*y(k,44)
         mat(k,61) = mat(k,61) + rxt(k,143)*y(k,15)

         mat(k,244) = -( rxt(k,149)*y(k,15) + rxt(k,150)*y(k,6) )
         mat(k,376) = - rxt(k,149)*y(k,53 )
         mat(k,495) = - rxt(k,150)*y(k,53 )

         mat(k,541) = mat(k,541) + .950*rxt(k,110)*y(k,32)
         mat(k,277) = mat(k,277) + .950*rxt(k,110)*y(k,1)
         mat(k,541) = mat(k,541) + .800*rxt(k,112)*y(k,33)
         mat(k,264) = mat(k,264) + .800*rxt(k,112)*y(k,1)
         mat(k,571) = mat(k,571) + .250*rxt(k,113)*y(k,34)
         mat(k,294) = mat(k,294) + .250*rxt(k,113)*y(k,4)
         mat(k,495) = mat(k,495) + .250*rxt(k,115)*y(k,34)
         mat(k,294) = mat(k,294) + .250*rxt(k,115)*y(k,6)
         mat(k,432) = mat(k,432) + .240*rxt(k,117)*y(k,34)
         mat(k,294) = mat(k,294) + .240*rxt(k,117)*y(k,11)
         mat(k,516) = mat(k,516) + .250*rxt(k,118)*y(k,34)
         mat(k,294) = mat(k,294) + .250*rxt(k,118)*y(k,23)
         mat(k,376) = mat(k,376) + rxt(k,144)*y(k,26)
         mat(k,28) = mat(k,28) + rxt(k,144)*y(k,15)
         mat(k,571) = mat(k,571) + .250*rxt(k,154)*y(k,58)
         mat(k,310) = mat(k,310) + .250*rxt(k,154)*y(k,4)
         mat(k,495) = mat(k,495) + .250*rxt(k,155)*y(k,58)
         mat(k,310) = mat(k,310) + .250*rxt(k,155)*y(k,6)
         mat(k,432) = mat(k,432) + .100*rxt(k,157)*y(k,58)
         mat(k,310) = mat(k,310) + .100*rxt(k,157)*y(k,11)
         mat(k,516) = mat(k,516) + .250*rxt(k,158)*y(k,58)
         mat(k,310) = mat(k,310) + .250*rxt(k,158)*y(k,23)
         mat(k,376) = mat(k,376) + rxt(k,166)*y(k,48)
         mat(k,256) = mat(k,256) + rxt(k,166)*y(k,15)

         mat(k,233) = -( rxt(k,146)*y(k,4) + rxt(k,147)*y(k,6) + rxt(k,148)*y(k,16) )
         mat(k,570) = - rxt(k,146)*y(k,56 )
         mat(k,494) = - rxt(k,147)*y(k,56 )
         mat(k,467) = - rxt(k,148)*y(k,56 )

         mat(k,494) = mat(k,494) + rxt(k,145)*y(k,19)
         mat(k,155) = mat(k,155) + rxt(k,145)*y(k,6)

         mat(k,202) = -( rxt(k,151)*y(k,15) + rxt(k,152)*y(k,6) )
         mat(k,372) = - rxt(k,151)*y(k,57 )
         mat(k,491) = - rxt(k,152)*y(k,57 )

         mat(k,568) = mat(k,568) + .080*rxt(k,104)*y(k,31)
         mat(k,325) = mat(k,325) + .080*rxt(k,104)*y(k,4)
         mat(k,568) = mat(k,568) + rxt(k,114)*y(k,34)
         mat(k,293) = mat(k,293) + rxt(k,114)*y(k,4)
         mat(k,568) = mat(k,568) + .794*rxt(k,146)*y(k,56)
         mat(k,232) = mat(k,232) + .794*rxt(k,146)*y(k,4)
         mat(k,491) = mat(k,491) + .794*rxt(k,147)*y(k,56)
         mat(k,232) = mat(k,232) + .794*rxt(k,147)*y(k,6)
         mat(k,465) = mat(k,465) + .794*rxt(k,148)*y(k,56)
         mat(k,232) = mat(k,232) + .794*rxt(k,148)*y(k,16)

         mat(k,313) = -( rxt(k,154)*y(k,4) + rxt(k,155)*y(k,6) + rxt(k,156)*y(k,16) &
                      + rxt(k,157)*y(k,11) + rxt(k,158)*y(k,23) )
         mat(k,577) = - rxt(k,154)*y(k,58 )
         mat(k,501) = - rxt(k,155)*y(k,58 )
         mat(k,473) = - rxt(k,156)*y(k,58 )
         mat(k,438) = - rxt(k,157)*y(k,58 )
         mat(k,522) = - rxt(k,158)*y(k,58 )

         mat(k,382) = mat(k,382) + rxt(k,153)*y(k,51)
         mat(k,71) = mat(k,71) + rxt(k,153)*y(k,15)
         mat(k,382) = mat(k,382) + rxt(k,159)*y(k,59)
         mat(k,25) = mat(k,25) + rxt(k,159)*y(k,15)
         mat(k,382) = mat(k,382) + .500*rxt(k,161)*y(k,60)
         mat(k,119) = mat(k,119) + .500*rxt(k,161)*y(k,15)

         mat(k,24) = -( (rxt(k,159) + rxt(k,160)) * y(k,15) )
         mat(k,345) = - (rxt(k,159) + rxt(k,160)) * y(k,59)

         mat(k,450) = mat(k,450) + rxt(k,148)*y(k,56)
         mat(k,231) = mat(k,231) + rxt(k,148)*y(k,16)
         mat(k,450) = mat(k,450) + rxt(k,156)*y(k,58)
         mat(k,307) = mat(k,307) + rxt(k,156)*y(k,16)

         mat(k,115) = -( rxt(k,161)*y(k,15) )
         mat(k,363) = - rxt(k,161)*y(k,60 )

         mat(k,459) = mat(k,459) + rxt(k,106)*y(k,31)
         mat(k,324) = mat(k,324) + rxt(k,106)*y(k,16)

      end do

      end subroutine IMP_NLNMAT04

      subroutine IMP_NLNMAT( mat, y, rxt )

      use MO_GRID,   only : pcnstm1
      use CHEM_MODS, only : rxntot, imp_nzcnt, clsze

      implicit none

!----------------------------------------------
!       ... Dummy args
!----------------------------------------------
      real, intent(in)    ::  y(clsze,pcnstm1)
      real, intent(in)    ::  rxt(clsze,rxntot)
      real, intent(inout) ::  mat(clsze,imp_nzcnt)

      call IMP_NLNMAT01( mat, y, rxt )
      call IMP_NLNMAT02( mat, y, rxt )
      call IMP_NLNMAT03( mat, y, rxt )
      call IMP_NLNMAT04( mat, y, rxt )

      end subroutine IMP_NLNMAT

      end module MO_IMP_NLN_MATRIX

      module MO_ROD_NLN_MATRIX

      CONTAINS

      subroutine ROD_NLNMAT( mat, y, rxt )

      use MO_GRID,   only : pcnstm1
      use CHEM_MODS, only : rxntot, rod_nzcnt, clsze

      implicit none

!----------------------------------------------
!       ... Dummy args
!----------------------------------------------
      real, intent(in)    ::  y(clsze,pcnstm1)
      real, intent(in)    ::  rxt(clsze,rxntot)
      real, intent(inout) ::  mat(clsze,rod_nzcnt)


      end subroutine ROD_NLNMAT

      end module MO_ROD_NLN_MATRIX

      module MO_IMP_FACTOR

      CONTAINS
                                                                        
      subroutine IMP_LU_FAC01( lu )
                                                                        
      use CHEM_MODS, only : imp_nzcnt, clsze
                                                                        
      implicit none
                                                                        
!-----------------------------------------------------------------------
!       ... Dummy args
!-----------------------------------------------------------------------
      real, intent(inout) ::   lu(clsze,imp_nzcnt)
                                                                        
!-----------------------------------------------------------------------
!       ... Local variables
!-----------------------------------------------------------------------
      integer :: k
                                                                        
      do k = 1,clsze
         lu(k,1) = 1. / lu(k,1)
         lu(k,2) = lu(k,2) * lu(k,1)
         lu(k,3) = lu(k,3) * lu(k,1)
         lu(k,586) = lu(k,586) - lu(k,2) * lu(k,558)
         lu(k,587) = lu(k,587) - lu(k,3) * lu(k,558)
                                                                        
         lu(k,4) = 1. / lu(k,4)
         lu(k,5) = lu(k,5) * lu(k,4)
         lu(k,6) = lu(k,6) * lu(k,4)
         lu(k,365) = lu(k,365) - lu(k,5) * lu(k,341)
         lu(k,384) = lu(k,384) - lu(k,6) * lu(k,341)
                                                                        
         lu(k,7) = 1. / lu(k,7)
         lu(k,8) = lu(k,8) * lu(k,7)
         lu(k,9) = lu(k,9) * lu(k,7)
         lu(k,369) = lu(k,369) - lu(k,8) * lu(k,342)
         lu(k,384) = lu(k,384) - lu(k,9) * lu(k,342)
                                                                        
         lu(k,10) = 1. / lu(k,10)
         lu(k,11) = lu(k,11) * lu(k,10)
         lu(k,12) = lu(k,12) * lu(k,10)
         lu(k,369) = lu(k,369) - lu(k,11) * lu(k,343)
         lu(k,384) = lu(k,384) - lu(k,12) * lu(k,343)
                                                                        
         lu(k,13) = 1. / lu(k,13)
         lu(k,14) = lu(k,14) * lu(k,13)
         lu(k,15) = lu(k,15) * lu(k,13)
         lu(k,384) = lu(k,384) - lu(k,14) * lu(k,344)
         lu(k,388) = lu(k,388) - lu(k,15) * lu(k,344)
         lu(k,475) = lu(k,475) - lu(k,14) * lu(k,449)
         lu(k,479) = lu(k,479) - lu(k,15) * lu(k,449)
                                                                        
         lu(k,16) = 1. / lu(k,16)
         lu(k,17) = lu(k,17) * lu(k,16)
         lu(k,18) = lu(k,18) * lu(k,16)
         lu(k,19) = lu(k,19) * lu(k,16)
         lu(k,396) = lu(k,396) - lu(k,17) * lu(k,393)
         lu(k,401) = lu(k,401) - lu(k,18) * lu(k,393)
         lu(k,405) = lu(k,405) - lu(k,19) * lu(k,393)
         lu(k,487) = lu(k,487) - lu(k,17) * lu(k,484)
         lu(k,504) = lu(k,504) - lu(k,18) * lu(k,484)
         lu(k,508) = lu(k,508) - lu(k,19) * lu(k,484)
                                                                        
         lu(k,20) = 1. / lu(k,20)
         lu(k,21) = lu(k,21) * lu(k,20)
         lu(k,22) = lu(k,22) * lu(k,20)
         lu(k,23) = lu(k,23) * lu(k,20)
         lu(k,75) = - lu(k,21) * lu(k,73)
         lu(k,76) = - lu(k,22) * lu(k,73)
         lu(k,78) = - lu(k,23) * lu(k,73)
         lu(k,565) = lu(k,565) - lu(k,21) * lu(k,559)
         lu(k,572) = lu(k,572) - lu(k,22) * lu(k,559)
         lu(k,583) = lu(k,583) - lu(k,23) * lu(k,559)
                                                                        
         lu(k,24) = 1. / lu(k,24)
         lu(k,25) = lu(k,25) * lu(k,24)
         lu(k,26) = lu(k,26) * lu(k,24)
         lu(k,237) = - lu(k,25) * lu(k,231)
         lu(k,238) = - lu(k,26) * lu(k,231)
         lu(k,313) = lu(k,313) - lu(k,25) * lu(k,307)
         lu(k,314) = - lu(k,26) * lu(k,307)
         lu(k,382) = lu(k,382) - lu(k,25) * lu(k,345)
         lu(k,384) = lu(k,384) - lu(k,26) * lu(k,345)
         lu(k,473) = lu(k,473) - lu(k,25) * lu(k,450)
         lu(k,475) = lu(k,475) - lu(k,26) * lu(k,450)
                                                                        
         lu(k,27) = 1. / lu(k,27)
         lu(k,28) = lu(k,28) * lu(k,27)
         lu(k,29) = lu(k,29) * lu(k,27)
         lu(k,30) = lu(k,30) * lu(k,27)
         lu(k,213) = - lu(k,28) * lu(k,209)
         lu(k,216) = lu(k,216) - lu(k,29) * lu(k,209)
         lu(k,217) = - lu(k,30) * lu(k,209)
         lu(k,376) = lu(k,376) - lu(k,28) * lu(k,346)
         lu(k,384) = lu(k,384) - lu(k,29) * lu(k,346)
         lu(k,385) = lu(k,385) - lu(k,30) * lu(k,346)
         lu(k,495) = lu(k,495) - lu(k,28) * lu(k,485)
         lu(k,503) = lu(k,503) - lu(k,29) * lu(k,485)
         lu(k,504) = lu(k,504) - lu(k,30) * lu(k,485)
                                                                        
         lu(k,31) = 1. / lu(k,31)
         lu(k,32) = lu(k,32) * lu(k,31)
         lu(k,33) = lu(k,33) * lu(k,31)
         lu(k,34) = lu(k,34) * lu(k,31)
         lu(k,132) = lu(k,132) - lu(k,32) * lu(k,128)
         lu(k,134) = - lu(k,33) * lu(k,128)
         lu(k,137) = lu(k,137) - lu(k,34) * lu(k,128)
         lu(k,374) = lu(k,374) - lu(k,32) * lu(k,347)
         lu(k,384) = lu(k,384) - lu(k,33) * lu(k,347)
         lu(k,388) = lu(k,388) - lu(k,34) * lu(k,347)
         lu(k,431) = lu(k,431) - lu(k,32) * lu(k,422)
         lu(k,440) = - lu(k,33) * lu(k,422)
         lu(k,444) = lu(k,444) - lu(k,34) * lu(k,422)
                                                                        
         lu(k,35) = 1. / lu(k,35)
         lu(k,36) = lu(k,36) * lu(k,35)
         lu(k,37) = lu(k,37) * lu(k,35)
         lu(k,38) = lu(k,38) * lu(k,35)
         lu(k,384) = lu(k,384) - lu(k,36) * lu(k,348)
         lu(k,385) = lu(k,385) - lu(k,37) * lu(k,348)
         lu(k,388) = lu(k,388) - lu(k,38) * lu(k,348)
         lu(k,400) = lu(k,400) - lu(k,36) * lu(k,394)
         lu(k,401) = lu(k,401) - lu(k,37) * lu(k,394)
         lu(k,404) = lu(k,404) - lu(k,38) * lu(k,394)
         lu(k,475) = lu(k,475) - lu(k,36) * lu(k,451)
         lu(k,476) = lu(k,476) - lu(k,37) * lu(k,451)
         lu(k,479) = lu(k,479) - lu(k,38) * lu(k,451)
                                                                        
         lu(k,39) = 1. / lu(k,39)
         lu(k,40) = lu(k,40) * lu(k,39)
         lu(k,41) = lu(k,41) * lu(k,39)
         lu(k,42) = lu(k,42) * lu(k,39)
         lu(k,43) = lu(k,43) * lu(k,39)
         lu(k,44) = lu(k,44) * lu(k,39)
         lu(k,356) = lu(k,356) - lu(k,40) * lu(k,349)
         lu(k,377) = lu(k,377) - lu(k,41) * lu(k,349)
         lu(k,384) = lu(k,384) - lu(k,42) * lu(k,349)
         lu(k,388) = lu(k,388) - lu(k,43) * lu(k,349)
         lu(k,391) = lu(k,391) - lu(k,44) * lu(k,349)
         lu(k,534) = - lu(k,40) * lu(k,533)
         lu(k,542) = lu(k,542) - lu(k,41) * lu(k,533)
         lu(k,549) = lu(k,549) - lu(k,42) * lu(k,533)
         lu(k,553) = lu(k,553) - lu(k,43) * lu(k,533)
         lu(k,556) = lu(k,556) - lu(k,44) * lu(k,533)
                                                                        
      end do
                                                                        
      end subroutine IMP_LU_FAC01
                                                                        
      subroutine IMP_LU_FAC02( lu )
                                                                        
      use CHEM_MODS, only : imp_nzcnt, clsze
                                                                        
      implicit none
                                                                        
!-----------------------------------------------------------------------
!       ... Dummy args
!-----------------------------------------------------------------------
      real, intent(inout) ::   lu(clsze,imp_nzcnt)
                                                                        
!-----------------------------------------------------------------------
!       ... Local variables
!-----------------------------------------------------------------------
      integer :: k
                                                                        
      do k = 1,clsze
         lu(k,45) = 1. / lu(k,45)
         lu(k,46) = lu(k,46) * lu(k,45)
         lu(k,47) = lu(k,47) * lu(k,45)
         lu(k,48) = lu(k,48) * lu(k,45)
         lu(k,49) = lu(k,49) * lu(k,45)
         lu(k,297) = lu(k,297) - lu(k,46) * lu(k,290)
         lu(k,299) = - lu(k,47) * lu(k,290)
         lu(k,301) = - lu(k,48) * lu(k,290)
         lu(k,303) = lu(k,303) - lu(k,49) * lu(k,290)
         lu(k,381) = lu(k,381) - lu(k,46) * lu(k,350)
         lu(k,384) = lu(k,384) - lu(k,47) * lu(k,350)
         lu(k,386) = lu(k,386) - lu(k,48) * lu(k,350)
         lu(k,388) = lu(k,388) - lu(k,49) * lu(k,350)
         lu(k,472) = lu(k,472) - lu(k,46) * lu(k,452)
         lu(k,475) = lu(k,475) - lu(k,47) * lu(k,452)
         lu(k,477) = lu(k,477) - lu(k,48) * lu(k,452)
         lu(k,479) = lu(k,479) - lu(k,49) * lu(k,452)
                                                                        
         lu(k,50) = 1. / lu(k,50)
         lu(k,51) = lu(k,51) * lu(k,50)
         lu(k,52) = lu(k,52) * lu(k,50)
         lu(k,53) = lu(k,53) * lu(k,50)
         lu(k,54) = lu(k,54) * lu(k,50)
         lu(k,131) = lu(k,131) - lu(k,51) * lu(k,129)
         lu(k,132) = lu(k,132) - lu(k,52) * lu(k,129)
         lu(k,134) = lu(k,134) - lu(k,53) * lu(k,129)
         lu(k,137) = lu(k,137) - lu(k,54) * lu(k,129)
         lu(k,365) = lu(k,365) - lu(k,51) * lu(k,351)
         lu(k,374) = lu(k,374) - lu(k,52) * lu(k,351)
         lu(k,384) = lu(k,384) - lu(k,53) * lu(k,351)
         lu(k,388) = lu(k,388) - lu(k,54) * lu(k,351)
         lu(k,461) = lu(k,461) - lu(k,51) * lu(k,453)
         lu(k,466) = - lu(k,52) * lu(k,453)
         lu(k,475) = lu(k,475) - lu(k,53) * lu(k,453)
         lu(k,479) = lu(k,479) - lu(k,54) * lu(k,453)
                                                                        
         lu(k,55) = 1. / lu(k,55)
         lu(k,56) = lu(k,56) * lu(k,55)
         lu(k,57) = lu(k,57) * lu(k,55)
         lu(k,58) = lu(k,58) * lu(k,55)
         lu(k,59) = lu(k,59) * lu(k,55)
         lu(k,166) = lu(k,166) - lu(k,56) * lu(k,165)
         lu(k,167) = lu(k,167) - lu(k,57) * lu(k,165)
         lu(k,171) = - lu(k,58) * lu(k,165)
         lu(k,174) = lu(k,174) - lu(k,59) * lu(k,165)
         lu(k,364) = lu(k,364) - lu(k,56) * lu(k,352)
         lu(k,369) = lu(k,369) - lu(k,57) * lu(k,352)
         lu(k,384) = lu(k,384) - lu(k,58) * lu(k,352)
         lu(k,388) = lu(k,388) - lu(k,59) * lu(k,352)
         lu(k,460) = - lu(k,56) * lu(k,454)
         lu(k,463) = lu(k,463) - lu(k,57) * lu(k,454)
         lu(k,475) = lu(k,475) - lu(k,58) * lu(k,454)
         lu(k,479) = lu(k,479) - lu(k,59) * lu(k,454)
                                                                        
         lu(k,60) = 1. / lu(k,60)
         lu(k,61) = lu(k,61) * lu(k,60)
         lu(k,62) = lu(k,62) * lu(k,60)
         lu(k,63) = lu(k,63) * lu(k,60)
         lu(k,64) = lu(k,64) * lu(k,60)
         lu(k,178) = lu(k,178) - lu(k,61) * lu(k,177)
         lu(k,179) = lu(k,179) - lu(k,62) * lu(k,177)
         lu(k,180) = - lu(k,63) * lu(k,177)
         lu(k,183) = lu(k,183) - lu(k,64) * lu(k,177)
         lu(k,370) = lu(k,370) - lu(k,61) * lu(k,353)
         lu(k,377) = lu(k,377) - lu(k,62) * lu(k,353)
         lu(k,384) = lu(k,384) - lu(k,63) * lu(k,353)
         lu(k,390) = lu(k,390) - lu(k,64) * lu(k,353)
         lu(k,464) = lu(k,464) - lu(k,61) * lu(k,455)
         lu(k,468) = lu(k,468) - lu(k,62) * lu(k,455)
         lu(k,475) = lu(k,475) - lu(k,63) * lu(k,455)
         lu(k,481) = lu(k,481) - lu(k,64) * lu(k,455)
                                                                        
         lu(k,65) = 1. / lu(k,65)
         lu(k,66) = lu(k,66) * lu(k,65)
         lu(k,67) = lu(k,67) * lu(k,65)
         lu(k,68) = lu(k,68) * lu(k,65)
         lu(k,69) = lu(k,69) * lu(k,65)
         lu(k,377) = lu(k,377) - lu(k,66) * lu(k,354)
         lu(k,384) = lu(k,384) - lu(k,67) * lu(k,354)
         lu(k,387) = lu(k,387) - lu(k,68) * lu(k,354)
         lu(k,388) = lu(k,388) - lu(k,69) * lu(k,354)
         lu(k,433) = lu(k,433) - lu(k,66) * lu(k,423)
         lu(k,440) = lu(k,440) - lu(k,67) * lu(k,423)
         lu(k,443) = lu(k,443) - lu(k,68) * lu(k,423)
         lu(k,444) = lu(k,444) - lu(k,69) * lu(k,423)
         lu(k,468) = lu(k,468) - lu(k,66) * lu(k,456)
         lu(k,475) = lu(k,475) - lu(k,67) * lu(k,456)
         lu(k,478) = lu(k,478) - lu(k,68) * lu(k,456)
         lu(k,479) = lu(k,479) - lu(k,69) * lu(k,456)
                                                                        
         lu(k,70) = 1. / lu(k,70)
         lu(k,71) = lu(k,71) * lu(k,70)
         lu(k,72) = lu(k,72) * lu(k,70)
         lu(k,204) = - lu(k,71) * lu(k,201)
         lu(k,205) = lu(k,205) - lu(k,72) * lu(k,201)
         lu(k,330) = - lu(k,71) * lu(k,322)
         lu(k,332) = - lu(k,72) * lu(k,322)
         lu(k,382) = lu(k,382) - lu(k,71) * lu(k,355)
         lu(k,384) = lu(k,384) - lu(k,72) * lu(k,355)
         lu(k,438) = lu(k,438) - lu(k,71) * lu(k,424)
         lu(k,440) = lu(k,440) - lu(k,72) * lu(k,424)
         lu(k,501) = lu(k,501) - lu(k,71) * lu(k,486)
         lu(k,503) = lu(k,503) - lu(k,72) * lu(k,486)
         lu(k,522) = lu(k,522) - lu(k,71) * lu(k,512)
         lu(k,524) = - lu(k,72) * lu(k,512)
         lu(k,577) = lu(k,577) - lu(k,71) * lu(k,560)
         lu(k,579) = lu(k,579) - lu(k,72) * lu(k,560)
                                                                        
         lu(k,74) = 1. / lu(k,74)
         lu(k,75) = lu(k,75) * lu(k,74)
         lu(k,76) = lu(k,76) * lu(k,74)
         lu(k,77) = lu(k,77) * lu(k,74)
         lu(k,78) = lu(k,78) * lu(k,74)
         lu(k,79) = lu(k,79) * lu(k,74)
         lu(k,367) = lu(k,367) - lu(k,75) * lu(k,356)
         lu(k,377) = lu(k,377) - lu(k,76) * lu(k,356)
         lu(k,385) = lu(k,385) - lu(k,77) * lu(k,356)
         lu(k,388) = lu(k,388) - lu(k,78) * lu(k,356)
         lu(k,392) = - lu(k,79) * lu(k,356)
         lu(k,535) = - lu(k,75) * lu(k,534)
         lu(k,542) = lu(k,542) - lu(k,76) * lu(k,534)
         lu(k,550) = lu(k,550) - lu(k,77) * lu(k,534)
         lu(k,553) = lu(k,553) - lu(k,78) * lu(k,534)
         lu(k,557) = lu(k,557) - lu(k,79) * lu(k,534)
         lu(k,565) = lu(k,565) - lu(k,75) * lu(k,561)
         lu(k,572) = lu(k,572) - lu(k,76) * lu(k,561)
         lu(k,580) = lu(k,580) - lu(k,77) * lu(k,561)
         lu(k,583) = lu(k,583) - lu(k,78) * lu(k,561)
         lu(k,587) = lu(k,587) - lu(k,79) * lu(k,561)
                                                                        
         lu(k,80) = 1. / lu(k,80)
         lu(k,81) = lu(k,81) * lu(k,80)
         lu(k,82) = lu(k,82) * lu(k,80)
         lu(k,83) = lu(k,83) * lu(k,80)
         lu(k,84) = lu(k,84) * lu(k,80)
         lu(k,377) = lu(k,377) - lu(k,81) * lu(k,357)
         lu(k,384) = lu(k,384) - lu(k,82) * lu(k,357)
         lu(k,387) = lu(k,387) - lu(k,83) * lu(k,357)
         lu(k,390) = lu(k,390) - lu(k,84) * lu(k,357)
         lu(k,411) = lu(k,411) - lu(k,81) * lu(k,409)
         lu(k,413) = - lu(k,82) * lu(k,409)
         lu(k,416) = lu(k,416) - lu(k,83) * lu(k,409)
         lu(k,419) = lu(k,419) - lu(k,84) * lu(k,409)
         lu(k,468) = lu(k,468) - lu(k,81) * lu(k,457)
         lu(k,475) = lu(k,475) - lu(k,82) * lu(k,457)
         lu(k,478) = lu(k,478) - lu(k,83) * lu(k,457)
         lu(k,481) = lu(k,481) - lu(k,84) * lu(k,457)
         lu(k,517) = lu(k,517) - lu(k,81) * lu(k,513)
         lu(k,524) = lu(k,524) - lu(k,82) * lu(k,513)
         lu(k,527) = lu(k,527) - lu(k,83) * lu(k,513)
         lu(k,530) = lu(k,530) - lu(k,84) * lu(k,513)
                                                                        
      end do
                                                                        
      end subroutine IMP_LU_FAC02
                                                                        
      subroutine IMP_LU_FAC03( lu )
                                                                        
      use CHEM_MODS, only : imp_nzcnt, clsze
                                                                        
      implicit none
                                                                        
!-----------------------------------------------------------------------
!       ... Dummy args
!-----------------------------------------------------------------------
      real, intent(inout) ::   lu(clsze,imp_nzcnt)
                                                                        
!-----------------------------------------------------------------------
!       ... Local variables
!-----------------------------------------------------------------------
      integer :: k
                                                                        
      do k = 1,clsze
         lu(k,85) = 1. / lu(k,85)
         lu(k,86) = lu(k,86) * lu(k,85)
         lu(k,87) = lu(k,87) * lu(k,85)
         lu(k,88) = lu(k,88) * lu(k,85)
         lu(k,89) = lu(k,89) * lu(k,85)
         lu(k,90) = lu(k,90) * lu(k,85)
         lu(k,91) = lu(k,91) * lu(k,85)
         lu(k,377) = lu(k,377) - lu(k,86) * lu(k,358)
         lu(k,384) = lu(k,384) - lu(k,87) * lu(k,358)
         lu(k,385) = lu(k,385) - lu(k,88) * lu(k,358)
         lu(k,387) = lu(k,387) - lu(k,89) * lu(k,358)
         lu(k,389) = lu(k,389) - lu(k,90) * lu(k,358)
         lu(k,390) = lu(k,390) - lu(k,91) * lu(k,358)
         lu(k,398) = - lu(k,86) * lu(k,395)
         lu(k,400) = lu(k,400) - lu(k,87) * lu(k,395)
         lu(k,401) = lu(k,401) - lu(k,88) * lu(k,395)
         lu(k,403) = - lu(k,89) * lu(k,395)
         lu(k,405) = lu(k,405) - lu(k,90) * lu(k,395)
         lu(k,406) = lu(k,406) - lu(k,91) * lu(k,395)
         lu(k,517) = lu(k,517) - lu(k,86) * lu(k,514)
         lu(k,524) = lu(k,524) - lu(k,87) * lu(k,514)
         lu(k,525) = lu(k,525) - lu(k,88) * lu(k,514)
         lu(k,527) = lu(k,527) - lu(k,89) * lu(k,514)
         lu(k,529) = - lu(k,90) * lu(k,514)
         lu(k,530) = lu(k,530) - lu(k,91) * lu(k,514)
                                                                        
         lu(k,92) = 1. / lu(k,92)
         lu(k,93) = lu(k,93) * lu(k,92)
         lu(k,94) = lu(k,94) * lu(k,92)
         lu(k,95) = lu(k,95) * lu(k,92)
         lu(k,133) = lu(k,133) - lu(k,93) * lu(k,130)
         lu(k,134) = lu(k,134) - lu(k,94) * lu(k,130)
         lu(k,137) = lu(k,137) - lu(k,95) * lu(k,130)
         lu(k,295) = lu(k,295) - lu(k,93) * lu(k,291)
         lu(k,299) = lu(k,299) - lu(k,94) * lu(k,291)
         lu(k,303) = lu(k,303) - lu(k,95) * lu(k,291)
         lu(k,311) = lu(k,311) - lu(k,93) * lu(k,308)
         lu(k,314) = lu(k,314) - lu(k,94) * lu(k,308)
         lu(k,318) = lu(k,318) - lu(k,95) * lu(k,308)
         lu(k,326) = lu(k,326) - lu(k,93) * lu(k,323)
         lu(k,332) = lu(k,332) - lu(k,94) * lu(k,323)
         lu(k,336) = lu(k,336) - lu(k,95) * lu(k,323)
         lu(k,377) = lu(k,377) - lu(k,93) * lu(k,359)
         lu(k,384) = lu(k,384) - lu(k,94) * lu(k,359)
         lu(k,388) = lu(k,388) - lu(k,95) * lu(k,359)
         lu(k,433) = lu(k,433) - lu(k,93) * lu(k,425)
         lu(k,440) = lu(k,440) - lu(k,94) * lu(k,425)
         lu(k,444) = lu(k,444) - lu(k,95) * lu(k,425)
                                                                        
         lu(k,96) = 1. / lu(k,96)
         lu(k,97) = lu(k,97) * lu(k,96)
         lu(k,98) = lu(k,98) * lu(k,96)
         lu(k,99) = lu(k,99) * lu(k,96)
         lu(k,100) = lu(k,100) * lu(k,96)
         lu(k,101) = lu(k,101) * lu(k,96)
         lu(k,102) = lu(k,102) * lu(k,96)
         lu(k,140) = lu(k,140) - lu(k,97) * lu(k,139)
         lu(k,141) = lu(k,141) - lu(k,98) * lu(k,139)
         lu(k,142) = lu(k,142) - lu(k,99) * lu(k,139)
         lu(k,143) = - lu(k,100) * lu(k,139)
         lu(k,144) = - lu(k,101) * lu(k,139)
         lu(k,146) = lu(k,146) - lu(k,102) * lu(k,139)
         lu(k,366) = lu(k,366) - lu(k,97) * lu(k,360)
         lu(k,374) = lu(k,374) - lu(k,98) * lu(k,360)
         lu(k,377) = lu(k,377) - lu(k,99) * lu(k,360)
         lu(k,378) = lu(k,378) - lu(k,100) * lu(k,360)
         lu(k,384) = lu(k,384) - lu(k,101) * lu(k,360)
         lu(k,388) = lu(k,388) - lu(k,102) * lu(k,360)
         lu(k,462) = lu(k,462) - lu(k,97) * lu(k,458)
         lu(k,466) = lu(k,466) - lu(k,98) * lu(k,458)
         lu(k,468) = lu(k,468) - lu(k,99) * lu(k,458)
         lu(k,469) = - lu(k,100) * lu(k,458)
         lu(k,475) = lu(k,475) - lu(k,101) * lu(k,458)
         lu(k,479) = lu(k,479) - lu(k,102) * lu(k,458)
                                                                        
         lu(k,103) = 1. / lu(k,103)
         lu(k,104) = lu(k,104) * lu(k,103)
         lu(k,105) = lu(k,105) * lu(k,103)
         lu(k,106) = lu(k,106) * lu(k,103)
         lu(k,225) = lu(k,225) - lu(k,104) * lu(k,223)
         lu(k,226) = - lu(k,105) * lu(k,223)
         lu(k,229) = lu(k,229) - lu(k,106) * lu(k,223)
         lu(k,245) = lu(k,245) - lu(k,104) * lu(k,243)
         lu(k,246) = - lu(k,105) * lu(k,243)
         lu(k,248) = lu(k,248) - lu(k,106) * lu(k,243)
         lu(k,252) = lu(k,252) - lu(k,104) * lu(k,250)
         lu(k,253) = - lu(k,105) * lu(k,250)
         lu(k,255) = lu(k,255) - lu(k,106) * lu(k,250)
         lu(k,384) = lu(k,384) - lu(k,104) * lu(k,361)
         lu(k,385) = lu(k,385) - lu(k,105) * lu(k,361)
         lu(k,389) = lu(k,389) - lu(k,106) * lu(k,361)
         lu(k,400) = lu(k,400) - lu(k,104) * lu(k,396)
         lu(k,401) = lu(k,401) - lu(k,105) * lu(k,396)
         lu(k,405) = lu(k,405) - lu(k,106) * lu(k,396)
         lu(k,503) = lu(k,503) - lu(k,104) * lu(k,487)
         lu(k,504) = lu(k,504) - lu(k,105) * lu(k,487)
         lu(k,508) = lu(k,508) - lu(k,106) * lu(k,487)
                                                                        
         lu(k,107) = 1. / lu(k,107)
         lu(k,108) = lu(k,108) * lu(k,107)
         lu(k,109) = lu(k,109) * lu(k,107)
         lu(k,110) = lu(k,110) * lu(k,107)
         lu(k,111) = lu(k,111) * lu(k,107)
         lu(k,112) = lu(k,112) * lu(k,107)
         lu(k,113) = lu(k,113) * lu(k,107)
         lu(k,114) = lu(k,114) * lu(k,107)
         lu(k,377) = lu(k,377) - lu(k,108) * lu(k,362)
         lu(k,378) = lu(k,378) - lu(k,109) * lu(k,362)
         lu(k,384) = lu(k,384) - lu(k,110) * lu(k,362)
         lu(k,385) = lu(k,385) - lu(k,111) * lu(k,362)
         lu(k,386) = lu(k,386) - lu(k,112) * lu(k,362)
         lu(k,388) = lu(k,388) - lu(k,113) * lu(k,362)
         lu(k,389) = lu(k,389) - lu(k,114) * lu(k,362)
         lu(k,398) = lu(k,398) - lu(k,108) * lu(k,397)
         lu(k,399) = - lu(k,109) * lu(k,397)
         lu(k,400) = lu(k,400) - lu(k,110) * lu(k,397)
         lu(k,401) = lu(k,401) - lu(k,111) * lu(k,397)
         lu(k,402) = lu(k,402) - lu(k,112) * lu(k,397)
         lu(k,404) = lu(k,404) - lu(k,113) * lu(k,397)
         lu(k,405) = lu(k,405) - lu(k,114) * lu(k,397)
         lu(k,411) = lu(k,411) - lu(k,108) * lu(k,410)
         lu(k,412) = - lu(k,109) * lu(k,410)
         lu(k,413) = lu(k,413) - lu(k,110) * lu(k,410)
         lu(k,414) = lu(k,414) - lu(k,111) * lu(k,410)
         lu(k,415) = lu(k,415) - lu(k,112) * lu(k,410)
         lu(k,417) = lu(k,417) - lu(k,113) * lu(k,410)
         lu(k,418) = lu(k,418) - lu(k,114) * lu(k,410)
                                                                        
         lu(k,115) = 1. / lu(k,115)
         lu(k,116) = lu(k,116) * lu(k,115)
         lu(k,117) = lu(k,117) * lu(k,115)
         lu(k,118) = lu(k,118) * lu(k,115)
         lu(k,119) = lu(k,119) * lu(k,115)
         lu(k,120) = lu(k,120) * lu(k,115)
         lu(k,121) = lu(k,121) * lu(k,115)
         lu(k,122) = lu(k,122) * lu(k,115)
         lu(k,326) = lu(k,326) - lu(k,116) * lu(k,324)
         lu(k,327) = lu(k,327) - lu(k,117) * lu(k,324)
         lu(k,328) = lu(k,328) - lu(k,118) * lu(k,324)
         lu(k,330) = lu(k,330) - lu(k,119) * lu(k,324)
         lu(k,331) = lu(k,331) - lu(k,120) * lu(k,324)
         lu(k,332) = lu(k,332) - lu(k,121) * lu(k,324)
         lu(k,336) = lu(k,336) - lu(k,122) * lu(k,324)
         lu(k,377) = lu(k,377) - lu(k,116) * lu(k,363)
         lu(k,379) = lu(k,379) - lu(k,117) * lu(k,363)
         lu(k,380) = lu(k,380) - lu(k,118) * lu(k,363)
         lu(k,382) = lu(k,382) - lu(k,119) * lu(k,363)
         lu(k,383) = lu(k,383) - lu(k,120) * lu(k,363)
         lu(k,384) = lu(k,384) - lu(k,121) * lu(k,363)
         lu(k,388) = lu(k,388) - lu(k,122) * lu(k,363)
         lu(k,468) = lu(k,468) - lu(k,116) * lu(k,459)
         lu(k,470) = lu(k,470) - lu(k,117) * lu(k,459)
         lu(k,471) = lu(k,471) - lu(k,118) * lu(k,459)
         lu(k,473) = lu(k,473) - lu(k,119) * lu(k,459)
         lu(k,474) = lu(k,474) - lu(k,120) * lu(k,459)
         lu(k,475) = lu(k,475) - lu(k,121) * lu(k,459)
         lu(k,479) = lu(k,479) - lu(k,122) * lu(k,459)
                                                                        
         lu(k,123) = 1. / lu(k,123)
         lu(k,124) = lu(k,124) * lu(k,123)
         lu(k,125) = lu(k,125) * lu(k,123)
         lu(k,126) = lu(k,126) * lu(k,123)
         lu(k,127) = lu(k,127) * lu(k,123)
         lu(k,168) = - lu(k,124) * lu(k,166)
         lu(k,171) = lu(k,171) - lu(k,125) * lu(k,166)
         lu(k,173) = lu(k,173) - lu(k,126) * lu(k,166)
         lu(k,175) = - lu(k,127) * lu(k,166)
         lu(k,186) = - lu(k,124) * lu(k,185)
         lu(k,193) = lu(k,193) - lu(k,125) * lu(k,185)
         lu(k,195) = - lu(k,126) * lu(k,185)
         lu(k,198) = - lu(k,127) * lu(k,185)
         lu(k,370) = lu(k,370) - lu(k,124) * lu(k,364)
         lu(k,384) = lu(k,384) - lu(k,125) * lu(k,364)
         lu(k,387) = lu(k,387) - lu(k,126) * lu(k,364)
         lu(k,390) = lu(k,390) - lu(k,127) * lu(k,364)
         lu(k,430) = - lu(k,124) * lu(k,426)
         lu(k,440) = lu(k,440) - lu(k,125) * lu(k,426)
         lu(k,443) = lu(k,443) - lu(k,126) * lu(k,426)
         lu(k,446) = lu(k,446) - lu(k,127) * lu(k,426)
         lu(k,464) = lu(k,464) - lu(k,124) * lu(k,460)
         lu(k,475) = lu(k,475) - lu(k,125) * lu(k,460)
         lu(k,478) = lu(k,478) - lu(k,126) * lu(k,460)
         lu(k,481) = lu(k,481) - lu(k,127) * lu(k,460)
         lu(k,567) = lu(k,567) - lu(k,124) * lu(k,562)
         lu(k,579) = lu(k,579) - lu(k,125) * lu(k,562)
         lu(k,582) = lu(k,582) - lu(k,126) * lu(k,562)
         lu(k,585) = lu(k,585) - lu(k,127) * lu(k,562)
                                                                        
      end do
                                                                        
      end subroutine IMP_LU_FAC03
                                                                        
      subroutine IMP_LU_FAC04( lu )
                                                                        
      use CHEM_MODS, only : imp_nzcnt, clsze
                                                                        
      implicit none
                                                                        
!-----------------------------------------------------------------------
!       ... Dummy args
!-----------------------------------------------------------------------
      real, intent(inout) ::   lu(clsze,imp_nzcnt)
                                                                        
!-----------------------------------------------------------------------
!       ... Local variables
!-----------------------------------------------------------------------
      integer :: k
                                                                        
      do k = 1,clsze
         lu(k,131) = 1. / lu(k,131)
         lu(k,132) = lu(k,132) * lu(k,131)
         lu(k,133) = lu(k,133) * lu(k,131)
         lu(k,134) = lu(k,134) * lu(k,131)
         lu(k,135) = lu(k,135) * lu(k,131)
         lu(k,136) = lu(k,136) * lu(k,131)
         lu(k,137) = lu(k,137) * lu(k,131)
         lu(k,138) = lu(k,138) * lu(k,131)
         lu(k,374) = lu(k,374) - lu(k,132) * lu(k,365)
         lu(k,377) = lu(k,377) - lu(k,133) * lu(k,365)
         lu(k,384) = lu(k,384) - lu(k,134) * lu(k,365)
         lu(k,385) = lu(k,385) - lu(k,135) * lu(k,365)
         lu(k,387) = lu(k,387) - lu(k,136) * lu(k,365)
         lu(k,388) = lu(k,388) - lu(k,137) * lu(k,365)
         lu(k,392) = lu(k,392) - lu(k,138) * lu(k,365)
         lu(k,431) = lu(k,431) - lu(k,132) * lu(k,427)
         lu(k,433) = lu(k,433) - lu(k,133) * lu(k,427)
         lu(k,440) = lu(k,440) - lu(k,134) * lu(k,427)
         lu(k,441) = lu(k,441) - lu(k,135) * lu(k,427)
         lu(k,443) = lu(k,443) - lu(k,136) * lu(k,427)
         lu(k,444) = lu(k,444) - lu(k,137) * lu(k,427)
         lu(k,448) = lu(k,448) - lu(k,138) * lu(k,427)
         lu(k,466) = lu(k,466) - lu(k,132) * lu(k,461)
         lu(k,468) = lu(k,468) - lu(k,133) * lu(k,461)
         lu(k,475) = lu(k,475) - lu(k,134) * lu(k,461)
         lu(k,476) = lu(k,476) - lu(k,135) * lu(k,461)
         lu(k,478) = lu(k,478) - lu(k,136) * lu(k,461)
         lu(k,479) = lu(k,479) - lu(k,137) * lu(k,461)
         lu(k,483) = lu(k,483) - lu(k,138) * lu(k,461)
         lu(k,569) = lu(k,569) - lu(k,132) * lu(k,563)
         lu(k,572) = lu(k,572) - lu(k,133) * lu(k,563)
         lu(k,579) = lu(k,579) - lu(k,134) * lu(k,563)
         lu(k,580) = lu(k,580) - lu(k,135) * lu(k,563)
         lu(k,582) = lu(k,582) - lu(k,136) * lu(k,563)
         lu(k,583) = lu(k,583) - lu(k,137) * lu(k,563)
         lu(k,587) = lu(k,587) - lu(k,138) * lu(k,563)
                                                                        
         lu(k,140) = 1. / lu(k,140)
         lu(k,141) = lu(k,141) * lu(k,140)
         lu(k,142) = lu(k,142) * lu(k,140)
         lu(k,143) = lu(k,143) * lu(k,140)
         lu(k,144) = lu(k,144) * lu(k,140)
         lu(k,145) = lu(k,145) * lu(k,140)
         lu(k,146) = lu(k,146) * lu(k,140)
         lu(k,147) = lu(k,147) * lu(k,140)
         lu(k,212) = lu(k,212) - lu(k,141) * lu(k,210)
         lu(k,214) = lu(k,214) - lu(k,142) * lu(k,210)
         lu(k,215) = - lu(k,143) * lu(k,210)
         lu(k,216) = lu(k,216) - lu(k,144) * lu(k,210)
         lu(k,217) = lu(k,217) - lu(k,145) * lu(k,210)
         lu(k,219) = lu(k,219) - lu(k,146) * lu(k,210)
         lu(k,222) = - lu(k,147) * lu(k,210)
         lu(k,374) = lu(k,374) - lu(k,141) * lu(k,366)
         lu(k,377) = lu(k,377) - lu(k,142) * lu(k,366)
         lu(k,378) = lu(k,378) - lu(k,143) * lu(k,366)
         lu(k,384) = lu(k,384) - lu(k,144) * lu(k,366)
         lu(k,385) = lu(k,385) - lu(k,145) * lu(k,366)
         lu(k,388) = lu(k,388) - lu(k,146) * lu(k,366)
         lu(k,392) = lu(k,392) - lu(k,147) * lu(k,366)
         lu(k,466) = lu(k,466) - lu(k,141) * lu(k,462)
         lu(k,468) = lu(k,468) - lu(k,142) * lu(k,462)
         lu(k,469) = lu(k,469) - lu(k,143) * lu(k,462)
         lu(k,475) = lu(k,475) - lu(k,144) * lu(k,462)
         lu(k,476) = lu(k,476) - lu(k,145) * lu(k,462)
         lu(k,479) = lu(k,479) - lu(k,146) * lu(k,462)
         lu(k,483) = lu(k,483) - lu(k,147) * lu(k,462)
         lu(k,569) = lu(k,569) - lu(k,141) * lu(k,564)
         lu(k,572) = lu(k,572) - lu(k,142) * lu(k,564)
         lu(k,573) = lu(k,573) - lu(k,143) * lu(k,564)
         lu(k,579) = lu(k,579) - lu(k,144) * lu(k,564)
         lu(k,580) = lu(k,580) - lu(k,145) * lu(k,564)
         lu(k,583) = lu(k,583) - lu(k,146) * lu(k,564)
         lu(k,587) = lu(k,587) - lu(k,147) * lu(k,564)
                                                                        
         lu(k,148) = 1. / lu(k,148)
         lu(k,149) = lu(k,149) * lu(k,148)
         lu(k,150) = lu(k,150) * lu(k,148)
         lu(k,151) = lu(k,151) * lu(k,148)
         lu(k,152) = lu(k,152) * lu(k,148)
         lu(k,295) = lu(k,295) - lu(k,149) * lu(k,292)
         lu(k,299) = lu(k,299) - lu(k,150) * lu(k,292)
         lu(k,301) = lu(k,301) - lu(k,151) * lu(k,292)
         lu(k,303) = lu(k,303) - lu(k,152) * lu(k,292)
         lu(k,311) = lu(k,311) - lu(k,149) * lu(k,309)
         lu(k,314) = lu(k,314) - lu(k,150) * lu(k,309)
         lu(k,316) = - lu(k,151) * lu(k,309)
         lu(k,318) = lu(k,318) - lu(k,152) * lu(k,309)
         lu(k,377) = lu(k,377) - lu(k,149) * lu(k,367)
         lu(k,384) = lu(k,384) - lu(k,150) * lu(k,367)
         lu(k,386) = lu(k,386) - lu(k,151) * lu(k,367)
         lu(k,388) = lu(k,388) - lu(k,152) * lu(k,367)
         lu(k,433) = lu(k,433) - lu(k,149) * lu(k,428)
         lu(k,440) = lu(k,440) - lu(k,150) * lu(k,428)
         lu(k,442) = lu(k,442) - lu(k,151) * lu(k,428)
         lu(k,444) = lu(k,444) - lu(k,152) * lu(k,428)
         lu(k,496) = lu(k,496) - lu(k,149) * lu(k,488)
         lu(k,503) = lu(k,503) - lu(k,150) * lu(k,488)
         lu(k,505) = lu(k,505) - lu(k,151) * lu(k,488)
         lu(k,507) = lu(k,507) - lu(k,152) * lu(k,488)
         lu(k,517) = lu(k,517) - lu(k,149) * lu(k,515)
         lu(k,524) = lu(k,524) - lu(k,150) * lu(k,515)
         lu(k,526) = lu(k,526) - lu(k,151) * lu(k,515)
         lu(k,528) = lu(k,528) - lu(k,152) * lu(k,515)
         lu(k,542) = lu(k,542) - lu(k,149) * lu(k,535)
         lu(k,549) = lu(k,549) - lu(k,150) * lu(k,535)
         lu(k,551) = lu(k,551) - lu(k,151) * lu(k,535)
         lu(k,553) = lu(k,553) - lu(k,152) * lu(k,535)
         lu(k,572) = lu(k,572) - lu(k,149) * lu(k,565)
         lu(k,579) = lu(k,579) - lu(k,150) * lu(k,565)
         lu(k,581) = lu(k,581) - lu(k,151) * lu(k,565)
         lu(k,583) = lu(k,583) - lu(k,152) * lu(k,565)
                                                                        
         lu(k,153) = 1. / lu(k,153)
         lu(k,154) = lu(k,154) * lu(k,153)
         lu(k,155) = lu(k,155) * lu(k,153)
         lu(k,156) = lu(k,156) * lu(k,153)
         lu(k,157) = lu(k,157) * lu(k,153)
         lu(k,158) = lu(k,158) * lu(k,153)
         lu(k,159) = lu(k,159) * lu(k,153)
         lu(k,160) = lu(k,160) * lu(k,153)
         lu(k,161) = lu(k,161) * lu(k,153)
         lu(k,162) = lu(k,162) * lu(k,153)
         lu(k,163) = lu(k,163) * lu(k,153)
         lu(k,164) = lu(k,164) * lu(k,153)
         lu(k,373) = lu(k,373) - lu(k,154) * lu(k,368)
         lu(k,375) = - lu(k,155) * lu(k,368)
         lu(k,377) = lu(k,377) - lu(k,156) * lu(k,368)
         lu(k,379) = lu(k,379) - lu(k,157) * lu(k,368)
         lu(k,380) = lu(k,380) - lu(k,158) * lu(k,368)
         lu(k,383) = lu(k,383) - lu(k,159) * lu(k,368)
         lu(k,384) = lu(k,384) - lu(k,160) * lu(k,368)
         lu(k,386) = lu(k,386) - lu(k,161) * lu(k,368)
         lu(k,388) = lu(k,388) - lu(k,162) * lu(k,368)
         lu(k,389) = lu(k,389) - lu(k,163) * lu(k,368)
         lu(k,391) = lu(k,391) - lu(k,164) * lu(k,368)
         lu(k,492) = lu(k,492) - lu(k,154) * lu(k,489)
         lu(k,494) = lu(k,494) - lu(k,155) * lu(k,489)
         lu(k,496) = lu(k,496) - lu(k,156) * lu(k,489)
         lu(k,498) = lu(k,498) - lu(k,157) * lu(k,489)
         lu(k,499) = lu(k,499) - lu(k,158) * lu(k,489)
         lu(k,502) = lu(k,502) - lu(k,159) * lu(k,489)
         lu(k,503) = lu(k,503) - lu(k,160) * lu(k,489)
         lu(k,505) = lu(k,505) - lu(k,161) * lu(k,489)
         lu(k,507) = lu(k,507) - lu(k,162) * lu(k,489)
         lu(k,508) = lu(k,508) - lu(k,163) * lu(k,489)
         lu(k,510) = lu(k,510) - lu(k,164) * lu(k,489)
         lu(k,538) = lu(k,538) - lu(k,154) * lu(k,536)
         lu(k,540) = - lu(k,155) * lu(k,536)
         lu(k,542) = lu(k,542) - lu(k,156) * lu(k,536)
         lu(k,544) = lu(k,544) - lu(k,157) * lu(k,536)
         lu(k,545) = lu(k,545) - lu(k,158) * lu(k,536)
         lu(k,548) = - lu(k,159) * lu(k,536)
         lu(k,549) = lu(k,549) - lu(k,160) * lu(k,536)
         lu(k,551) = lu(k,551) - lu(k,161) * lu(k,536)
         lu(k,553) = lu(k,553) - lu(k,162) * lu(k,536)
         lu(k,554) = lu(k,554) - lu(k,163) * lu(k,536)
         lu(k,556) = lu(k,556) - lu(k,164) * lu(k,536)
                                                                        
         lu(k,167) = 1. / lu(k,167)
         lu(k,168) = lu(k,168) * lu(k,167)
         lu(k,169) = lu(k,169) * lu(k,167)
         lu(k,170) = lu(k,170) * lu(k,167)
         lu(k,171) = lu(k,171) * lu(k,167)
         lu(k,172) = lu(k,172) * lu(k,167)
         lu(k,173) = lu(k,173) * lu(k,167)
         lu(k,174) = lu(k,174) * lu(k,167)
         lu(k,175) = lu(k,175) * lu(k,167)
         lu(k,176) = lu(k,176) * lu(k,167)
         lu(k,370) = lu(k,370) - lu(k,168) * lu(k,369)
         lu(k,374) = lu(k,374) - lu(k,169) * lu(k,369)
         lu(k,377) = lu(k,377) - lu(k,170) * lu(k,369)
         lu(k,384) = lu(k,384) - lu(k,171) * lu(k,369)
         lu(k,385) = lu(k,385) - lu(k,172) * lu(k,369)
         lu(k,387) = lu(k,387) - lu(k,173) * lu(k,369)
         lu(k,388) = lu(k,388) - lu(k,174) * lu(k,369)
         lu(k,390) = lu(k,390) - lu(k,175) * lu(k,369)
         lu(k,392) = lu(k,392) - lu(k,176) * lu(k,369)
         lu(k,430) = lu(k,430) - lu(k,168) * lu(k,429)
         lu(k,431) = lu(k,431) - lu(k,169) * lu(k,429)
         lu(k,433) = lu(k,433) - lu(k,170) * lu(k,429)
         lu(k,440) = lu(k,440) - lu(k,171) * lu(k,429)
         lu(k,441) = lu(k,441) - lu(k,172) * lu(k,429)
         lu(k,443) = lu(k,443) - lu(k,173) * lu(k,429)
         lu(k,444) = lu(k,444) - lu(k,174) * lu(k,429)