00001
00002
00003
00004 #include "cddefines.h"
00005 #include "taulines.h"
00006 #include "dense.h"
00007 #include "ionbal.h"
00008 #include "thermal.h"
00009 #include "phycon.h"
00010 #include "hmi.h"
00011 #include "dynamics.h"
00012 #include "conv.h"
00013 #include "trace.h"
00014 #include "coolheavy.h"
00015 #include "timesc.h"
00016 #include "thirdparty.h"
00017 #include "mole.h"
00018 #include "mole_co_priv.h"
00019 #include "mole_co_atom.h"
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036 void CO_solve(
00037
00038 bool *lgNegPop,
00039
00040 bool *lgZerPop )
00041 {
00042
00043 int32 merror;
00044 long int i, j, k, n,
00045 nelem , ion , ion2;
00046 double
00047 co_denominator;
00048 float cartot_mol, oxytot_mol;
00049 float abundan;
00050 struct chem_element_s *element;
00051 double sum;
00052
00053 CO_step();
00054
00055
00056 for(i=0;i<mole.num_comole_calc;i++)
00057 {
00058 sum = 0.;
00059 for(j=0;j<mole.num_comole_calc;j++)
00060 {
00061 sum = sum+fabs(mole.c[i][j]);
00062 }
00063 if(sum < SMALLFLOAT && fabs(mole.b[i]) < SMALLFLOAT)
00064 {
00065 mole.c[i][i] = 1.;
00066 }
00067 }
00068
00069 cartot_mol = dense.xMolecules[ipCARBON] + findspecies("C")->hevmol + findspecies("C+")->hevmol;
00070 oxytot_mol = dense.xMolecules[ipOXYGEN] + findspecies("O")->hevmol + findspecies("O+")->hevmol;
00071
00072 ASSERT( cartot_mol >= 0. && oxytot_mol >= 0.);
00073
00074 *lgNegPop = false;
00075 *lgZerPop = false;
00076
00077
00078 for(nelem=ipLITHIUM; nelem < LIMELM; ++nelem)
00079 {
00080 for(ion=0; ion < nelem+2; ++ion)
00081 {
00082
00083 mole.sink[nelem][ion] = 0.;
00084 mole.source[nelem][ion] = 0.;
00085
00086 for(ion2=0; ion2< nelem+2; ++ion2)
00087 {
00088 mole.xMoleChTrRate[nelem][ion][ion2] = 0.;
00089 }
00090 }
00091 }
00092
00093
00094
00095
00096
00097
00098 if( iteration >= 2 && dynamics.lgAdvection )
00099 {
00100 for(i=0;i<mole.num_comole_calc;i++)
00101 {
00102 if(COmole[i]->n_nuclei > 1) {
00103 mole.c[i][i] -= dynamics.Rate;
00104
00105
00106
00107
00108 mole.b[i] -= (COmole[i]->hevmol*dynamics.Rate-dynamics.CO_molec[i]);
00109 }
00110 }
00111 }
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140 for(i=0;i<mole.num_elements; ++i)
00141 {
00142 element = chem_element[i];
00143 mole.sink[element->ipCl][0] -= (mole.c[element->ipMl][element->ipMl] + mole.c[element->ipMl][element->ipMlP])*dense.lgElmtOn[element->ipCl];
00144 mole.sink[element->ipCl][1] -= (mole.c[element->ipMlP][element->ipMlP] + mole.c[element->ipMlP][element->ipMl] )*dense.lgElmtOn[element->ipCl];
00145 }
00146
00147
00148
00149 for(j=0;j < mole.num_elements; ++j)
00150 {
00151 element = chem_element[j];
00152 for(i=0; i < mole.num_comole_calc; ++i)
00153 {
00154 mole.source[element->ipCl][1] += COmole[i]->hevmol*mole.c[i][element->ipMlP]*dense.lgElmtOn[element->ipCl];
00155 mole.source[element->ipCl][0] += COmole[i]->hevmol*mole.c[i][element->ipMl]*dense.lgElmtOn[element->ipCl];
00156 }
00157 }
00158
00159
00160
00161
00162 for(i=0;i < mole.num_elements; ++i)
00163 {
00164 element = chem_element[i];
00165 mole.source[element->ipCl][1] -= ( mole.c[element->ipMlP][element->ipMlP]*COmole[element->ipMlP]->hevmol +
00166 mole.c[element->ipMl][element->ipMlP]*COmole[element->ipMl]->hevmol)*dense.lgElmtOn[element->ipCl];
00167
00168 mole.source[element->ipCl][0] -= ( mole.c[element->ipMl][element->ipMl]*COmole[element->ipMl]->hevmol +
00169 mole.c[element->ipMlP][element->ipMl]*COmole[element->ipMlP]->hevmol)*dense.lgElmtOn[element->ipCl];
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182 mole.source[element->ipCl][1] = mole.source[element->ipCl][1] - mole.b[element->ipMlP];
00183 mole.source[element->ipCl][0] = mole.source[element->ipCl][0] - mole.b[element->ipMl];
00184
00185 }
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195 for(i=2; i < LIMELM; ++i)
00196 {
00197 for(j=0; j < 2; ++j)
00198 {
00199
00200 if(mole.source[i][j] < 0)
00201 {
00202 mole.sink[i][j] -= (mole.source[i][j]/SDIV(dense.xIonDense[i][j]));
00203 mole.source[i][j] = 0;
00204 }
00205 }
00206 }
00207
00208
00209
00210 for(i=0;i < mole.num_elements; ++i)
00211 {
00212 element = chem_element[i];
00213 mole.xMoleChTrRate[element->ipCl][1][0] = (float)(mole.c[element->ipMlP][element->ipMl] -
00214 ionbal.RateRecomTot[element->ipCl][0])*dense.lgElmtOn[element->ipCl];
00215
00216 mole.xMoleChTrRate[element->ipCl][0][1] = (float)(mole.c[element->ipMl][element->ipMlP] -
00217 ionbal.RateIonizTot[element->ipCl][0])*dense.lgElmtOn[element->ipCl];
00218 }
00219
00220
00221 for(i=2; i < LIMELM; ++i)
00222 {
00223 for(j=0; j < 2; ++j)
00224 {
00225 for(k=0; k< 2; ++k)
00226 {
00227
00228 if(j != k)
00229 {
00230 if( mole.xMoleChTrRate[i][j][k] < 0. )
00231 {
00232 mole.xMoleChTrRate[i][k][j] -= mole.xMoleChTrRate[i][j][k];
00233 mole.xMoleChTrRate[i][j][k] = 0.;
00234 }
00235 }
00236 }
00237 }
00238 }
00239
00240 for(n=0;n<mole.num_elements;n++)
00241 {
00242 tot_ion[n] = dense.gas_phase[chem_element[n]->ipCl];
00243 for( i=2; i < chem_element[n]->ipCl+2; i++ )
00244 {
00245 tot_ion[n] -= dense.xIonDense[chem_element[n]->ipCl][i];
00246 }
00247 }
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262 if( iteration >= 2 && dynamics.lgAdvection)
00263 fixit();
00264
00265 for(n=0;n<mole.num_elements;n++)
00266 {
00267 mole.b[chem_element[n]->ipMl] = tot_ion[n];
00268 }
00269
00270
00271
00272
00273
00274
00275
00276
00277 for( j=0; j < mole.num_comole_calc; j++ )
00278 {
00279 for(i=0;i<mole.num_elements;i++)
00280 {
00281 element = chem_element[i];
00282 mole.c[j][element->ipMl] = COmole[j]->nElem[element->ipCl];
00283 }
00284 }
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309 if( trace.lgTrace && trace.lgTr_CO_Mole )
00310 {
00311 fprintf( ioQQQ, " COMOLE matrix\n " );
00312
00313 # if 0
00314 for( i=0; i < MIN2(mole.num_comole_calc,8); i++ )
00315 {
00316 fprintf( ioQQQ, "%4.4s", COmole[i]->label );
00317 }
00318 fprintf( ioQQQ, " \n" );
00319
00320 for( j=0; j < mole.num_comole_calc; j++ )
00321 {
00322 fprintf( ioQQQ, " %4.4s", COmole[j]->label );
00323 fprintf( ioQQQ, " " );
00324 for( i=0; i < MIN2(mole.num_comole_calc,8); i++ )
00325 {
00326 fprintf( ioQQQ, "%9.1e", mole.c[i][j] );
00327 }
00328 fprintf( ioQQQ, "\n" );
00329 }
00330 # endif
00331
00332 fprintf( ioQQQ, " COMOLE matrix\n " );
00333
00334 for( i=0; i < mole.num_comole_calc; i++ )
00335 {
00336 fprintf( ioQQQ, "%4.4s", COmole[i]->label );
00337 }
00338 fprintf( ioQQQ, " \n" );
00339
00340 for( j=0; j < mole.num_comole_calc; j++ )
00341 {
00342 fprintf( ioQQQ, " %4.4s", COmole[j]->label );
00343 fprintf( ioQQQ, " " );
00344 for( i=0; i < mole.num_comole_calc; i++ )
00345 {
00346
00347 fprintf( ioQQQ, ", %.15e", mole.c[i][j] );
00348 }
00349 fprintf( ioQQQ, "\n" );
00350 }
00351
00352 fprintf( ioQQQ, " COMOLE b\n " );
00353 for( j=0; j < mole.num_comole_calc; j++ )
00354 {
00355 fprintf( ioQQQ, " %4.4s", COmole[j]->label );
00356 fprintf( ioQQQ, ", %.15e\n",mole.b[j] );
00357 }
00358
00359 #if 0
00360 fprintf( ioQQQ,
00361 " COMOLE H2 den:%10.3e, H2,3+=%10.2e%10.2e Carb sum=%10.3e Oxy sum=%10.3e\n",
00362 hmi.H2_total,
00363 hmi.Hmolec[ipMH2p],
00364 hmi.Hmolec[ipMH3p],
00365 mole.c[mole.num_heavy_molec][findspecies("C")->index],
00366 mole.c[mole.num_heavy_molec][findspecies("O")->index] );
00367
00368 #endif
00369 }
00370
00371
00372 for( j=0; j < mole.num_comole_calc; j++ )
00373 {
00374 if( -mole.c[j][j] > SMALLFLOAT )
00375 {
00376
00377 timesc.AgeCOMoleDest[j] = -1./mole.c[j][j];
00378
00379
00380 }
00381 else
00382 {
00383 timesc.AgeCOMoleDest[j] = 0.;
00384 }
00385 }
00386
00387
00388 for( j=0; j < mole.num_comole_calc; j++ )
00389 {
00390 for( i=0; i < mole.num_comole_calc; i++ )
00391 {
00392 mole.amat[i][j] = mole.c[i][j];
00393 }
00394 }
00395
00396 merror = 0;
00397
00398
00399 getrf_wrapper(mole.num_comole_calc,mole.num_comole_calc,mole.amat[0],mole.num_comole_calc, ipiv,&merror);
00400
00401
00402 if( merror != 0 )
00403 {
00404 fprintf( ioQQQ, " CO_solve getrf_wrapper error\n" );
00405 puts( "[Stop in CO_solve]" );
00406 cdEXIT(EXIT_FAILURE);
00407 }
00408
00409 getrs_wrapper('N',mole.num_comole_calc,1,mole.amat[0],mole.num_comole_calc,ipiv,mole.b,mole.num_comole_calc,&merror);
00410
00411
00412 if( merror != 0 )
00413 {
00414 fprintf( ioQQQ, " CO_solve: dgetrs finds singular or ill-conditioned matrix\n" );
00415 puts( "[Stop in CO_solve]" );
00416 cdEXIT(EXIT_FAILURE);
00417 }
00418
00419
00420 *lgNegPop = false;
00421 *lgZerPop = false;
00422 co_denominator = 0;
00423 # if 0
00424 if(fnzone > 1.02)
00425 {
00426 for( i=0; i < mole.num_comole_calc; i++)
00427 {
00428
00429 for( j=61; j < 62; j++ )
00430 {
00431
00432 printf("ZONE\t%.5e\tSPECIES_1\t%s\tSPECIES_2\t%s\tDEST_RATE\t%.3e\tbvec\t%.3e\n",
00433 fnzone, COmole[i]->label,COmole[j]->label,mole.c[i][j],mole.b[i]);
00434
00435 }
00436 }
00437 }
00438 # endif
00439 for( i=0; i < mole.num_comole_calc; i++ )
00440 {
00441
00442 if( mole.b[i] < 0. )
00443 {
00444
00445
00446
00447
00448
00449
00450
00451
00452
00453
00454
00455
00456
00457
00458
00459
00460
00461
00462
00463
00464
00465
00466
00467
00468
00469
00470 if( dense.lgElmtOn[COmole[i]->nelem_hevmol] )
00471 {
00472 co_denominator = dense.gas_phase[COmole[i]->nelem_hevmol];
00473 }
00474 else
00475 {
00476
00477 co_denominator = 1e10;
00478 }
00479
00480
00481
00482 ASSERT( co_denominator > SMALLFLOAT || !dense.lgElmtOn[COmole[i]->nelem_hevmol] );
00483
00484
00485
00486
00487
00488
00489
00490 if( mole.b[i] / SDIV(co_denominator) > -1e-5 )
00491 {
00492
00493
00494
00495 mole.b[i] *= -1.;
00496 }
00497 else
00498 {
00499
00500 *lgNegPop = true;
00501
00502
00503
00504
00505
00506 if( nzone>co.co_nzone )
00507 {
00508 static long int nzFail=-2;
00509 long int nzFail2 = (long)fnzone*100;
00510
00511 if( !conv.lgSearch )
00512 fprintf(ioQQQ," PROBLEM ");
00513 fprintf(ioQQQ,
00514 " CO_solve neg pop for species %li %s, value is %.2e rel value is %.2e zone %.2f Te %.4e Search?%c\n",
00515 i ,
00516 COmole[i]->label ,
00517
00518 mole.b[i],
00519
00520 mole.b[i] / SDIV(co_denominator) ,
00521 fnzone,
00522 phycon.te,TorF( conv.lgSearch ) );
00523
00524
00525 if( nzFail2 !=nzFail )
00526 {
00527 nzFail = nzFail2;
00528 ConvFail( "pops" , "CO");
00529 }
00530 }
00531 mole.b[i] *= -1;
00532
00533
00534
00535
00536
00537 }
00538 }
00539 else if( mole.b[i] == 0. )
00540 {
00541
00542
00543 *lgZerPop = true;
00544
00545
00546 # if 0
00547 if( CODEBUG>1 )
00548 {
00549 fprintf(ioQQQ," PROBLEM ");
00550 fprintf(ioQQQ,
00551 " CO_solve zero pop for species %li %s, value is %.2e zone %li\n",
00552 i , COmole[i]->label , mole.b[i], nzone);
00553 }
00554 # endif
00555 mole.b[i] = SMALLFLOAT;
00556 }
00557 COmole[i]->hevmol = (float)mole.b[i];
00558
00559 ASSERT( !isnan( COmole[i]->hevmol ) );
00560 }
00561
00562
00563
00564
00565
00566
00567
00568
00569 if( *lgNegPop && (nzone>co.co_nzone &&!conv.lgSearch) )
00570 {
00571 conv.lgConvPops = false;
00572 fprintf(ioQQQ," CO network negative population occurred, Te=%.4e, calling ConvFail. ",
00573 phycon.te);
00574 fprintf( ioQQQ, " CO/C=%9.1e", findspecies("CO")->hevmol/SDIV(dense.gas_phase[ipCARBON]) );
00575 fprintf( ioQQQ, "\n" );
00576
00577 *lgNegPop = false;
00578 }
00579
00580
00581 if( 0 && *lgNegPop )
00582 {
00583
00584 for(j=0; j<2; ++j )
00585 {
00586 double sumcar = 0.;
00587 double sumoxy = 0.;
00588 double renorm;
00589 for( i=0; i < mole.num_comole_calc; i++ )
00590 {
00591
00592 sumcar += COmole[i]->hevmol*COmole[i]->nElem[ipCARBON];
00593 sumoxy += COmole[i]->hevmol*COmole[i]->nElem[ipOXYGEN];
00594 }
00595 renorm = (tot_ion[0] + tot_ion[1]) / SDIV( sumcar+sumoxy);
00596 if(j)
00597 fprintf(ioQQQ,"\t%f\n", renorm);
00598 else
00599 fprintf(ioQQQ,"renormco\t%.2f\t%f", fnzone, renorm);
00600 for( i=0; i < mole.num_comole_calc; i++ )
00601 {
00602 COmole[i]->hevmol *= (float)renorm;
00603
00604 ASSERT( !isnan( COmole[i]->hevmol ) );
00605 }
00606 }
00607 }
00608
00609 if( merror != 0 )
00610 {
00611 fprintf( ioQQQ, " COMOLE matrix inversion error, MERROR=%5ld zone=%5ld\n",
00612 (long)merror, nzone );
00613 ShowMe();
00614 fprintf( ioQQQ, " Product matrix\n " );
00615
00616 for( i=0; i < MIN2(mole.num_comole_calc,8); i++ )
00617 {
00618 fprintf( ioQQQ, "%4.4s", COmole[i]->label );
00619 }
00620 fprintf( ioQQQ, " \n" );
00621
00622 for( j=0; j < mole.num_comole_calc; j++ )
00623 {
00624 fprintf( ioQQQ, " %4.4s", COmole[j]->label );
00625 fprintf( ioQQQ, " " );
00626
00627 for( i=0; i < MIN2(mole.num_comole_calc,8); i++ )
00628 {
00629 fprintf( ioQQQ, "%9.1e", mole.amat[i][j]*
00630 COmole[i]->hevmol );
00631 }
00632 fprintf( ioQQQ, "\n" );
00633 }
00634
00635 if( mole.num_comole_calc >= 9 )
00636 {
00637 fprintf( ioQQQ, " COMOLE matrix\n " );
00638 for( i=8; i < mole.num_comole_calc; i++ )
00639 {
00640 fprintf( ioQQQ, "%4.4s", COmole[i]->label );
00641 }
00642 fprintf( ioQQQ, " \n" );
00643
00644 for( j=0; j < mole.num_comole_calc; j++ )
00645 {
00646 fprintf( ioQQQ, " %4.4s", COmole[j]->label );
00647 fprintf( ioQQQ, " " );
00648 for( i=8; i < mole.num_comole_calc; i++ )
00649 {
00650 fprintf( ioQQQ, "%9.1e",
00651 mole.amat[i][j]* COmole[i]->hevmol );
00652 }
00653 fprintf( ioQQQ, "\n" );
00654 }
00655 }
00656
00657 fprintf( ioQQQ, " Mole dens:" );
00658 for( j=0; j < mole.num_comole_calc; j++ )
00659 {
00660 fprintf( ioQQQ, " %4.4s:%.2e", COmole[j]->label
00661 , COmole[j]->hevmol );
00662 }
00663 fprintf( ioQQQ, " \n" );
00664
00665 puts( "[Stop in CO_solve]" );
00666 cdEXIT(EXIT_FAILURE);
00667 }
00668
00669 if( trace.lgTrace && trace.lgTr_CO_Mole )
00670 {
00671 fprintf( ioQQQ, " Product matrix\n " );
00672
00673 for( i=0; i < MIN2(mole.num_comole_calc,8); i++ )
00674 {
00675 fprintf( ioQQQ, "%4.4s", COmole[i]->label );
00676 }
00677 fprintf( ioQQQ, " \n" );
00678
00679 for( j=0; j < mole.num_comole_calc; j++ )
00680 {
00681 fprintf( ioQQQ, " %4.4s", COmole[j]->label );
00682 fprintf( ioQQQ, " " );
00683 for( i=0; i < MIN2(mole.num_comole_calc,8); i++ )
00684 {
00685 fprintf( ioQQQ, "%9.1e", mole.amat[i][j]*
00686 COmole[i]->hevmol );
00687 }
00688 fprintf( ioQQQ, "\n" );
00689
00690 }
00691
00692 if( mole.num_comole_calc >= 9 )
00693 {
00694 fprintf( ioQQQ, " COMOLE matrix\n " );
00695 for( i=8; i < mole.num_comole_calc; i++ )
00696 {
00697 fprintf( ioQQQ, "%4.4s", COmole[i]->label );
00698 }
00699 fprintf( ioQQQ, " \n" );
00700
00701 for( j=0; j < mole.num_comole_calc; j++ )
00702 {
00703 fprintf( ioQQQ, " %4.4s", COmole[j]->label );
00704 fprintf( ioQQQ, " " );
00705
00706 for( i=8; i < mole.num_comole_calc; i++ )
00707 {
00708 fprintf( ioQQQ, "%9.1e", mole.amat[i][j]* COmole[i]->hevmol );
00709 }
00710 fprintf( ioQQQ, "\n" );
00711 }
00712 }
00713
00714 fprintf( ioQQQ, " Mole dens:" );
00715 for( j=0; j < mole.num_comole_calc; j++ )
00716 {
00717 fprintf( ioQQQ, " %4.4s:%.2e", COmole[j]->label
00718 , COmole[j]->hevmol );
00719 }
00720 fprintf( ioQQQ, " \n" );
00721 }
00722
00723
00724 co.CODissHeat = (float)(CO_findrate("PHOTON,CO=>C,O")*1e-12);
00725
00726 thermal.heating[0][9] = co.CODissHeat;
00727
00728
00729 abundan = findspecies("CO")->hevmol;
00730
00731 ASSERT( C12O16Rotate[0].IonStg < LIMELM+2 );
00732 ASSERT( C12O16Rotate[0].nelem-1 < LIMELM+2 );
00733 dense.xIonDense[ C12O16Rotate[0].nelem-1][C12O16Rotate[0].IonStg-1] = abundan;
00734
00735 CO_PopsEmisCool(&C12O16Rotate, nCORotate,abundan , "12CO",
00736 &CoolHeavy.C12O16Rot,&CoolHeavy.dC12O16Rot );
00737
00738
00739
00740
00741 abundan = findspecies("CO")->hevmol/co.RatioC12O16_2_C13O16;
00742
00743 ASSERT( C13O16Rotate[0].IonStg < LIMELM+2 );
00744 ASSERT( C13O16Rotate[0].nelem-1 < LIMELM+2 );
00745 dense.xIonDense[ C13O16Rotate[0].nelem-1][C13O16Rotate[0].IonStg-1] = abundan;
00746
00747 CO_PopsEmisCool(&C13O16Rotate, nCORotate,abundan ,"13CO",
00748 &CoolHeavy.C13O16Rot,&CoolHeavy.dC13O16Rot );
00749
00750
00751 for( i=0;i<LIMELM; ++i )
00752 {
00753 dense.xMolecules[i] = 0.;
00754 }
00755
00756
00757
00758 for(i=0;i<N_H_MOLEC;i++)
00759 {
00760 dense.xMolecules[ipHYDROGEN] += hmi.Hmolec[i]*hmi.nProton[i];
00761 }
00762 dense.xMolecules[ipHYDROGEN] -= (hmi.Hmolec[ipMH] + hmi.Hmolec[ipMHp]);
00763
00764
00765
00766
00767
00768
00769
00770 dense.H_sum_in_CO = 0;
00771 for(i=0; i<mole.num_comole_calc; ++i)
00772 {
00773 if(COmole[i]->n_nuclei <= 1)
00774 continue;
00775
00776 dense.H_sum_in_CO += COmole[i]->nElem[ipHYDROGEN]*COmole[i]->hevmol;
00777
00778 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
00779 {
00780 if( mole.lgElem_in_chemistry[nelem] )
00781 {
00782 dense.xMolecules[nelem] +=COmole[i]->nElem[nelem]*COmole[i]->hevmol;
00783 }
00784 }
00785 }
00786
00787
00788
00789
00790
00791
00792
00793
00794
00795
00796
00797
00798
00799 for(i=0;i<mole.num_elements; ++i)
00800 {
00801 element = chem_element[i];
00802
00803 COmole[element->ipMlP]->hevmol = ((dense.xIonDense[element->ipCl][1]+ COmole[element->ipMlP]->hevmol)/2)*dense.lgElmtOn[element->ipCl];
00804
00805 COmole[element->ipMl]->hevmol = ((dense.xIonDense[element->ipCl][0]+ COmole[element->ipMl]->hevmol)/2)*dense.lgElmtOn[element->ipCl];
00806
00807
00808 ASSERT( !isnan( COmole[element->ipMlP]->hevmol ) );
00809 ASSERT( !isnan( COmole[element->ipMl]->hevmol ) );
00810 }
00811
00812
00813 if( dense.lgElmtOn[ipCARBON] &&
00814 fabs(dense.xIonDense[ipCARBON][0]- findspecies("C")->hevmol)/SDIV(dense.gas_phase[ipCARBON]) >0.001 )
00815 {
00816 conv.lgConvIoniz = false;
00817 sprintf( conv.chConvIoniz, "CO C0 con");
00818 conv.BadConvIoniz[0] = dense.xIonDense[ipCARBON][0];
00819 conv.BadConvIoniz[1] = findspecies("C")->hevmol;
00820 }
00821 else if( dense.lgElmtOn[ipCARBON] &&
00822 fabs(dense.xIonDense[ipCARBON][1]- findspecies("C+")->hevmol)/SDIV(dense.gas_phase[ipCARBON]) >0.001 )
00823 {
00824 conv.lgConvIoniz = false;
00825 sprintf( conv.chConvIoniz, "CO C1 con");
00826 conv.BadConvIoniz[0] = dense.xIonDense[ipCARBON][1];
00827 conv.BadConvIoniz[1] = findspecies("C+")->hevmol;
00828 }
00829 else if( dense.lgElmtOn[ipOXYGEN] &&
00830 fabs(dense.xIonDense[ipOXYGEN][0]- findspecies("O")->hevmol)/SDIV(dense.gas_phase[ipOXYGEN]) >0.001 )
00831 {
00832 conv.lgConvIoniz = false;
00833 sprintf( conv.chConvIoniz, "CO O0 con");
00834 conv.BadConvIoniz[0] = dense.xIonDense[ipOXYGEN][0];
00835 conv.BadConvIoniz[1] = findspecies("O")->hevmol;
00836 }
00837 else if( dense.lgElmtOn[ipOXYGEN] &&
00838 fabs(dense.xIonDense[ipOXYGEN][1]- findspecies("O+")->hevmol)/SDIV(dense.gas_phase[ipOXYGEN]) >0.001 )
00839 {
00840 conv.lgConvIoniz = false;
00841 sprintf( conv.chConvIoniz, "CO O1 con");
00842 conv.BadConvIoniz[0] = dense.xIonDense[ipOXYGEN][1];
00843 conv.BadConvIoniz[1] = findspecies("O+")->hevmol;
00844 }
00845
00846
00847 for(i=0;i<mole.num_elements; ++i)
00848 {
00849 element = chem_element[i];
00850 dense.xIonDense[element->ipCl][0] = COmole[element->ipMl]->hevmol*dense.lgElmtOn[element->ipCl];
00851 dense.xIonDense[element->ipCl][1] = COmole[element->ipMlP]->hevmol*dense.lgElmtOn[element->ipCl];
00852 dense.IonLow[element->ipCl] = 0;
00853 dense.IonHigh[element->ipCl] = MAX2( dense.IonHigh[element->ipCl] , 1 );
00854 }
00855
00856
00857 # define EPS_MOLE 0.1
00858 if( dense.lgElmtOn[ipHYDROGEN] &&
00859 dense.xMolecules[ipHYDROGEN] > dense.gas_phase[ipHYDROGEN]*(1.+EPS_MOLE) )
00860 {
00861 conv.lgConvIoniz = false;
00862 sprintf( conv.chConvIoniz, "COcon%2i",ipHYDROGEN );
00863 conv.BadConvIoniz[0] = dense.xMolecules[ipHYDROGEN];
00864 conv.BadConvIoniz[1] = dense.gas_phase[ipHYDROGEN];
00865 }
00866 else if( dense.lgElmtOn[ipCARBON] &&
00867 dense.xMolecules[ipCARBON] > dense.gas_phase[ipCARBON]*(1.+EPS_MOLE) )
00868 {
00869 conv.lgConvIoniz = false;
00870 sprintf( conv.chConvIoniz, "COcon%2i",ipCARBON );
00871 conv.BadConvIoniz[0] = dense.xMolecules[ipCARBON];
00872 conv.BadConvIoniz[1] = dense.gas_phase[ipCARBON];
00873 }
00874 else if( dense.lgElmtOn[ipOXYGEN] &&
00875 dense.xMolecules[ipOXYGEN] > dense.gas_phase[ipOXYGEN]*(1.+EPS_MOLE) )
00876 {
00877 conv.lgConvIoniz = false;
00878 sprintf( conv.chConvIoniz, "COcon%2i",ipOXYGEN );
00879 conv.BadConvIoniz[0] = dense.xMolecules[ipOXYGEN];
00880 conv.BadConvIoniz[1] = dense.gas_phase[ipOXYGEN];
00881 }
00882 else if( dense.lgElmtOn[ipSILICON] &&
00883 dense.xMolecules[ipSILICON] > dense.gas_phase[ipSILICON]*(1.+EPS_MOLE) )
00884 {
00885 conv.lgConvIoniz = false;
00886 sprintf( conv.chConvIoniz, "COcon%2i",ipSILICON );
00887 conv.BadConvIoniz[0] = dense.xMolecules[ipSILICON];
00888 conv.BadConvIoniz[1] = dense.gas_phase[ipSILICON];
00889 }
00890 else if( dense.lgElmtOn[ipSULPHUR] &&
00891 dense.xMolecules[ipSULPHUR] > dense.gas_phase[ipSULPHUR]*(1.+EPS_MOLE) )
00892 {
00893 conv.lgConvIoniz = false;
00894 sprintf( conv.chConvIoniz, "COcon%2i",ipSULPHUR );
00895 conv.BadConvIoniz[0] = dense.xMolecules[ipSULPHUR];
00896 conv.BadConvIoniz[1] = dense.gas_phase[ipSULPHUR];
00897 }
00898 else if( dense.lgElmtOn[ipNITROGEN] &&
00899 dense.xMolecules[ipNITROGEN] > dense.gas_phase[ipNITROGEN]*(1.+EPS_MOLE) )
00900 {
00901 conv.lgConvIoniz = false;
00902 sprintf( conv.chConvIoniz, "COcon%2i",ipNITROGEN );
00903 conv.BadConvIoniz[0] = dense.xMolecules[ipNITROGEN];
00904 conv.BadConvIoniz[1] = dense.gas_phase[ipNITROGEN];
00905 }
00906
00907 else if( dense.lgElmtOn[ipCHLORINE] &&
00908 dense.xMolecules[ipCHLORINE] > dense.gas_phase[ipCHLORINE]*(1.+EPS_MOLE) )
00909 {
00910 conv.lgConvIoniz = false;
00911 sprintf( conv.chConvIoniz, "COcon%2i",ipCHLORINE );
00912 conv.BadConvIoniz[0] = dense.xMolecules[ipCHLORINE];
00913 conv.BadConvIoniz[1] = dense.gas_phase[ipCHLORINE];
00914 }
00915
00916
00917 co.nitro_dissoc_rate = (float)(CO_dissoc_rate("N"));
00918
00919 DEBUG_EXIT( "CO_solve()" );
00920 return;
00921
00922
00923
00924 }
00925
00926 # undef EPS_MOLE
00927