/*codriv - public routine, calls comole to converge molecules */ /*comole fills in matrix for heavy elements molecular routines */ #include #include "cddefines.h" #define NMOLE 15 #include "showme.h" #include "physconst.h" #include "cmfrac.h" #include "secondaries.h" #include "h.h" #include "tepowers.h" #include "hmi.h" #include "heating.h" #include "trace.h" #include "nomole.h" #include "he.h" #include "phycon.h" #include "coatom.h" #include "hevmolec.h" #include "moldrv.h" #include "called.h" #include "ghabng.h" #include "h2ozer.h" #include "codish.h" #include "covib.h" #include "timesc.h" #include "zonecnt.h" #include "ionfracs.h" #include "destcrt.h" #include "typmatrx.h" #include "linpack.h" #include "sexp.h" #include "matin1.h" #include "veclib.h" #include "hmrate.h" #include "lotsco.h" #include "comneg.h" #include "colimt.h" #include "comole.h" /*molav average old and new molecular equilibrium balance from comole */ void molav(long int JobMolAv); void comole(int *lgNegPop, int *lgZerPop) { char chLab[NMOLE][5]; long int i, icatm, ich, ich2p, ichp, ico, icop, icp, ih2o, ih2op, ih3op, io2, io2p, ioatm, ioh, iohp, ipiv[NMOLE], j, job, loop, merror; double cartot, cch, cf, ch2pc, ch2pch, ch2pco, ch2ph3, chch2, chchp, chhc, chhchp, choco, chohwa, chooh, chpe, chph, chphch, chplte, co2co, coch, cocn, cocop, cocp, coo2, cooh, copco, coplte, coprm, cpch2p, cpchp, crden, csav[NMOLE][NMOLE], e3opoh, e3opwa, eo2o, expc2, expchp, expco, expcop, expo2, h2lim, h2pchp, h3oh3p, h3opwa, h3owai, h3pchp, h3pcop, h3pohp, h3pwai, h3pwap, ho2p, o2co, o2cop, o2o2p, o2oh, o2rm, oh2wat, ohco, ohcop, oho, oho2, ohohp, ohohwa, ohph3p, ohpoh, ohprm, ohpwap, oohp, opohp, oxytot, pc2lte, pcolte, phodco, po2lte, ratio, wa2oh, wachoh, wacop, wah3op, waiwai, waoh, waohp, waph3p, wapoh, wapohp, waprm, wapwa, wawap; double amat[NMOLE][NMOLE], bvec[NMOLE], c[NMOLE + 1][NMOLE + 1], rcond, work[NMOLE]; float oldcarb; # ifdef DEBUG_FUN fputs( "<+>comole()\n", debug_fp ); # endif /* take averages of molecular abundances, this routine follows below */ molav(1); loop = 0; while( moldrv.lgDriveMole || (loop < 1) ) { loop += 1; /* this block intentionally killed */ /*if( moldrv.lgDriveMole ) { // option to allow user to enter oh abundances fprintf( ioQQQ, " enter log of OH abundance\n" ); readf( 5, "(a80)", "%S ", chCard,81 ); i = 1; // option to set OH here fac = FFmtRead(chCard,&i,76,&lgEOL); if( fac == 1. ) { # ifdef DEBUG_FUN fputs( " <->comole()\n", debug_fp ); # endif return; } else if( fac == 0. ) { molav(2); } else { hevmolec.hevmol[IPOH-1] = (float)pow(10.,fac); } fprintf( ioQQQ, " enter log of H2O abundance\n" ); readf( 5, "(a80)", "%S ", chCard,81 ); i = 1; // option to set H2O here fac = FFmtRead(chCard,&i,76,&lgEOL); if( fac != 0. ) hevmolec.hevmol[IPH2O-1] = (float)pow(10.,fac); }*/ hevmolec.hevmol[IPCTWO-1] = 0.; hevmolec.hevmol[IPC2P-1] = 0.; hevmolec.hevmol[IPOTWO-1] = 0.; hevmolec.hevmol[IPO2P-1] = 0.; hevmolec.hevmol[IPH3P-1] = 0.; /* H2LIM is smallest ratio of H2/HDEN for which we will * even try to invert heavy element network */ if( nomole.lgNoH2Mole ) { /* H2 network is turned off, ignore this limit */ h2lim = 0.; } else { h2lim = 1e-4; } /* do we want to try to calculate the CO? */ if( nomole.lgNoCOMole || (phycon.te > 3800.) || (hmi.htwo/phycon.hden < h2lim) || /* these are the atomic carbon and oxygen abundances */ (xIonFracs[1][5]*xIonFracs[1][7] == 0.) ) { coatom.catmic = xIonFracs[1][5]; coatom.oatmic = xIonFracs[1][7]; hevmolec.hevmol[IPCH-1] = 0.; hevmolec.hevmol[IPCHP-1] = 0.; hevmolec.hevmol[IPOH-1] = 0.; hevmolec.hevmol[IPOHP-1] = 0.; hevmolec.hevmol[IPOTWO-1] = 0.; hevmolec.hevmol[IPCO-1] = 0.; hevmolec.hevmol[IPCOP-1] = 0.; hevmolec.hevmol[IPH2O-1] = 0.; hevmolec.hevmol[IPH2OP-1] = 0.; hevmolec.hevmol[IPO2P-1] = 0.; hevmolec.hevmol[IPC2P-1] = 0.; hevmolec.hevmol[IPH3OP-1] = 0.; hevmolec.hevmol[IPCH2P-1] = 0.; hevmolec.hevmol[IPCH2-1] = 0.; hevmolec.hevmol[IPCH3-1] = 0.; hevmolec.hevmol[IPCP-1] = 0.; /* heating due to CO photodissociation */ HeatingCom.heating[9][0] = 0.; covib.COVibCool = 0.; covib.dCOVdT = 0.; codish.CODissHeat = 0.; # ifdef DEBUG_FUN fputs( " <->comole()\n", debug_fp ); # endif return; } /* dissociation energies; * CO 11.09, CO+ 8.3eV, CH+ 3.8eV * LTE population of CO+ */ /* following, pc2lte, chplte, po2lte, coplte, pcolte, are not ever * used for anything - first document what they are relative to, and * then use them to define departure coefficients */ /* turn this error off for now */ /*lint -e550 */ expcop = sexp(9.632e4/phycon.te); if( expcop > 0. ) { coplte = SAHA/(TePowers.te32*expcop)*(1./(6.*9.)); } else { coplte = 0.; } expco = sexp(1.287e5/phycon.te); if( expco > 0. ) { pcolte = SAHA/(TePowers.te32*expco)*(1./(9.*9.)); } else { pcolte = 0.; } expchp = sexp(4.4098e4/phycon.te); if( expchp > 0. ) { chplte = SAHA/(TePowers.te32*expchp)*(1./(1.*6.)); } else { chplte = 0.; } expc2 = sexp(4.4098e4/phycon.te); if( expc2 > 0. ) { pc2lte = SAHA/(TePowers.te32*expc2)*(1./(9.*9.)); } else { pc2lte = 0.; } expo2 = sexp(4.4098e4/phycon.te); if( expo2 > 0. ) { po2lte = SAHA/(TePowers.te32*expo2)*(3./(9.*9.)); } else { po2lte = 0.; } /* molecule formation */ for( i=0; i < NMOLE; i++ ) { c[NMOLE][i] = 0.; for( j=0; j < NMOLE; j++ ) { c[j][i] = 0.; } } /* set up pointers */ icatm = 1; strcpy( chLab[icatm-1], "C " ); ioatm = 2; strcpy( chLab[ioatm-1], "O " ); ich = 3; strcpy( chLab[ich-1], "CH " ); ichp = 4; strcpy( chLab[ichp-1], "CH+ " ); ioh = 5; strcpy( chLab[ioh-1], "OH " ); iohp = 6; strcpy( chLab[iohp-1], "OH+ " ); ich2p = 7; strcpy( chLab[ich2p-1], "CH2+" ); icp = 8; strcpy( chLab[icp-1], "C+ " ); ico = 9; strcpy( chLab[ico-1], "CO " ); icop = 10; strcpy( chLab[icop-1], "CO+ " ); ih2o = 11; strcpy( chLab[ih2o-1], "H2O " ); ih2op = 12; strcpy( chLab[ih2op-1], "H2O+" ); ih3op = 13; strcpy( chLab[ih3op-1], "H3O+" ); io2 = 14; strcpy( chLab[io2-1], "O2 " ); io2p = 15; strcpy( chLab[io2p-1], "O2+ " ); /* carbon conservation */ cartot = 0.; for( i=3; i <= 7; i++ ) { cartot += xIonFracs[i][5]; } cartot = xIonFracs[0][5] - cartot; c[NMOLE][icatm-1] = cartot; /* oxygen conservation */ if( coatom.oatmic == 0. ) coatom.oatmic = xIonFracs[1][7]; oxytot = 0.; for( i=2; i <= 9; i++ ) { oxytot += xIonFracs[i][7]; } oxytot = xIonFracs[0][7] - oxytot; c[NMOLE][ioatm-1] = oxytot; /*-------------------------------------------------------------------- * * first equation is carbon conservation */ c[icatm-1][icatm-1] = 1.; c[ioatm-1][icatm-1] = 0.; c[ich-1][icatm-1] = 1.; c[ichp-1][icatm-1] = 1.; c[ioh-1][icatm-1] = 0.; c[iohp-1][icatm-1] = 0.; c[io2-1][icatm-1] = 0.; c[io2p-1][icatm-1] = 0.; c[ico-1][icatm-1] = 1.; c[icop-1][icatm-1] = 1.; c[ih2o-1][icatm-1] = 0.; c[ih2op-1][icatm-1] = 0.; c[ih3op-1][icatm-1] = 0.; c[ich2p-1][icatm-1] = 1.; c[icp-1][icatm-1] = 1.; /*-------------------------------------------------------------------- * * second equation is oxygen conservation */ c[icatm-1][ioatm-1] = 0.; c[ioatm-1][ioatm-1] = 1.; c[ich-1][ioatm-1] = 0.; c[ichp-1][ioatm-1] = 0.; c[ioh-1][ioatm-1] = 1.; c[iohp-1][ioatm-1] = 1.; c[io2-1][ioatm-1] = 2.; c[io2p-1][ioatm-1] = 2.; c[ico-1][ioatm-1] = 1.; c[icop-1][ioatm-1] = 1.; c[ih2o-1][ioatm-1] = 1.; c[ih2op-1][ioatm-1] = 1.; c[ih3op-1][ioatm-1] = 1.; c[ich2p-1][ioatm-1] = 0.; c[icp-1][ioatm-1] = 0.; /*-------------------------------------------------------------------- * * balance equations; * (1,n) is C, called ICATM * (2,n) is O, called IOATM * (3,n) is CH, called ICH * (4,n) is CHP, called ICHP * (5,n) is OH, called IOH * (6,n) is OHP, called IOHP * (7,n) is CH2P, called ICH2P * (8,n) is CP, called ICP * (9,n) is CO, called ICO * (10,n) is COP, called ICOP * (11,n) is H2O, called IH2O * (12,n) is H2OP, called IH2OP * (13,n) is H3OP, called IH3OP * (14,n) is OTWO, IO2 * (15,n) is O2P, called IO2P * *-------------------------------------------------------------------- * * CH balance equations * * sum of CH + HII => ch+ and both H+ and Ho * also electrons to 2e */ chchp = hmrate(8.83e-14,0.5,1.29e5)*(h.hii + h.hi + hmi.htwo) + 1.90e-9*h.hii + hmrate(5.2e-10,0.5,1.50e5)*phycon.eden; c[ich-1][ich-1] -= chchp; /* o + ch => co + h */ choco = hmrate(9.53e-11,0.5,0.)*coatom.oatmic; c[ich-1][ich-1] -= choco; /* O + CH => C + OH */ chooh = hmrate(1.73e-11,0.50,4000.)*coatom.oatmic; c[ich-1][ich-1] -= chooh; /* CH + H => C + H2 */ chhc = hmrate(1.80e-11,0.50,4000.)*h.hi; /* there is no carbon balance equation */ c[ich-1][ich-1] -= chhc; /* CH + H => CH+ + H + e */ chhchp = hmrate(8.83e-14,0.5,1.29e5)*h.hi; c[ich-1][ich-1] -= chhchp; /* C + H2 => CH + H */ cch = hmrate(4.50e-11,0.5,1.56e4)*hmi.htwo; c[icatm-1][ich-1] += cch; /* CO + H => CH + O */ coch = hmrate(9.53e-11,0.5,8.83e4)*h.hi; c[ico-1][ich-1] += coch; /* H2O + C => CH + OH */ wachoh = hmrate(1.43e-10,0.50,2.40e4)*coatom.catmic; c[ih2o-1][ich-1] += wachoh; /* e- + CH2+ => CH + H */ ch2pch = hmrate(2.50e-7,-0.5,0.)*phycon.eden; c[ich2p-1][ich-1] += ch2pch; /* H2 + CH2+ => CH + H3+ */ ch2ph3 = hmrate(1.20e-9,0.,4.22e4)*hmi.htwo; c[ich2p-1][ich-1] += ch2ph3; /* CH + H2 => CH2 + H */ chch2 = hmrate(3.60e-10,0.,3.90e3)*hmi.htwo; c[ich-1][ich-1] -= chch2; /* OH + CH => H2O + C */ chohwa = 1e-10*hevmolec.hevmol[IPOH-1]; c[ich-1][ich-1] -= chohwa; /*-------------------------------------------------------------------- * * CH+ balance equation * * CH+ + e => C + H */ chpe = hmrate(3.00e-7,-0.40,0.)*phycon.eden; c[ichp-1][ichp-1] -= chpe; /* CH+ + H => H2 + C+ */ chph = hmrate(6.00e-10,-0.25,0.)*h.hi; c[ichp-1][ichp-1] -= chph; /* make CH+ from CH via collisions with protons and electrons, H */ c[ich-1][ichp-1] += chchp + chhchp; /* H2 + C+ => H + CH+ */ cpchp = hmrate(2.00e-10,0.,4640.)*hmi.htwo; c[icp-1][ichp-1] += cpchp; /* C + H3+ => H2 + CH+ */ h3pchp = 2.00e-9*hmi.h3plus; c[icatm-1][ichp-1] += h3pchp; /* H + CH2+ => H2 + CH+ */ h2pchp = hmrate(1.00e-9,0.,1.20e4)*h.hi; c[ich2p-1][ichp-1] += h2pchp; /* H2 + CH+ => H + CH2+ */ chphch = 1.00e-9*hmi.htwo; c[ichp-1][ichp-1] -= chphch; /*-------------------------------------------------------------------- * * OH balance equations * * various processes which change OH to OH+ */ ohohp = 2.10e-10*h.hii + hmrate(8.83e-14,0.5,1.53e5)*(h.hii + h.hi + hmi.htwo) + 7.60e-10*hmi.h2plus + hmrate(5.20e-10, 0.5,1.5e5)*phycon.eden; c[ioh-1][ioh-1] -= ohohp; /* this is first major step to making CO * OH + C+ => H + CO+ */ ohcop = 7.70e-10*hevmolec.hevmol[IPCP-1]; c[ioh-1][ioh-1] -= ohcop; /* OH + CH => H2O + C */ chohwa = 1e-10*hevmolec.hevmol[IPCH-1]; c[ioh-1][ioh-1] -= chohwa; /* C + OH => CO + H */ ohco = hmrate(1.11e-10,0.5,0.)*coatom.catmic; c[ioh-1][ioh-1] -= ohco; /* O + OH => O2 + H */ oho2 = hmrate(4.33e-11,-0.5,30.0)*coatom.oatmic; c[ioh-1][ioh-1] -= oho2; /* OH + H => O + H2 */ oho = hmrate(6.60e-13,1.53,2970.)*h.hi; /* this is really OH + cr => O + H, enhance from Gredel et al. ApJ 347, 289 */ oho += Secondaries.csupra*500.; c[ioh-1][ioh-1] -= oho; /* two ways to make water, H2 and OH * first is OH + H2 => H2O + H */ oh2wat = hmrate(8.80e-13,1.95,1420.)*hmi.htwo; /* second is OH + OH => H2O + O */ ohohwa = hmrate(4.20e-12,0.,242.)*hevmolec.hevmol[IPOH-1]; c[ioh-1][ioh-1] += -oh2wat - ohohwa*2.; /* O2 + H => OH + O */ o2oh = hmrate(1.63e-9,-0.9,8750.)*h.hi; c[io2-1][ioh-1] += o2oh; /* OH+ + stuff goes to OH */ ohpoh = hmrate(2.10e-9,0.,2800.)*h.hi + hmrate(7.60e-10,0., 2.36e4)*hmi.htwo; c[iohp-1][ioh-1] += ohpoh; /* O + CH => C + OH */ c[ich-1][ioh-1] += chooh; /* CO + H => C + OH */ cooh = hmrate(1.11e-10,0.5,7.77e4)*h.hi; c[ico-1][ioh-1] += cooh; /* water into OH */ waoh = 2.03e-10*he.heii + hmrate(9.20e-9,-0.5,6.96e4)*phycon.eden + hmrate(7.44e-12,1.57,9140.)*h.hi; c[ih2o-1][ioh-1] += waoh; /* H2O + O => OH + OH */ wa2oh = hmrate(1.35e-12,1.75,7860.)*coatom.oatmic; c[ih2o-1][ioh-1] += wa2oh*2.; /* H2O + C => CH + OH */ c[ih2o-1][ioh-1] += wachoh; /* ionized water into OH */ wapoh = hmrate(7.60e-10,0.,4e4)*h.hi + hmrate(1.30e-9,0.,2.03e4)* hmi.htwo + hmrate(2e-7,-0.5,0.)*phycon.eden; c[ih2op-1][ioh-1] += wapoh; /* H2 + H2O+ => OH + H3+ */ waph3p = hmrate(1.30e-9,0.,2.03e4)*hmi.htwo; c[ih2op-1][ioh-1] += waph3p; /* e- + H3O+ => OH + H3 */ e3opoh = hmrate(3.00e-7,-0.5,0.)*phycon.eden; c[ih3op-1][ioh-1] += e3opoh; /* OH + H3P => H2 + H2O+ */ h3pwap = 1.0e-9*hmi.h3plus; c[ioh-1][ioh-1] -= h3pwap; /*-------------------------------------------------------------------- * * OH+ balance equations * * following remove OH+ from network * both collide with H, go to H2 or H2+ */ ohprm = hmrate(1.60e-9,0.,2.26e4)*h.hi + hmrate(4.90e-10,-0.03, 1970.)*h.hi + hmrate(2.00e-7,-0.5,0.)*phycon.eden; c[iohp-1][iohp-1] -= ohprm; /* OH+ + H => OH + H+ and OH+ + H2 => OH + H2+ */ c[iohp-1][iohp-1] -= ohpoh; /* h2 + oh+ => h + h2o+ */ ohpwap = 1.10e-9*hmi.htwo; c[iohp-1][iohp-1] -= ohpwap; /* several processes converting OH into OH+ */ c[ioh-1][iohp-1] += ohohp; /* water into OH+ */ waohp = 2.86e-10*he.heii; c[ih2o-1][iohp-1] += waohp; /* H + H2O+ => H2 + OH+ */ wapohp = hmrate(1.70e-9,0.29,1.40e4)*h.hi; c[ih2op-1][iohp-1] += wapohp; /* O+ + H2 => H + OH+ */ if( xIonFracs[2][7] > 0. ) { ratio = xIonFracs[2][7]/xIonFracs[1][7]; opohp = 1.60e-9*ratio*hmi.htwo; c[ioatm-1][iohp-1] += opohp; } /* O + H2+ => H + OH+ */ oohp = 1.50e-9*hmi.h2plus; c[ioatm-1][iohp-1] += oohp; /* O + H3P => H2 + OHP */ h3pohp = 8.00e-10*hmi.h3plus; c[ioatm-1][iohp-1] += h3pohp; /* H2 + OH+ => O + H3P */ ohph3p = hmrate(8.00e-10,0.,2900.)*hmi.htwo; c[iohp-1][iohp-1] -= ohph3p; /*-------------------------------------------------------------------- * * O2 balance equation * * reactions that remove O2 from network */ o2o2p = 1.20e-9*h.hii + hmrate(2.37e-15,1.04,1.40e5)*h.hii + 2.70e-9*hmi.h2plus + hmrate(1.40e-11,1.04,1.40e5)*phycon.eden + hmrate(2.37e-15,1.04,1.40e5)*(h.hi + hmi.htwo); c[io2-1][io2-1] -= o2o2p; o2rm = 1.10e-9*he.heii + hmrate(2.30e-9,-0.5,6.96e4)*phycon.eden; /* cr destruction with enhancement from Gredel et al. */ o2rm += Secondaries.csupra*730.; c[io2-1][io2-1] -= o2rm; /* O2 + C+ => O + CO+ */ o2cop = 3.76e-10*hevmolec.hevmol[IPCP-1]; c[io2-1][io2-1] -= o2cop; /* O2 + C+ => CO + O+ */ o2co = 6.14e-10*hevmolec.hevmol[IPCP-1]; c[io2-1][io2-1] -= o2co; /* C + O2 => CO + O */ co2co = hmrate(1.80e-11,0.5,0.)*coatom.catmic; c[io2-1][io2-1] -= co2co; /* O2 + H => OH + O */ c[io2-1][io2-1] -= o2oh; /* o + oh => o2 + h */ c[ioh-1][io2-1] += oho2; /* o + co => c + o2 */ coo2 = hmrate(2.90e-11,0.5,6.93e4)*coatom.oatmic; c[ico-1][io2-1] += coo2; /* h + o2+ => o2 + h+ */ ho2p = hmrate(2.80e-10,-0.04,1.78e4)*h.hi; c[io2p-1][io2-1] += ho2p; /*-------------------------------------------------------------------- * * O2+ balance * * e + O2+ => O + O */ eo2o = hmrate(2.00e-7,-0.5,0.)*phycon.eden; c[io2p-1][io2p-1] -= eo2o; /* H + O2+ => O2 + H+ */ c[io2p-1][io2p-1] -= ho2p; /* various processes creating O2+ from O2 */ c[io2-1][io2p-1] += o2o2p; /*------------------------------------------------------------------ * * CO balance equation carbon monoxide * * CO + H+ => H + CO+, also CO + H+ => CO+ + */ cocop = hmrate(1.90e-10,0.,4660.)*h.hii + hmrate(1.13e-13, 0.6,1.63e5)*(h.hii + h.hi + hmi.htwo) + 2.80e-9*hmi.h2plus + hmrate(6.70e-10,0.6,1.63e5)*phycon.eden; c[ico-1][ico-1] -= cocop; /* co + hE+ => o + c+ + hE */ cocp = 1.70e-9*he.heii*1.1; c[ico-1][ico-1] -= cocp; /* CO + e => C + O + e */ cocn = hmrate(8.1e-10,-0.50,1.14e5)*phycon.eden; /* estra factor is enhancement due to light */ cocn += Secondaries.csupra*20.; c[ico-1][ico-1] -= cocn; /* CO + hnu => C + O, from Teilens and Hollenback, Table 5 */ phodco = ghabng.GammaHabing*1.4e-11; c[ico-1][ico-1] -= phodco; /* O + CO => C + O2 */ c[ico-1][ico-1] -= coo2; /* CO + H => C + OH */ c[ico-1][ico-1] -= cooh; /* CO + H => CH + O */ c[ico-1][ico-1] -= coch; /* O + CH => CO + H */ c[ich-1][ico-1] += choco; /* C + OH => CO + H */ c[ioh-1][ico-1] += ohco; /* O2 + C+ => CO + O+ */ o2co = 6.14e-10*hevmolec.hevmol[IPOTWO-1]; c[icp-1][ico-1] += o2co; /* C + O2 => CO + O */ c[io2-1][ico-1] += co2co; /* CO+ + H => CO + H+ */ copco = 1.90e-10*h.hi; c[icop-1][ico-1] += copco; /* CO + H3+ => H2 + CO+ */ h3pcop = 1.70e-9*hmi.h3plus; c[ico-1][ico-1] -= h3pcop; /* remember longest CO timescale */ if( -c[ico-1][ico-1] > SMALLFLOAT ) { /* this is rate CO is destroyed, equal to formation rate in equilibrium */ timesc.AgeCOMoleDest = -1./c[ico-1][ico-1]; /* moved to radinc */ /*timesc.BigCOMoleForm = (float)MAX2(timesc.BigCOMoleForm,timesc.AgeCOMoleForm);*/ } else { timesc.AgeCOMoleDest = 0.; } /*-------------------------------------------------------------------- * * CO+ balance equation * * e + CO+ => C + O */ coprm = hmrate(1.80e-7,-0.5,0.)*phycon.eden; c[icop-1][icop-1] -= coprm; /* CO+ + H => CO + H+ */ c[icop-1][icop-1] -= copco; /* co + h+ => h + co+, ALSO co + h+ => co++ */ c[ico-1][icop-1] += cocop; /* OH + C+ => H + CO+ */ ohcop = 7.70e-10*hevmolec.hevmol[IPOH-1]; c[icp-1][icop-1] += ohcop; /* H2O + C+ => H2 + CO+ */ wacop = 2.70e-9*hevmolec.hevmol[IPCP-1]; c[ih2o-1][icop-1] += wacop; /* CO + H3+ => H2 + CO+ */ c[ico-1][icop-1] += h3pcop; /* o + ch2+ => h + co+ */ ch2pco = 7.50e-10*coatom.oatmic; c[ich2p-1][icop-1] += ch2pco; /* O2 + C+ => O + CO+ */ o2cop = 3.76e-10*hevmolec.hevmol[IPCP-1]; c[io2-1][icop-1] += o2cop; /*-------------------------------------------------------------------- * * water H2O balance * * processes converting wat to wap */ wawap = 8.20e-9*h.hii + hmrate(7.55e-14,0.45,1.47e5)*(h.hi + hmi.htwo + h.hii) + 3.90e-9*hmi.h2plus + hmrate(4.40e-10, 0.45,1.47e5)*phycon.eden + 6.05e-11*he.heii; c[ih2o-1][ih2o-1] -= wawap; /* h2o => oh+ */ c[ih2o-1][ih2o-1] -= waohp; /* water into OH */ c[ih2o-1][ih2o-1] -= waoh; /* H2O + C => CH + OH */ c[ih2o-1][ih2o-1] -= wachoh; /* H2O + O => OH + OH */ c[ih2o-1][ih2o-1] -= wa2oh; /* h2o + c+ => h2 + co+ */ c[ih2o-1][ih2o-1] -= wacop; if( c[ih2o-1][ih2o-1] == 0. ) { /* possible for zero water destruction, this is an error */ c[ih2o-1][ih2o-1] = -1e-35; h2ozer.lgH2Ozer = TRUE; } /* end of water destruction mechanisms * CH + OH => H2O + C */ chohwa = 1e-10*hevmolec.hevmol[IPOH-1]; c[ich-1][ih2o-1] += chohwa; /* OH + stuff => water */ c[ioh-1][ih2o-1] += oh2wat + ohohwa; /* ionized water into water */ wapwa = hmrate(8.20e-9,0.,1.15e4)*h.hi + hmrate(3.90e-9,0., 3.27e4)*hmi.htwo; c[ih2op-1][ih2o-1] += wapwa; /* H2 + H3O => H2O + H3+ */ h3oh3p = hmrate(5.90e-9,0.,3.28e4)*hmi.htwo; c[ih3op-1][ih2o-1] += h3oh3p; /* H2O + H2+ => H + H3O+ */ wah3op = 3.40e-9*hmi.h2plus; c[ih2o-1][ih2o-1] -= wah3op; /* H + H3O+ => H2O + H2+ */ h3opwa = hmrate(5.90e-9,0.,3.28e4)*hmi.htwo; c[ih3op-1][ih2o-1] += h3opwa; /* e- + H3O+ => H2O + H */ e3opwa = hmrate(1.00e-6,-0.5,0.)*phycon.eden; c[ih3op-1][ih2o-1] += e3opwa; /* H2O + H3P => H2 + H3OP */ h3pwai = 5.90e-9*hmi.h3plus; c[ih2o-1][ih2o-1] -= h3pwai; /*-------------------------------------------------------------------- * H2O+ balance * * removes ionized water from network */ waprm = hmrate(2e-7,-0.5,0.)*phycon.eden; c[ih2op-1][ih2op-1] -= waprm; /* H + H2O+ => H2 + OH+ */ c[ih2op-1][ih2op-1] -= wapohp; /* ionized water into water */ c[ih2op-1][ih2op-1] -= wapwa; /* ionized water into OH */ c[ih2op-1][ih2op-1] -= wapoh; /* H2 + OH+ => H + H2O */ c[iohp-1][ih2op-1] += ohpwap; /* processes converting wat to wap */ c[ih2o-1][ih2op-1] += wawap; /* OH + H3P => H2 + H2O+ */ h3pwap = 1.0e-9*hmi.h3plus; c[ioh-1][ih2op-1] += h3pwap; /* H2 + H2O+ => OH + H3+ */ c[ih2op-1][ih2op-1] -= waph3p; /* H2 + H2O+ => H + H3O+ */ waiwai = 6.10e-10*hmi.htwo; c[ih2op-1][ih2op-1] -= waiwai; /* H + H3O+ => H2 + H2O+ */ h3owai = hmrate(6.00e-9,0.390,1.98e4)*h.hi; c[ih3op-1][ih2op-1] += h3owai; /*-------------------------------------------------------------------- * H3O+ balance * H2O + H3P => H2 + H3OP */ c[ih2o-1][ih3op-1] += h3pwai; /* H2 + H3O+ => H2O + H3+ */ c[ih3op-1][ih3op-1] -= h3oh3p; /* H2O + H2+ => H + H3O+ */ c[ih2o-1][ih3op-1] += wah3op; /* H2 + H2O+ => H + H3O+ */ c[ih2op-1][ih3op-1] += waiwai; /* H + H3O+ => H2 + H2O+ */ c[ih3op-1][ih3op-1] -= h3owai; /* H + H3O+ => H2O + H2+ */ c[ih3op-1][ih3op-1] -= h3opwa; /* e- + H3O+ => H2O + H */ c[ih3op-1][ih3op-1] -= e3opwa; /* e- + H3O+ => OH + H3 */ c[ih3op-1][ih3op-1] -= e3opoh; /*-------------------------------------------------------------------- * CH2+ balance * e- + CH2+ => C + H2, there is no carbon balance eqn */ ch2pc = hmrate(2.50e-7,-0.5,0.)*phycon.eden; c[ich2p-1][ich2p-1] -= ch2pc; /* e- + CH2+ => CH + H */ c[ich2p-1][ich2p-1] -= ch2pch; /* H + CH2+ => H2 + CH+ */ c[ich2p-1][ich2p-1] -= h2pchp; /* H2 + CH2+ => CH + H3+ */ c[ich2p-1][ich2p-1] -= ch2ph3; /* O + CH2+ => H + CO+ */ c[ich2p-1][ich2p-1] -= ch2pco; /* H2 + C+ => CH2+ */ cpch2p = hmrate(4.00e-16,-0.20,0.)*hmi.htwo; c[icp-1][ich2p-1] += cpch2p; /* H2 + CH+ => H + CH2+ */ c[ichp-1][ich2p-1] += chphch; /*-------------------------------------------------------------------- * C+ balance * atomic processes changing C+ into Co, calculated in car bal */ c[icp-1][icp-1] -= DestCrt.create[0][5]; /* atomic processes changing Co into C+, calculated in car bal */ c[icatm-1][icp-1] += DestCrt.destroy[0][5]; /* CO + He+ => C + C+ + He */ c[ico-1][icp-1] += cocp; /* H2 + C+ => H + CH+ */ c[icp-1][icp-1] -= cpchp; /* H2 + C+ => CH2+ */ c[icp-1][icp-1] -= cpch2p; /* OH + C+ => H + CO+ */ ohcop = 7.70e-10*hevmolec.hevmol[IPOH-1]; c[icp-1][icp-1] -= ohcop; /* CH+ + H => H2 + C+ */ c[ichp-1][icp-1] += chph; /* O2 + C+ => O + CO+ */ o2cop = 3.76e-10*hevmolec.hevmol[IPOTWO-1]; c[icp-1][icp-1] -= o2cop; /* O2 + C+ => CO + O+ */ o2co = 6.14e-10*hevmolec.hevmol[IPOTWO-1]; c[icp-1][icp-1] -= o2co; /*-------------------------------------------------------------------- * */ if( trace.lgTrace && trace.lgTrMole ) { fprintf( ioQQQ, " COMOLE matrix\n " ); /* turn off lint flag of following evaluation of constant value boolean */ /*lint -e506 */ for( i=0; i < MIN2(NMOLE,8); i++ ) { fprintf( ioQQQ, "%4.4s", chLab[i] ); } fprintf( ioQQQ, " \n" ); for( j=0; j < NMOLE; j++ ) { fprintf( ioQQQ, " %4.4s", chLab[j] ); fprintf( ioQQQ, " " ); for( i=0; i < MIN2(NMOLE,8); i++ ) { fprintf( ioQQQ, "%9.1e", c[i][j] ); } fprintf( ioQQQ, "\n" ); } if( NMOLE >= 9 ) { fprintf( ioQQQ, " COMOLE matrix\n " ); for( i=8; i < NMOLE; i++ ) { fprintf( ioQQQ, "%4.4s", chLab[i] ); } fprintf( ioQQQ, " \n" ); for( j=0; j < NMOLE; j++ ) { fprintf( ioQQQ, " %4.4s", chLab[j] ); fprintf( ioQQQ, " " ); for( i=8; i < NMOLE; i++ ) { fprintf( ioQQQ, "%9.1e", c[i][j] ); } fprintf( ioQQQ, "\n" ); } } fprintf( ioQQQ, " COMOLE H2 den:%10.3e, H2,3+=%10.2e%10.2e Carb sum=%10.3e Oxy sum=%10.3e\n", hmi.htwo, hmi.h2plus, hmi.h3plus, c[NMOLE][icatm-1], c[NMOLE][ioatm-1] ); /* turn back on lint flag of evaluation of constant value boolean */ /*lint +e506 */ } /* invert matrix */ for( i=0; i < NMOLE; i++ ) { for( j=0; j < NMOLE; j++ ) { csav[j][i] = c[j][i]; } } /* which matrix solver? */ if( strcmp(TypMatrx.chMatrix,"matin1 ") == 0 ) { /*matin1();*/ /* this should be set upon error return to non-zero val*/ merror = 1 ; } else if( strcmp(TypMatrx.chMatrix,"linpack") == 0 ) { /* this one may be more robust */ for( j=0; j < NMOLE; j++ ) { for( i=0; i < NMOLE; i++ ) { amat[i][j] = c[i][j]; } bvec[j] = c[NMOLE][j]; } DGETRF(NMOLE,NMOLE,(double*)amat,NMOLE, ipiv,&merror); DGETRS('N',NMOLE,1,(double*)amat,NMOLE,ipiv,bvec,NMOLE, &merror); if( merror != 0 ) { fprintf( ioQQQ, " comole dgetrs error\n" ); puts( "[Stop in comole]" ); exit(1); } /* now put results back into z so rest of code treates only * one case - as if matin1 had been used */ for( i=0; i < NMOLE; i++ ) { c[NMOLE][i] = bvec[i]; } } else if( strcmp(TypMatrx.chMatrix,"veclib ") == 0 ) { /* Jason found this one on the Exemplar, distributed source just stops */ fprintf( ioQQQ, " this has not been checked since H atom conv\n" ); for( j=0; j < NMOLE; j++ ) { for( i=0; i < NMOLE; i++ ) { amat[i][j] = c[i][j]; } bvec[j] = c[NMOLE][j]; } job = 0; rcond = 0.; dgeco((double*)amat,NMOLE,NMOLE,ipiv,rcond,work); dgesl((double*)amat,NMOLE,NMOLE,ipiv,bvec,job); /* if this matrix inversion routine is ever used * put a definition of the error return in the following variable*/ merror = 0; /* now put results back into z so rest of code treates only * one case - as if matin1 had been used */ for( i=0; i < NMOLE; i++ ) { c[NMOLE][i] = bvec[i]; } /* code not in place, so this does not work */ puts( "[Stop in comole]" ); exit(1); } else { fprintf( ioQQQ, " chMatrix type insane in hydrogen, was%7.7s\n", TypMatrx.chMatrix ); puts( "[Stop in comole]" ); exit(1); } /* check for negative populations, which happens when 100% co */ *lgNegPop = FALSE; *lgZerPop = FALSE; for( i=0; i < NMOLE; i++ ) { if( c[NMOLE][i] < 0. ) { *lgNegPop = TRUE; } else if( c[NMOLE][i] == 0. ) { *lgZerPop = TRUE; } } if( merror != 0 ) { fprintf( ioQQQ, " COMOLE matrix inversion error, MERROR=%5ld zone=%5ld\n", merror, ZoneCnt.nzone ); ShowMe(); fprintf( ioQQQ, " Product matrix\n " ); /* turn off lint flag of evaluation of constant value boolean */ /*lint -e506 */ for( i=0; i < MIN2(NMOLE,8); i++ ) { fprintf( ioQQQ, "%4.4s", chLab[i] ); } fprintf( ioQQQ, " \n" ); for( j=0; j < NMOLE; j++ ) { fprintf( ioQQQ, " %4.4s", chLab[j] ); fprintf( ioQQQ, " " ); for( i=0; i < MIN2(NMOLE,8); i++ ) { fprintf( ioQQQ, "%9.1e", csav[i][j]* c[NMOLE][i] ); } fprintf( ioQQQ, "\n" ); } if( NMOLE >= 9 ) { fprintf( ioQQQ, " COMOLE matrix\n " ); for( i=8; i < NMOLE; i++ ) { fprintf( ioQQQ, "%4.4s", chLab[i] ); } fprintf( ioQQQ, " \n" ); for( j=0; j < NMOLE; j++ ) { fprintf( ioQQQ, " %4.4s", chLab[j] ); fprintf( ioQQQ, " " ); for( i=8; i < NMOLE; i++ ) { fprintf( ioQQQ, "%9.1e", csav[i][j]* c[NMOLE][i] ); } fprintf( ioQQQ, "\n" ); } } fprintf( ioQQQ, " Molecular densities relative to H2\n " ); for( j=0; j < NMOLE; j++ ) { fprintf( ioQQQ, "%4.4s:%10.2e", chLab[j] , c[NMOLE][j] ); } fprintf( ioQQQ, " \n" ); puts( "[Stop in comole]" ); exit(1); } if( trace.lgTrace && trace.lgTrMole ) { fprintf( ioQQQ, " Product matrix\n " ); for( i=0; i < MIN2(NMOLE,8); i++ ) { fprintf( ioQQQ, "%4.4s", chLab[i] ); } fprintf( ioQQQ, " \n" ); for( j=0; j < NMOLE; j++ ) { fprintf( ioQQQ, " %4.4s", chLab[j] ); fprintf( ioQQQ, " " ); for( i=0; i < MIN2(NMOLE,8); i++ ) { fprintf( ioQQQ, "%9.1e", csav[i][j]* c[NMOLE][i] ); } fprintf( ioQQQ, "\n" ); } if( NMOLE >= 9 ) { fprintf( ioQQQ, " COMOLE matrix\n " ); for( i=8; i < NMOLE; i++ ) { fprintf( ioQQQ, "%4.4s", chLab[i] ); } fprintf( ioQQQ, " \n" ); for( j=0; j < NMOLE; j++ ) { fprintf( ioQQQ, " %4.4s", chLab[j] ); fprintf( ioQQQ, " " ); for( i=8; i < NMOLE; i++ ) { fprintf( ioQQQ, "%9.1e", csav[i][j]* c[NMOLE][i] ); } fprintf( ioQQQ, "\n" ); } } fprintf( ioQQQ, " Molecular densities relative to H2\n " ); for( j=0; j < NMOLE; j++ ) { fprintf( ioQQQ, "%4.4s:%10.2e", chLab[j] , c[NMOLE][j] ); } fprintf( ioQQQ, " \n" ); } /* turn back on lint flag of evaluation of constant value boolean */ /*lint +e506 */ /* (1,n) is C atomic */ oldcarb = coatom.catmic; coatom.catmic = (float)c[NMOLE][icatm-1]; /* (2,n) is O atomic */ coatom.oatmic = (float)c[NMOLE][ioatm-1]; /* (3,n) is CH */ hevmolec.hevmol[IPCH-1] = (float)c[NMOLE][ich-1]; /* (4,n) is CHP */ hevmolec.hevmol[IPCHP-1] = (float)c[NMOLE][ichp-1]; /* (5,n) is OH */ hevmolec.hevmol[IPOH-1] = (float)c[NMOLE][ioh-1]; /* (6,n) is OHP */ hevmolec.hevmol[IPOHP-1] = (float)c[NMOLE][iohp-1]; /* (9,n) is CO */ hevmolec.hevmol[IPCO-1] = (float)c[NMOLE][ico-1]; /* (10,n) is COP */ hevmolec.hevmol[IPCOP-1] = (float)c[NMOLE][icop-1]; /* (11,n) is H2O */ hevmolec.hevmol[IPH2O-1] = (float)c[NMOLE][ih2o-1]; /* (12,n) is H2OP */ hevmolec.hevmol[IPH2OP-1] = (float)c[NMOLE][ih2op-1]; /* (14,n) is H3OP */ hevmolec.hevmol[IPH3OP-1] = (float)c[NMOLE][ih3op-1]; /* (15,n) is CH2P */ hevmolec.hevmol[IPCH2P-1] = (float)c[NMOLE][ich2p-1]; /* (16,n) is CP */ hevmolec.hevmol[IPCP-1] = (float)c[NMOLE][icp-1]; /* (7,n) is OTWO */ hevmolec.hevmol[IPOTWO-1] = (float)c[NMOLE][io2-1]; /* (8,n) is O2P */ hevmolec.hevmol[IPO2P-1] = (float)c[NMOLE][io2p-1]; /* rememeber largest fraction of carbon in molecules */ cf = (hevmolec.hevmol[IPCH-1] + hevmolec.hevmol[IPCHP-1] + hevmolec.hevmol[IPCO-1] + hevmolec.hevmol[IPCOP-1] + hevmolec.hevmol[IPCH2P-1])/ xIonFracs[0][5]; #if 0 if( hevmolec.hevmol[IPCO-1] / xIonFracs[0][5] > 1. ) { fprintf(ioQQQ," co/c>1, fracs are %g %g %g %g %g\n", hevmolec.hevmol[IPCH-1]/xIonFracs[0][5], hevmolec.hevmol[IPCHP-1]/xIonFracs[0][5], hevmolec.hevmol[IPCO-1]/xIonFracs[0][5], hevmolec.hevmol[IPCOP-1]/xIonFracs[0][5], hevmolec.hevmol[IPCH2P-1]/xIonFracs[0][5] ); } #endif cmfrac.CarMolFrac = (float)MAX2(cmfrac.CarMolFrac,cf); /* heating due to CO photodissociation */ codish.CODissHeat = (float)(phodco*hevmolec.hevmol[IPCO-1]*1e-12); /*666 write(6,'('' codish phodco, GammaHabing'',1p,3e10.2)') CODissHeat, * 1 phodco , GammaHabing */ HeatingCom.heating[9][0] = codish.CODissHeat; /* cooling due to coll excit of vib-rot levels * from Hollenbach and McKee 1979 */ crden = 1.9e3*TePowers.sqrte; /* this is their cooling coefficent, per H2, and CO */ covib.COVibCool = (float)(1.94e-23*phycon.te*phycon.te/(hmi.htwo + crden + 1.5*pow(hmi.htwo*crden,0.5))); /* this is total cooling per unit vol */ /* multiply by smaller of co or total carbon abundance (all phases). There are models * where ALL c is in co and matrix is unstable, returning co/c > 1. other molecules * have negative abundances. this will eventually be detected by the code, which * will bail out, but we don't get to that stage if co cooling is oscillating. */ covib.COVibCool *= (float)(hmi.htwo*MIN2(hevmolec.hevmol[IPCO-1],xIonFracs[0][5])*0.3); /* its derivative */ covib.dCOVdT = (float)(covib.COVibCool*2./phycon.te); /* check for negative or zero atomic carbon populations */ if( coatom.catmic <= 0. ) { coatom.catmic = (float)(fabs(coatom.catmic/2.)); coatom.catmic = oldcarb; if( hevmolec.hevmol[IPCO-1] > hevmolec.hevmol[IPCH-1] ) { hevmolec.hevmol[IPCO-1] = (float)MAX2(hevmolec.hevmol[IPCO-1]/ 2.,hevmolec.hevmol[IPCO-1]-coatom.catmic); } else { hevmolec.hevmol[IPCH-1] = (float)MAX2(hevmolec.hevmol[IPCH-1]/ 2.,hevmolec.hevmol[IPCH-1]-coatom.catmic); } } if( coatom.oatmic <= 0. ) { /* WRITE(QQ,'('' WARNING; O population is negative or 0.'',1P, * 1 E10.2)') OATMIC */ coatom.oatmic = 0.; } if( hevmolec.hevmol[IPCH-1] <= 0. ) { /* WRITE(QQ,'('' WARNING; CH population is negative or 0.'',1P, * 1 E10.2)') CH */ hevmolec.hevmol[IPCH-1] = -hevmolec.hevmol[IPCH-1]/2.f; } if( hevmolec.hevmol[IPCHP-1] <= 0. ) { /* WRITE(QQ,'('' WARNING; CH+ population is negative or 0.'',1P, * 1 E10.2)') CHP */ hevmolec.hevmol[IPCHP-1] = -hevmolec.hevmol[IPCHP-1]/2.f; } if( hevmolec.hevmol[IPOH-1] <= 0. ) { /* WRITE(QQ,'('' WARNING; OH population is negative or 0.'',1P, * 1 E10.2)') OH */ hevmolec.hevmol[IPOH-1] = 0.; } if( hevmolec.hevmol[IPOHP-1] <= 0. ) { /* WRITE(QQ,'('' WARNING; OH+ population is negative or 0.'',1P, * 1 E10.2)') OHP */ hevmolec.hevmol[IPOHP-1] = 0.; } if( hevmolec.hevmol[IPCO-1] <= 0. ) { /* WRITE(QQ,'('' WARNING; CO population is negative or 0.'',1P, * 1 E10.2)') CO */ hevmolec.hevmol[IPCO-1] = 0.; } if( hevmolec.hevmol[IPCOP-1] <= 0. ) { /* WRITE(QQ,'('' WARNING; CO+ population is negative or 0.'',1P, * 1 E10.2)') COP */ hevmolec.hevmol[IPCOP-1] = 0.; } if( hevmolec.hevmol[IPH2O-1] <= 0. ) { hevmolec.hevmol[IPH2O-1] = 0.; } if( hevmolec.hevmol[IPH2OP-1] <= 0. ) { hevmolec.hevmol[IPH2OP-1] = 0.; } if( hevmolec.hevmol[IPH3OP-1] <= 0. ) { hevmolec.hevmol[IPH3OP-1] = 0.; } if( hevmolec.hevmol[IPCH2P-1] <= 0. ) { hevmolec.hevmol[IPCH2P-1] = -hevmolec.hevmol[IPCH2P-1]/2.f; } if( hevmolec.hevmol[IPCP-1] <= 0. ) { hevmolec.hevmol[IPCP-1] = -hevmolec.hevmol[IPCP-1]/2.f; } } # ifdef DEBUG_FUN fputs( " <->comole()\n", debug_fp ); # endif return; } /*lint +e550 */ /*molav average old and new molecular equilibrium balance from comole */ void molav(long int JobMolAv) { long int i; double damper; static double hevsav[NHEVML]; # ifdef DEBUG_FUN fputs( "<+>molav()\n", debug_fp ); # endif /* when called with IX=1, set saver */ if( JobMolAv == 1 ) { /* this is first call - remember initial abundances */ for( i=0; i < NHEVML; i++ ) { hevsav[i] = hevmolec.hevmol[i]; } if( trace.lgTrace && trace.lgTrMole ) { fprintf( ioQQQ, " MOLAV pop1=" ); for(i=0; i < NHEVML; i++) fprintf( ioQQQ, "%9.1e", hevmolec.hevmol[i] ); fprintf( ioQQQ, "\n" ); } } else if( JobMolAv < 0 ) { for( i=0; i < NHEVML; i++ ) { hevmolec.hevmol[i] = 0.; hevsav[i] = 0.; } } else { /* get new numbers - take average of old and new */ for( i=0; i < NHEVML; i++ ) { /* parameter to mix old and new */ damper = 0.2; hevmolec.hevmol[i] = (float)(hevmolec.hevmol[i]*(1. - damper) + hevsav[i]*damper); hevsav[i] = hevmolec.hevmol[i]; } if( trace.lgTrace && trace.lgTrMole ) { fprintf( ioQQQ, " MOLAV pop2=" ); for(i=0; i < NHEVML; i++) fprintf( ioQQQ, "%9.1e", hevmolec.hevmol[i] ); fprintf( ioQQQ, "\n" ); } } # ifdef DEBUG_FUN fputs( " <->molav()\n", debug_fp ); # endif return; } /*=================================================================*/ /*codriv main driver for heavy molecular equilibrium routines */ #define COTOLER 0.10 #define LUPMAX 10 void codriv(void) { int lgNegPop, lgZerPop; static long int loop; static float ohold; # ifdef DEBUG_FUN fputs( "<+>codriv()\n", debug_fp ); # endif ohold = hevmolec.hevmol[IPOH-1]; comole(&lgNegPop,&lgZerPop); molav(1); /* >>chng 97 july 29, turn on flag for too much CO here too, there were * models in which this test never took place since oh went zero */ if( hevmolec.hevmol[IPOH-1] == 0. ) { if( hevmolec.hevmol[IPCO-1]/xIonFracs[0][5] > colimt.COLimit ) lotsco.lgLotsCO = TRUE; # ifdef DEBUG_FUN fputs( " <->codriv()\n", debug_fp ); # endif return; } loop = 1; while( ((fabs(ohold/hevmolec.hevmol[IPOH-1]-1.) > COTOLER || lgNegPop) && (loop < LUPMAX)) && (hevmolec.hevmol[IPCO-1]/xIonFracs[0][5] < colimt.COLimit) ) { molav(2); ohold = hevmolec.hevmol[IPOH-1]; comole(&lgNegPop,&lgZerPop); if( hevmolec.hevmol[IPOH-1] == 0. ) { # ifdef DEBUG_FUN fputs( " <->codriv()\n", debug_fp ); # endif return; } loop += 1; } /* this sets flag if there is too much CO, a precursor to negative pops, * and so a stopping condition */ if( hevmolec.hevmol[IPCO-1]/xIonFracs[0][5] > colimt.COLimit ) { lotsco.lgLotsCO = TRUE; } /* did the molecule network have negative pops? */ else if( lgNegPop ) { lotsco.lgLotsCO = TRUE; lgNegPop = FALSE; if( called.lgTalk ) { fprintf( ioQQQ, " CODRIV failed to converge, zone=%4ld, CO/C=%10.2e negpop?%1c\n", ZoneCnt.nzone, hevmolec.hevmol[IPCO-1]/xIonFracs[0][5], TorF(lgNegPop) ); } } /* this test, hit upper limit to number of iterations */ else if( loop == LUPMAX ) { if( called.lgTalk ) { fprintf( ioQQQ, " CODRIV failed to converge, zone=%4ld, CO/C=%10.2e negpop?%1c\n", ZoneCnt.nzone, hevmolec.hevmol[IPCO-1]/xIonFracs[0][5], TorF(lgNegPop) ); } } if( lgNegPop ) comneg.lgComNeg = TRUE; # ifdef DEBUG_FUN fputs( " <->codriv()\n", debug_fp ); # endif return; }