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
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035 void CO_solve(
00036
00037 bool *lgNegPop,
00038
00039 bool *lgZerPop )
00040 {
00041
00042 int32 merror;
00043 long int i, j, k, n,
00044 nelem , ion , ion2;
00045 double
00046 co_denominator;
00047 realnum cartot_mol, oxytot_mol;
00048 struct chem_element_s *element;
00049 double sum;
00050
00051 DEBUG_ENTRY( "CO_solve()" );
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 > dynamics.n_initial_relax+1 && dynamics.lgAdvection && ! dynamics.lgEquilibrium && dynamics.Rate != 0.0 )
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] = (realnum)(mole.c[element->ipMlP][element->ipMl] -
00214 ionbal.RateRecomTot[element->ipCl][0])*dense.lgElmtOn[element->ipCl];
00215
00216 mole.xMoleChTrRate[element->ipCl][0][1] = (realnum)(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 > dynamics.n_initial_relax+1 && 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 getrf_wrapper(mole.num_comole_calc,mole.num_comole_calc,mole.amat[0],mole.num_comole_calc, ipiv,&merror);
00399
00400 if( merror != 0 )
00401 {
00402 fprintf( ioQQQ, " CO_solve getrf_wrapper error\n" );
00403 cdEXIT(EXIT_FAILURE);
00404 }
00405
00406 getrs_wrapper('N',mole.num_comole_calc,1,mole.amat[0],mole.num_comole_calc,ipiv,mole.b,mole.num_comole_calc,&merror);
00407
00408 if( merror != 0 )
00409 {
00410 fprintf( ioQQQ, " CO_solve: dgetrs finds singular or ill-conditioned matrix\n" );
00411 cdEXIT(EXIT_FAILURE);
00412 }
00413
00414
00415 *lgNegPop = false;
00416 *lgZerPop = false;
00417 co_denominator = 0;
00418 # if 0
00419 if(fnzone > 1.02)
00420 {
00421 for( i=0; i < mole.num_comole_calc; i++)
00422 {
00423
00424 for( j=61; j < 62; j++ )
00425 {
00426
00427 printf("ZONE\t%.5e\tSPECIES_1\t%s\tSPECIES_2\t%s\tDEST_RATE\t%.3e\tbvec\t%.3e\n",
00428 fnzone, COmole[i]->label,COmole[j]->label,mole.c[i][j],mole.b[i]);
00429
00430 }
00431 }
00432 }
00433 # endif
00434 for( i=0; i < mole.num_comole_calc; i++ )
00435 {
00436
00437 if( mole.b[i] < 0. )
00438 {
00439
00440
00441
00442
00443
00444
00445
00446
00447
00448
00449
00450
00451
00452
00453
00454
00455
00456
00457
00458
00459
00460
00461
00462
00463
00464
00465 if( dense.lgElmtOn[COmole[i]->nelem_hevmol] )
00466 {
00467 co_denominator = dense.gas_phase[COmole[i]->nelem_hevmol];
00468 }
00469 else
00470 {
00471
00472 co_denominator = 1e10;
00473 }
00474
00475
00476
00477 ASSERT( co_denominator > SMALLFLOAT || !dense.lgElmtOn[COmole[i]->nelem_hevmol] );
00478
00479
00480
00481
00482
00483
00484
00485 if( mole.b[i] / SDIV(co_denominator) > -1e-5 )
00486 {
00487
00488
00489
00490 mole.b[i] *= -1.;
00491 }
00492 else
00493 {
00494
00495 *lgNegPop = true;
00496
00497
00498
00499
00500
00501 if( nzone>co.co_nzone )
00502 {
00503 static long int nzFail=-2;
00504 long int nzFail2 = (long)fnzone*100;
00505
00506 if( !conv.lgSearch )
00507 fprintf(ioQQQ," PROBLEM ");
00508 fprintf(ioQQQ,
00509 " CO_solve neg pop for species %li %s, value is %.2e rel value is %.2e zone %.2f Te %.4e Search?%c\n",
00510 i ,
00511 COmole[i]->label ,
00512
00513 mole.b[i],
00514
00515 mole.b[i] / SDIV(co_denominator) ,
00516 fnzone,
00517 phycon.te,TorF( conv.lgSearch ) );
00518
00519
00520 if( nzFail2 !=nzFail )
00521 {
00522 nzFail = nzFail2;
00523 ConvFail( "pops" , "CO");
00524 }
00525 }
00526 mole.b[i] *= -1;
00527
00528
00529
00530
00531
00532 }
00533 }
00534 else if( mole.b[i] == 0. )
00535 {
00536
00537
00538 *lgZerPop = true;
00539
00540
00541 # if 0
00542 if( CODEBUG>1 )
00543 {
00544 fprintf(ioQQQ," PROBLEM ");
00545 fprintf(ioQQQ,
00546 " CO_solve zero pop for species %li %s, value is %.2e zone %li\n",
00547 i , COmole[i]->label , mole.b[i], nzone);
00548 }
00549 # endif
00550 mole.b[i] = SMALLFLOAT;
00551 }
00552 COmole[i]->hevmol = (realnum)mole.b[i];
00553
00554 ASSERT( !isnan( COmole[i]->hevmol ) );
00555 }
00556
00557
00558
00559
00560
00561
00562
00563
00564 if( *lgNegPop && (nzone>co.co_nzone &&!conv.lgSearch) )
00565 {
00566 conv.lgConvPops = false;
00567 fprintf(ioQQQ," CO network negative population occurred, Te=%.4e, calling ConvFail. ",
00568 phycon.te);
00569 fprintf( ioQQQ, " CO/C=%9.1e", findspecies("CO")->hevmol/SDIV(dense.gas_phase[ipCARBON]) );
00570 fprintf( ioQQQ, "\n" );
00571
00572 *lgNegPop = false;
00573 }
00574
00575
00576 if( 0 && *lgNegPop )
00577 {
00578
00579 for(j=0; j<2; ++j )
00580 {
00581 double sumcar = 0.;
00582 double sumoxy = 0.;
00583 double renorm;
00584 for( i=0; i < mole.num_comole_calc; i++ )
00585 {
00586
00587 sumcar += COmole[i]->hevmol*COmole[i]->nElem[ipCARBON];
00588 sumoxy += COmole[i]->hevmol*COmole[i]->nElem[ipOXYGEN];
00589 }
00590 renorm = (tot_ion[0] + tot_ion[1]) / SDIV( sumcar+sumoxy);
00591 if(j)
00592 fprintf(ioQQQ,"\t%f\n", renorm);
00593 else
00594 fprintf(ioQQQ,"renormco\t%.2f\t%f", fnzone, renorm);
00595 for( i=0; i < mole.num_comole_calc; i++ )
00596 {
00597 COmole[i]->hevmol *= (realnum)renorm;
00598
00599 ASSERT( !isnan( COmole[i]->hevmol ) );
00600 }
00601 }
00602 }
00603
00604 if( merror != 0 )
00605 {
00606 fprintf( ioQQQ, " COMOLE matrix inversion error, MERROR=%5ld zone=%5ld\n",
00607 (long)merror, nzone );
00608 ShowMe();
00609 fprintf( ioQQQ, " Product matrix\n " );
00610
00611 for( i=0; i < MIN2(mole.num_comole_calc,8); i++ )
00612 {
00613 fprintf( ioQQQ, "%4.4s", COmole[i]->label );
00614 }
00615 fprintf( ioQQQ, " \n" );
00616
00617 for( j=0; j < mole.num_comole_calc; j++ )
00618 {
00619 fprintf( ioQQQ, " %4.4s", COmole[j]->label );
00620 fprintf( ioQQQ, " " );
00621
00622 for( i=0; i < MIN2(mole.num_comole_calc,8); i++ )
00623 {
00624 fprintf( ioQQQ, "%9.1e", mole.amat[i][j]*
00625 COmole[i]->hevmol );
00626 }
00627 fprintf( ioQQQ, "\n" );
00628 }
00629
00630 if( mole.num_comole_calc >= 9 )
00631 {
00632 fprintf( ioQQQ, " COMOLE matrix\n " );
00633 for( i=8; i < mole.num_comole_calc; i++ )
00634 {
00635 fprintf( ioQQQ, "%4.4s", COmole[i]->label );
00636 }
00637 fprintf( ioQQQ, " \n" );
00638
00639 for( j=0; j < mole.num_comole_calc; j++ )
00640 {
00641 fprintf( ioQQQ, " %4.4s", COmole[j]->label );
00642 fprintf( ioQQQ, " " );
00643 for( i=8; i < mole.num_comole_calc; i++ )
00644 {
00645 fprintf( ioQQQ, "%9.1e",
00646 mole.amat[i][j]* COmole[i]->hevmol );
00647 }
00648 fprintf( ioQQQ, "\n" );
00649 }
00650 }
00651
00652 fprintf( ioQQQ, " Mole dens:" );
00653 for( j=0; j < mole.num_comole_calc; j++ )
00654 {
00655 fprintf( ioQQQ, " %4.4s:%.2e", COmole[j]->label
00656 , COmole[j]->hevmol );
00657 }
00658 fprintf( ioQQQ, " \n" );
00659
00660 cdEXIT(EXIT_FAILURE);
00661 }
00662
00663 if( trace.lgTrace && trace.lgTr_CO_Mole )
00664 {
00665 fprintf( ioQQQ, " Product matrix\n " );
00666
00667 for( i=0; i < MIN2(mole.num_comole_calc,8); i++ )
00668 {
00669 fprintf( ioQQQ, "%4.4s", COmole[i]->label );
00670 }
00671 fprintf( ioQQQ, " \n" );
00672
00673 for( j=0; j < mole.num_comole_calc; j++ )
00674 {
00675 fprintf( ioQQQ, " %4.4s", COmole[j]->label );
00676 fprintf( ioQQQ, " " );
00677 for( i=0; i < MIN2(mole.num_comole_calc,8); i++ )
00678 {
00679 fprintf( ioQQQ, "%9.1e", mole.amat[i][j]*
00680 COmole[i]->hevmol );
00681 }
00682 fprintf( ioQQQ, "\n" );
00683
00684 }
00685
00686 if( mole.num_comole_calc >= 9 )
00687 {
00688 fprintf( ioQQQ, " COMOLE matrix\n " );
00689 for( i=8; i < mole.num_comole_calc; i++ )
00690 {
00691 fprintf( ioQQQ, "%4.4s", COmole[i]->label );
00692 }
00693 fprintf( ioQQQ, " \n" );
00694
00695 for( j=0; j < mole.num_comole_calc; j++ )
00696 {
00697 fprintf( ioQQQ, " %4.4s", COmole[j]->label );
00698 fprintf( ioQQQ, " " );
00699
00700 for( i=8; i < mole.num_comole_calc; i++ )
00701 {
00702 fprintf( ioQQQ, "%9.1e", mole.amat[i][j]* COmole[i]->hevmol );
00703 }
00704 fprintf( ioQQQ, "\n" );
00705 }
00706 }
00707
00708 fprintf( ioQQQ, " Mole dens:" );
00709 for( j=0; j < mole.num_comole_calc; j++ )
00710 {
00711 fprintf( ioQQQ, " %4.4s:%.2e", COmole[j]->label
00712 , COmole[j]->hevmol );
00713 }
00714 fprintf( ioQQQ, " \n" );
00715 }
00716
00717
00718 co.CODissHeat = (realnum)(CO_findrate("PHOTON,CO=>C,O")*1e-12);
00719
00720 thermal.heating[0][9] = co.CODissHeat;
00721
00722
00723 for( i=0;i<LIMELM; ++i )
00724 {
00725 dense.xMolecules[i] = 0.;
00726 }
00727
00728
00729
00730 for(i=0;i<N_H_MOLEC;i++)
00731 {
00732 dense.xMolecules[ipHYDROGEN] += hmi.Hmolec[i]*hmi.nProton[i];
00733 }
00734 dense.xMolecules[ipHYDROGEN] -= (hmi.Hmolec[ipMH] + hmi.Hmolec[ipMHp]);
00735
00736
00737
00738
00739
00740
00741
00742 dense.H_sum_in_CO = 0;
00743 for(i=0; i<mole.num_comole_calc; ++i)
00744 {
00745 if(COmole[i]->n_nuclei <= 1)
00746 continue;
00747
00748 dense.H_sum_in_CO += COmole[i]->nElem[ipHYDROGEN]*COmole[i]->hevmol;
00749
00750 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
00751 {
00752 if( mole.lgElem_in_chemistry[nelem] )
00753 {
00754 if( COmole[i]->hevmol > BIGFLOAT )
00755 {
00756 fprintf(ioQQQ, "PROBLEM DISASTER mole_co_solve has found "
00757 "a CO-network molecule with an insane abundance.\n");
00758 fprintf(ioQQQ, "Species number %li with abundance %.2e\n",
00759 i , COmole[i]->hevmol);
00760 TotalInsanity();
00761 }
00762 else if( conv.lgSearch &&
00763 COmole[i]->nElem[nelem]*COmole[i]->hevmol > dense.gas_phase[nelem])
00764 {
00765 if( !conv.lgSearch )
00766 fprintf(ioQQQ,
00767 "PROBLEM COmole limit trip call %li Seatch? %c nMole %li "
00768 "n(mole) %.2e n(gas) %.2e)\n",
00769 conv.nTotalIoniz,
00770 TorF(conv.lgSearch),i ,
00771 COmole[i]->nElem[nelem]*COmole[i]->hevmol,
00772 dense.gas_phase[nelem]);
00773 COmole[i]->hevmol = dense.gas_phase[nelem]/COmole[i]->nElem[nelem]/2.f;
00774 }
00775 dense.xMolecules[nelem] +=COmole[i]->nElem[nelem]*COmole[i]->hevmol;
00776 }
00777 }
00778 }
00779
00780
00781
00782
00783
00784
00785
00786
00787
00788
00789
00790
00791
00792 for(i=0;i<mole.num_elements; ++i)
00793 {
00794 element = chem_element[i];
00795
00796 COmole[element->ipMlP]->hevmol = ((dense.xIonDense[element->ipCl][1]+ COmole[element->ipMlP]->hevmol)/2)*dense.lgElmtOn[element->ipCl];
00797
00798 COmole[element->ipMl]->hevmol = ((dense.xIonDense[element->ipCl][0]+ COmole[element->ipMl]->hevmol)/2)*dense.lgElmtOn[element->ipCl];
00799
00800
00801 ASSERT( !isnan( COmole[element->ipMlP]->hevmol ) );
00802 ASSERT( !isnan( COmole[element->ipMl]->hevmol ) );
00803 }
00804
00805 findspecies("^13CO")->hevmol = findspecies("CO")->hevmol / co.C12_C13_isotope_ratio;
00806
00807
00808 if( dense.lgElmtOn[ipCARBON] &&
00809 fabs(dense.xIonDense[ipCARBON][0]- findspecies("C")->hevmol)/SDIV(dense.gas_phase[ipCARBON]) >0.001 )
00810 {
00811 conv.lgConvIoniz = false;
00812 sprintf( conv.chConvIoniz, "CO C0 con");
00813 conv.BadConvIoniz[0] = dense.xIonDense[ipCARBON][0];
00814 conv.BadConvIoniz[1] = findspecies("C")->hevmol;
00815 }
00816 else if( dense.lgElmtOn[ipCARBON] &&
00817 fabs(dense.xIonDense[ipCARBON][1]- findspecies("C+")->hevmol)/SDIV(dense.gas_phase[ipCARBON]) >0.001 )
00818 {
00819 conv.lgConvIoniz = false;
00820 sprintf( conv.chConvIoniz, "CO C1 con");
00821 conv.BadConvIoniz[0] = dense.xIonDense[ipCARBON][1];
00822 conv.BadConvIoniz[1] = findspecies("C+")->hevmol;
00823 }
00824 else if( dense.lgElmtOn[ipOXYGEN] &&
00825 fabs(dense.xIonDense[ipOXYGEN][0]- findspecies("O")->hevmol)/SDIV(dense.gas_phase[ipOXYGEN]) >0.001 )
00826 {
00827 conv.lgConvIoniz = false;
00828 sprintf( conv.chConvIoniz, "CO O0 con");
00829 conv.BadConvIoniz[0] = dense.xIonDense[ipOXYGEN][0];
00830 conv.BadConvIoniz[1] = findspecies("O")->hevmol;
00831 }
00832 else if( dense.lgElmtOn[ipOXYGEN] &&
00833 fabs(dense.xIonDense[ipOXYGEN][1]- findspecies("O+")->hevmol)/SDIV(dense.gas_phase[ipOXYGEN]) >0.001 )
00834 {
00835 conv.lgConvIoniz = false;
00836 sprintf( conv.chConvIoniz, "CO O1 con");
00837 conv.BadConvIoniz[0] = dense.xIonDense[ipOXYGEN][1];
00838 conv.BadConvIoniz[1] = findspecies("O+")->hevmol;
00839 }
00840
00841
00842 for(i=0;i<mole.num_elements; ++i)
00843 {
00844 element = chem_element[i];
00845 dense.xIonDense[element->ipCl][0] = COmole[element->ipMl]->hevmol*dense.lgElmtOn[element->ipCl];
00846 dense.xIonDense[element->ipCl][1] = COmole[element->ipMlP]->hevmol*dense.lgElmtOn[element->ipCl];
00847 dense.IonLow[element->ipCl] = 0;
00848 dense.IonHigh[element->ipCl] = MAX2( dense.IonHigh[element->ipCl] , 1 );
00849 }
00850
00851
00852 # define EPS_MOLE 0.1
00853 if( dense.lgElmtOn[ipHYDROGEN] &&
00854 dense.xMolecules[ipHYDROGEN] > dense.gas_phase[ipHYDROGEN]*(1.+EPS_MOLE) )
00855 {
00856 conv.lgConvIoniz = false;
00857 sprintf( conv.chConvIoniz, "COcon%2i",ipHYDROGEN );
00858 conv.BadConvIoniz[0] = dense.xMolecules[ipHYDROGEN];
00859 conv.BadConvIoniz[1] = dense.gas_phase[ipHYDROGEN];
00860 }
00861 else if( dense.lgElmtOn[ipCARBON] &&
00862 dense.xMolecules[ipCARBON] > dense.gas_phase[ipCARBON]*(1.+EPS_MOLE) )
00863 {
00864 conv.lgConvIoniz = false;
00865 sprintf( conv.chConvIoniz, "COcon%2i",ipCARBON );
00866 conv.BadConvIoniz[0] = dense.xMolecules[ipCARBON];
00867 conv.BadConvIoniz[1] = dense.gas_phase[ipCARBON];
00868 }
00869 else if( dense.lgElmtOn[ipOXYGEN] &&
00870 dense.xMolecules[ipOXYGEN] > dense.gas_phase[ipOXYGEN]*(1.+EPS_MOLE) )
00871 {
00872 conv.lgConvIoniz = false;
00873 sprintf( conv.chConvIoniz, "COcon%2i",ipOXYGEN );
00874 conv.BadConvIoniz[0] = dense.xMolecules[ipOXYGEN];
00875 conv.BadConvIoniz[1] = dense.gas_phase[ipOXYGEN];
00876 }
00877 else if( dense.lgElmtOn[ipSILICON] &&
00878 dense.xMolecules[ipSILICON] > dense.gas_phase[ipSILICON]*(1.+EPS_MOLE) )
00879 {
00880 conv.lgConvIoniz = false;
00881 sprintf( conv.chConvIoniz, "COcon%2i",ipSILICON );
00882 conv.BadConvIoniz[0] = dense.xMolecules[ipSILICON];
00883 conv.BadConvIoniz[1] = dense.gas_phase[ipSILICON];
00884 }
00885 else if( dense.lgElmtOn[ipSULPHUR] &&
00886 dense.xMolecules[ipSULPHUR] > dense.gas_phase[ipSULPHUR]*(1.+EPS_MOLE) )
00887 {
00888 conv.lgConvIoniz = false;
00889 sprintf( conv.chConvIoniz, "COcon%2i",ipSULPHUR );
00890 conv.BadConvIoniz[0] = dense.xMolecules[ipSULPHUR];
00891 conv.BadConvIoniz[1] = dense.gas_phase[ipSULPHUR];
00892 }
00893 else if( dense.lgElmtOn[ipNITROGEN] &&
00894 dense.xMolecules[ipNITROGEN] > dense.gas_phase[ipNITROGEN]*(1.+EPS_MOLE) )
00895 {
00896 conv.lgConvIoniz = false;
00897 sprintf( conv.chConvIoniz, "COcon%2i",ipNITROGEN );
00898 conv.BadConvIoniz[0] = dense.xMolecules[ipNITROGEN];
00899 conv.BadConvIoniz[1] = dense.gas_phase[ipNITROGEN];
00900 }
00901
00902 else if( dense.lgElmtOn[ipCHLORINE] &&
00903 dense.xMolecules[ipCHLORINE] > dense.gas_phase[ipCHLORINE]*(1.+EPS_MOLE) )
00904 {
00905 conv.lgConvIoniz = false;
00906 sprintf( conv.chConvIoniz, "COcon%2i",ipCHLORINE );
00907 conv.BadConvIoniz[0] = dense.xMolecules[ipCHLORINE];
00908 conv.BadConvIoniz[1] = dense.gas_phase[ipCHLORINE];
00909 }
00910
00911
00912 co.nitro_dissoc_rate = (realnum)(CO_dissoc_rate("N"));
00913 return;
00914
00915
00916
00917 }
00918
00919 # undef EPS_MOLE