00001
00002
00003
00004 #include "cddefines.h"
00005 #include "trace.h"
00006 #include "stopcalc.h"
00007 #include "struc.h"
00008 #include "rfield.h"
00009 #include "mole.h"
00010 #include "dense.h"
00011 #include "heavy.h"
00012 #include "wind.h"
00013 #include "geometry.h"
00014 #include "thermal.h"
00015 #include "radius.h"
00016 #include "phycon.h"
00017 #include "pressure.h"
00018 #include "ionbal.h"
00019 #include "conv.h"
00020
00021 #define XMULTP 3.0
00022
00023 int ConvInitSolution(void)
00024 {
00025 long int i,
00026 ionlup,
00027 isear,
00028 j,
00029 looplimit,
00030 nelem ,
00031 nflux_old,
00032 nelem_reflux ,
00033 ion_reflux;
00034
00035 double cnew,
00036 cold,
00037 dt,
00038 fac2,
00039 factor,
00040 hnew,
00041 hold,
00042 Cool_min_Heat,
00043 t1,
00044 t2,
00045 tinc,
00046 tmax,
00047 tmin,
00048 tnew,
00049 told;
00050
00051 double derrdt;
00052
00053 DEBUG_ENTRY( "ConvInitSolution()" );
00054
00055 tnew = DBL_MAX;
00056 cnew = DBL_MAX;
00057 hnew = DBL_MAX;
00058
00059
00060 conv.nPres2Ioniz = 0;
00061
00062
00063 conv.nTotalIoniz = 0;
00064
00065 conv.lgOscilOTS = false;
00066
00067 lgAbort = false;
00068 dense.lgEdenBad = false;
00069 dense.nzEdenBad = 0;
00070
00071
00072 conv.BigEdenError = 0.;
00073 conv.AverEdenError = 0.;
00074 conv.BigHeatCoolError = 0.;
00075 conv.AverHeatCoolError = 0.;
00076 conv.BigPressError = 0.;
00077 conv.AverPressError = 0.;
00078
00079
00080
00081 if( radius.sdrmin == radius.sdrmax )
00082 {
00083 radius.drad = MIN2( radius.sdrmax , radius.drad );
00084 radius.drad_x_fillfac = radius.drad * geometry.FillFac;
00085 }
00086
00087
00088
00089
00090
00091
00092 ASSERT( dense.xIonDense[ipHYDROGEN][0] >0 && dense.xIonDense[ipHYDROGEN][1]>= 0.);
00093
00094 if( trace.nTrConvg )
00095 {
00096 fprintf( ioQQQ, " \n" );
00097 fprintf( ioQQQ, " ConvInitSolution entered \n" );
00098 }
00099
00100
00101
00102
00103
00104
00105 if( thermal.ConstTemp != 0. )
00106 {
00107
00108 if( trace.lgTrace )
00109 {
00110 fprintf( ioQQQ, " ConvInitSolution called, forcing TE to%10.2e\n",
00111 thermal.ConstTemp );
00112 }
00113 phycon.te = thermal.ConstTemp;
00114 tfidle(false);
00115
00116
00117
00118 PresTotCurrent();
00119
00120
00121 conv.lgSearch = true;
00122
00123 ionlup = 0;
00124 conv.lgConvIoniz = false;
00125
00126
00127
00128
00129
00130
00131
00132 while( (ionlup<3) || (!lgAbort && (ionlup <= 3) && !conv.lgConvIoniz && (ionlup <= 7) ) )
00133 {
00134
00135
00136 PresTotCurrent();
00137 if( ConvEdenIoniz() )
00138 {
00139 DEBUG_EXIT( "ConvInitSolution()" );
00140
00141 return(1);
00142 }
00143 if( ionlup<3 )
00144 {
00145
00146
00147
00148 radius_first();
00149 }
00150
00151 ++ionlup;
00152 }
00153
00154 if( trace.lgTrace || trace.nTrConvg )
00155 {
00156 fprintf( ioQQQ, " ========================================================================================================================\n");
00157 fprintf( ioQQQ, " ConvInitSolution: search set false 1 Te=%.3e================================================================================\n" , phycon.te);
00158 fprintf( ioQQQ, " ========================================================================================================================\n");
00159 }
00160 conv.lgSearch = false;
00161 }
00162
00163
00164
00165
00166
00167
00168 else if( iteration != 1 )
00169 {
00170
00171 if( trace.lgTrace || trace.nTrConvg )
00172 {
00173 fprintf( ioQQQ, " ConvInitSolution called, ITER=%2ld resetting Te to %10.2e\n",
00174 iteration, struc.testr[0] );
00175 }
00176
00177 if( trace.lgTrace || trace.nTrConvg )
00178 {
00179 fprintf( ioQQQ, " search set true\n" );
00180 }
00181
00182
00183
00184
00185 conv.lgSearch = true;
00186
00187
00188
00189
00190 phycon.te = struc.testr[0];
00191 tfidle(false);
00192
00193
00194
00195 PresTotCurrent();
00196
00197
00198
00199
00200 if( ConvPresTempEdenIoniz() )
00201 {
00202 DEBUG_EXIT( "ConvInitSolution()" );
00203
00204 return(1);
00205 }
00206
00207 if( trace.lgTrace || trace.nTrConvg )
00208 {
00209 fprintf( ioQQQ, " ========================================================================================================================\n");
00210 fprintf( ioQQQ, " ConvInitSolution: search set false 2 Te=%.3e========================================================================================\n" , phycon.te);
00211 fprintf( ioQQQ, " ========================================================================================================================\n");
00212 }
00213 conv.lgSearch = false;
00214
00215
00216 if( ConvTempEdenIoniz() )
00217 {
00218 DEBUG_EXIT( "ConvInitSolution()" );
00219
00220 return(1);
00221 }
00222
00223
00224
00225 PresTotCurrent();
00226
00227 }
00228
00229 else
00230 {
00231
00232
00233
00234
00235
00236
00237
00238 conv.lgSearch = true;
00239
00240 if( trace.lgTrace )
00241 {
00242 fprintf( ioQQQ, " ConvInitSolution called, new temperature.\n" );
00243 }
00244
00245
00246 if( thermal.lgTeHigh )
00247 {
00248
00249 phycon.te = 1e6;
00250 tfidle(false);
00251 ionbal.lgNoDec = false;
00252 factor = 1./XMULTP;
00253 }
00254 else
00255 {
00256
00257 ionbal.lgNoDec = true;
00258 phycon.te = 4000.;
00259 tfidle(false);
00260 factor = XMULTP;
00261 }
00262
00263
00264 thermal.lgColNeg = false;
00265
00266
00267 PresTotCurrent();
00268
00269 thermal.htot = 1.;
00270 thermal.ctot = 1.;
00271
00272
00273
00274 if( trace.nTrConvg )
00275 {
00276 fprintf( ioQQQ, " ConvInitSolution calling ConvEdenIoniz1 with te=%10.3e\n",
00277 phycon.te );
00278 }
00279
00280 for( ionlup=0; ionlup<2; ++ionlup )
00281 {
00282
00283 if( ConvEdenIoniz() )
00284 {
00285 DEBUG_EXIT( "ConvInitSolution()" );
00286
00287 return(1);
00288 }
00289
00290
00291
00292 radius_first();
00293 }
00294
00295 if( lgAbort )
00296 {
00297
00298 fprintf( ioQQQ, " Search for initial conditions aborted - lgAbort set true.\n" );
00299 ShowMe();
00300
00301 DEBUG_EXIT( "ConvInitSolution()" );
00302
00303 return(1);
00304 }
00305
00306
00307 told = phycon.te;
00308 cold = thermal.ctot;
00309 hold = thermal.htot;
00310
00311 if( trace.lgTrace )
00312 {
00313 fprintf( ioQQQ, " return to ConvInitSolution 1, HTOT:%10.3e CTOT:%10.3e Te:%11.4e EDEN%11.4e<<<<<<<<<<<<<<<<<<<<<<<<<\n",
00314 thermal.htot, thermal.ctot, phycon.te, dense.eden );
00315 }
00316
00317
00318 Cool_min_Heat = fabs(thermal.ctot-thermal.htot)/(thermal.ctot - thermal.htot);
00319
00320
00321 if( !thermal.lgTeHigh && Cool_min_Heat >= 0. )
00322 {
00323
00324
00325 if( trace.nTrConvg )
00326 {
00327 fprintf( ioQQQ, " ConvInitSolution entering T<ConvInitSolution loop, te htot ctot=%10.3e%10.3e%10.3e Cool_min_Heat=%9.1e\n",
00328 phycon.te, thermal.htot, thermal.ctot, Cool_min_Heat );
00329 }
00330
00331
00332
00333
00334
00335
00336 fac2 = 0.8;
00337 while( phycon.te*fac2 > StopCalc.TeLowest && Cool_min_Heat >= 0. && !lgAbort )
00338 {
00339 double oxy_in_grains;
00340 long int ion;
00341
00342
00343
00344 told = phycon.te;
00345 cold = thermal.ctot;
00346 hold = thermal.htot;
00347
00348
00349
00350
00351
00352 float frac_mole_max = 0.;
00353 for( i=0; i<mole.num_comole_calc; ++i )
00354 {
00355 if( dense.lgElmtOn[COmole[i]->nelem_hevmol] )
00356 {
00357 float frac = COmole[i]->hevmol / SDIV(dense.gas_phase[COmole[i]->nelem_hevmol]);
00358 frac *= (float) COmole[i]->nElem[COmole[i]->nelem_hevmol];
00359 frac_mole_max = MAX2( frac_mole_max , frac );
00360 }
00361 }
00362
00363
00364
00365 oxy_in_grains = 0;
00366 for(i=0;i<mole.num_comole_calc;++i)
00367 {
00368 oxy_in_grains += (1 - COmole[i]->lgGas_Phase)*COmole[i]->hevmol*COmole[i]->nElem[ipOXYGEN];
00369 }
00370
00371 oxy_in_grains /= SDIV(dense.gas_phase[ipOXYGEN]);
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381 if( oxy_in_grains > 0.99 )
00382 {
00383 fac2 = 0.999;
00384 }
00385 else if( oxy_in_grains > 0.9 )
00386 {
00387 fac2 = 0.99;
00388 }
00389
00390 else if( phycon.te < 5. ||
00391
00392 oxy_in_grains > 0.1 )
00393 {
00394 fac2 = 0.97;
00395 }
00396
00397 else if( phycon.te < 20. || frac_mole_max > 0.1 )
00398 {
00399
00400
00401 fac2 = 0.95;
00402 }
00403 else
00404 {
00405
00406 fac2 = 0.8;
00407 }
00408
00409
00410 phycon.te = phycon.te * fac2;
00411
00412 tfidle(true);
00413
00414
00415 PresTotCurrent();
00416
00417
00418 if( trace.nTrConvg )
00419 {
00420 fprintf( ioQQQ, "\n ConvInitSolution calling ConvEdenIoniz2, te:%.2e fac2:%.3e\n",
00421 phycon.te , fac2);
00422 }
00423
00424 if( ConvEdenIoniz() )
00425 {
00426 DEBUG_EXIT( "ConvInitSolution()" );
00427
00428 return(1);
00429 }
00430
00431 tnew = phycon.te;
00432 cnew = thermal.ctot;
00433 hnew = thermal.htot;
00434
00435 if( trace.lgTrace || trace.nTrConvg )
00436 {
00437 fprintf( ioQQQ, " return to ConvInitSolution 2 Te:%.4e Htot:%.3e Ctot:%.3e eden:%.4e\n",
00438 phycon.te, thermal.htot, thermal.ctot, dense.eden );
00439 {
00440
00441 enum {DEBUG_LOC=true};
00442
00443 if( DEBUG_LOC )
00444 {
00445 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
00446 {
00447 if( dense.lgElmtOn[nelem] )
00448 {
00449 fprintf(ioQQQ,"DEBUG nelem=%li" , nelem);
00450 for( ion=0; ion<=dense.IonHigh[nelem]; ++ion )
00451 {
00452 fprintf(ioQQQ,"\t%.2e", dense.xIonDense[nelem][ion]/SDIV(dense.gas_phase[nelem]) );
00453 }
00454 fprintf(ioQQQ,"\n");
00455 }
00456 }
00457 fflush(ioQQQ);
00458 }
00459 }
00460 {
00461 FILE *test = stderr;
00462 if( ioQQQ == test )
00463 fflush( ioQQQ );
00464 }
00465 }
00466
00467 Cool_min_Heat = fabs(thermal.ctot-thermal.htot)/(thermal.ctot - thermal.htot);
00468 }
00469
00470 if( lgAbort )
00471 {
00472
00473 fprintf( ioQQQ, " Search for initial conditions aborted - lgAbort set true.\n" );
00474 ShowMe();
00475
00476 DEBUG_EXIT( "ConvInitSolution()" );
00477
00478 return(1);
00479 }
00480 else if( thermal.ctot > thermal.htot )
00481 {
00482
00483 fprintf( ioQQQ, " Equilibrium temperature below T=%8.2e. If this is OK, then reset with LOWEST command. This is ConvInitSolution. TeLowest was%9.2e\n",
00484 phycon.te, StopCalc.TeLowest );
00485
00486 DEBUG_EXIT( "ConvInitSolution()" );
00487
00488 return(1);
00489 }
00490
00491
00492 phycon.te = phycon.te / fac2;
00493 tfidle(true);
00494 factor = 1./(0.5*(1. + fac2));
00495 Cool_min_Heat = fabs(cold-hold)/(cold - hold);
00496 if( trace.nTrConvg )
00497 {
00498 fprintf( ioQQQ, " <4000K loop returns, te=%9.2e Cool_min_Heat%9.2e factor%9.2e\n",
00499 phycon.te, Cool_min_Heat, factor );
00500 }
00501 }
00502
00503 tmin = 1e9;
00504 tmax = 0.;
00505
00506
00507 t1 = StopCalc.TeLowest/1.001;
00508
00509
00510 t2 = StopCalc.TeHighest*1.001/factor;
00511 if( trace.nTrConvg )
00512 {
00513 fprintf( ioQQQ, " ConvInitSolution entering second loop te=%9.2e Cool_min_Heat%9.2e factor%9.2e t1, t2=%10.2e%10.2e\n",
00514 phycon.te, Cool_min_Heat, factor, t1, t2 );
00515 }
00516
00517 while( (thermal.ctot - thermal.htot)*Cool_min_Heat > 0. && phycon.te > t1 && phycon.te < t2 )
00518 {
00519
00520 told = phycon.te;
00521 cold = thermal.ctot;
00522 hold = thermal.htot;
00523 phycon.te = phycon.te * factor;
00524 tfidle(true);
00525
00526
00527
00528
00529
00530 if( factor > 1. )
00531 {
00532 for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
00533 {
00534 if( dense.lgElmtOn[nelem] )
00535 {
00536 dense.IonHigh[nelem] = nelem + 1;
00537 }
00538 }
00539 }
00540
00541
00542 PresTotCurrent();
00543
00544
00545
00546
00547
00548
00549
00550 tmin = MIN2(tmin,phycon.te);
00551 tmax = MAX2(tmax,phycon.te);
00552 conv.lgConvIoniz = false;
00553 strcpy( conv.chConvIoniz, "ConvInitSolution srh" );
00554 isear = 1;
00555 while( isear <= 6 && !conv.lgConvIoniz && !lgAbort )
00556 {
00557
00558 if( trace.nTrConvg )
00559 {
00560 fprintf( ioQQQ, "\n calling ConvEdenIoniz3, te=%9.2e\n",
00561 phycon.te );
00562 }
00563
00564
00565
00566 PresTotCurrent();
00567 if( ConvEdenIoniz() )
00568 {
00569 DEBUG_EXIT( "ConvInitSolution()" );
00570
00571 return(1);
00572 }
00573 if( trace.nTrConvg )
00574 {
00575 fprintf( ioQQQ, " ConvEdenIoniz3 returns, te=%9.2e htot%9.2e ctot%9.2e\n",
00576 phycon.te, thermal.htot, thermal.ctot );
00577 }
00578 if( trace.lgTrace )
00579 {
00580 fprintf( ioQQQ, " return to ConvInitSolution 3 HTOT:%10.3e ctot:%10.3e Te:%11.4e eden%11.4e<<<<<<<<<<<<<<<<<<<<<<<<<\n",
00581 thermal.htot, thermal.ctot, phycon.te, dense.eden );
00582 }
00583 isear += 1;
00584 }
00585 tnew = phycon.te;
00586 cnew = thermal.ctot;
00587 hnew = thermal.htot;
00588 }
00589
00590
00591 if( phycon.te <= t1 || phycon.te >= t2 )
00592 {
00593 if( phycon.te <= t1 )
00594 {
00595 fprintf( ioQQQ, " ConvInitSolution finds TE below TE = %.3eK, HTOT=%12.4e CTOT=%12.4e\n",
00596 t1, thermal.htot, thermal.ctot );
00597 }
00598 else if( phycon.te >= t2 )
00599 {
00600 fprintf( ioQQQ, " ConvInitSolution finds TE above TE = %.3eK, HTOT=%12.4e last CTOT=%12.4e\n",
00601 t2, thermal.htot, thermal.ctot );
00602 fprintf( ioQQQ,
00603 " The electron temperature is too high; the electrons are relativistic and little optical/UV radiation will result.\n");
00604 fprintf( ioQQQ,
00605 " This is often due to the continuum shape being wrong - unphysically strong radio or gamma ray emission will do this.\n");
00606 }
00607 else
00608 {
00609 fprintf( ioQQQ, " ConvInitSolution finds insanity \n" );
00610 ShowMe();
00611 }
00612
00613 DEBUG_EXIT( "ConvInitSolution()" );
00614
00615 return(1);
00616 }
00617
00618
00619
00620 ionbal.lgNoDec = false;
00621 thermal.lgColNeg = true;
00622
00623
00624 tinc = 1./pow(factor,0.2);
00625
00626 for( j=0; j < 2; j++ )
00627 {
00628 if( (thermal.ctot - thermal.htot > 1e-30) && thermal.ctot > 0. && !lgAbort )
00629 {
00630 Cool_min_Heat = fabs(thermal.ctot-thermal.htot)/(thermal.ctot - thermal.htot);
00631 if( trace.nTrConvg )
00632 {
00633 fprintf( ioQQQ, "\n\n >>>entering second loop, te reset to%9.2e\n",
00634 phycon.te );
00635 }
00636
00637
00638 PresTotCurrent();
00639
00640
00641 isear = 1;
00642 conv.lgConvIoniz = false;
00643
00644 if( trace.lgTrace )
00645 {
00646 fprintf( ioQQQ, " return to ConvInitSolution 4 HTOT:%10.3e CTOT:%10.3e Te:%11.4e EDEN%11.4e<<<<<<<<<<<<<<<<<<<<<<<<<\n",
00647 thermal.htot, thermal.ctot, phycon.te, dense.eden );
00648 }
00649 i = 0;
00650
00651
00652 looplimit = 12;
00653 told = phycon.te;
00654 cold = thermal.ctot;
00655 hold = thermal.htot;
00656
00657
00658 while( i < looplimit && (thermal.ctot - thermal.htot)*Cool_min_Heat > 0. &&
00659 !lgAbort )
00660 {
00661 i += 1;
00662 told = phycon.te;
00663 cold = thermal.ctot;
00664 hold = thermal.htot;
00665 phycon.te = phycon.te *tinc;
00666 tfidle(true);
00667
00668
00669 PresTotCurrent();
00670
00671
00672 if( ConvEdenIoniz() )
00673 {
00674 DEBUG_EXIT( "ConvInitSolution()" );
00675
00676 return(1);
00677 }
00678
00679 if( trace.nTrConvg )
00680 {
00681 fprintf( ioQQQ, " ConvEdenIoniz returned, te=%9.2e htot%9.2e ctot%9.2e\n",
00682 phycon.te, thermal.htot, thermal.ctot );
00683 }
00684
00685 tnew = phycon.te;
00686 cnew = thermal.ctot;
00687 hnew = thermal.htot;
00688
00689 if( trace.lgTrace )
00690 {
00691 fprintf( ioQQQ, " return to ConvInitSolution 4 HTOT:%10.3e CTOT:%10.3e Te:%11.4e EDEN%11.4e<<<<<<<<<<<<<<<<<<<<<<<<<\n",
00692 thermal.htot, thermal.ctot, phycon.te, dense.eden );
00693 }
00694 }
00695
00696
00697 tinc = 1./pow(tinc,0.2);
00698 }
00699 }
00700
00701 if( trace.lgTrace || trace.nTrConvg )
00702 {
00703 fprintf( ioQQQ, " ========================================================================================================================\n");
00704 fprintf( ioQQQ, " ConvInitSolution: search set false 3 Te=%.3e========================================================================================\n" , phycon.te);
00705 fprintf( ioQQQ, " ========================================================================================================================\n");
00706 }
00707 conv.lgSearch = false;
00708
00709
00710
00711
00712 if( tnew != told && tnew > 0. && told > 0. )
00713
00714 {
00715 derrdt = ((cnew-hnew)-(cold-hold))/(tnew-told);
00716 dt = -(cold-hold)/derrdt;
00717 phycon.te = (told + dt);
00718 phycon.te = MAX2(phycon.te,MIN2(told,tnew));
00719 phycon.te = MIN2(phycon.te,MAX2(told,tnew));
00720
00721 struc.DenMass[0] = dense.xMassDensity;
00722
00723
00724 PresTotCurrent();
00725
00726 }
00727
00728 if( trace.lgTrace )
00729 {
00730 fprintf( ioQQQ, " ConvInitSolution calls IONTE after interpolating a temperature of%10.3e Told=%10.2e Tnew=%10.2e dC-Hold%10.2e dC-Hnew%10.2e\n",
00731 phycon.te, told, tnew, cold - hold, cnew - hnew );
00732 }
00733
00734 if( trace.nTrConvg )
00735 {
00736 fprintf( ioQQQ, "\n\n ConvInitSolution calling ConvTempEdenIoniz with interpolated te=%10.3e\n",
00737 phycon.te );
00738 }
00739
00740 if( ConvTempEdenIoniz() )
00741 {
00742 DEBUG_EXIT( "ConvInitSolution()" );
00743
00744 return(1);
00745 }
00746
00747 if( trace.lgTrace )
00748 {
00749 fprintf( ioQQQ, " ConvInitSolution return, TE:%10.2e==================\n",
00750 phycon.te );
00751 }
00752 }
00753
00754
00755
00756 conv.nPres2Ioniz = 0;
00757
00758 dense.lgEdenBad = false;
00759 dense.nzEdenBad = 0;
00760
00761
00762
00763 conv.nTotalIoniz_start = conv.nTotalIoniz;
00764
00765
00766
00767 conv.BigEdenError = 0.;
00768 conv.AverEdenError = 0.;
00769 conv.BigHeatCoolError = 0.;
00770 conv.AverHeatCoolError = 0.;
00771 conv.BigPressError = 0.;
00772 conv.AverPressError = 0.;
00773
00774
00775 if( iteration == 1 )
00776 {
00777
00778
00779 struc.testr[0] = (float)phycon.te;
00780
00781 struc.hden[0] = dense.gas_phase[ipHYDROGEN];
00782
00783 struc.DenParticles[0] = dense.pden;
00784 struc.heatstr[0] = thermal.htot;
00785 struc.coolstr[0] = thermal.ctot;
00786 struc.volstr[0] = (float)radius.dVeff;
00787 struc.drad_x_fillfac[0] = (float)radius.drad_x_fillfac;
00788 struc.histr[0] = dense.xIonDense[ipHYDROGEN][0];
00789 struc.hiistr[0] = dense.xIonDense[ipHYDROGEN][1];
00790 struc.ednstr[0] = (float)dense.eden;
00791 struc.o3str[0] = dense.xIonDense[ipOXYGEN][2];
00792 struc.DenMass[0] = dense.xMassDensity;
00793 struc.drad[0] = (float)radius.drad;
00794 }
00795
00796
00797
00798
00799
00800
00801
00802
00803 nflux_old = rfield.nflux;
00804 nelem_reflux = -1;
00805 ion_reflux = -1;
00806 for( nelem=2; nelem < LIMELM; nelem++ )
00807 {
00808
00809 for( i=0; i < nelem; i++ )
00810 {
00811 if( dense.xIonDense[nelem][i+1] > 0. )
00812 {
00813 if( Heavy.ipHeavy[nelem][i] > rfield.nflux )
00814 {
00815 rfield.nflux = Heavy.ipHeavy[nelem][i];
00816 nelem_reflux = nelem;
00817 ion_reflux = i;
00818 }
00819 }
00820 }
00821 }
00822
00823
00824 if( nflux_old != rfield.nflux )
00825 {
00826
00827 rfield_opac_zero( nflux_old-1 , rfield.nflux );
00828
00829
00830
00831
00832
00833 PresTotCurrent();
00834
00835
00836 if( ConvBase(1) )
00837 {
00838
00839 DEBUG_EXIT( "ConvInitSolution()" );
00840 return( 1 );
00841 }
00842
00843
00844 if( trace.lgTrace )
00845 {
00846 fprintf(ioQQQ," nflux updated from %li to %li, anu from %g to %g \n",
00847 nflux_old , rfield.nflux ,
00848 rfield.anu[nflux_old] , rfield.anu[rfield.nflux] );
00849 fprintf(ioQQQ," caused by element %li ion %li \n",
00850 nelem_reflux ,ion_reflux );
00851 }
00852 }
00853
00854
00855
00856 if( (strcmp(dense.chDenseLaw,"WIND") == 0) && wind.windv < 0. )
00857 {
00858
00859
00860
00861 # define PCHNG 0.98
00862
00863 if( ConvPresTempEdenIoniz() )
00864 {
00865 DEBUG_EXIT( "ConvInitSolution()" );
00866
00867 return(1);
00868 }
00869
00870 pressure.PresTotlInit *= PCHNG;
00871 if( ConvPresTempEdenIoniz() )
00872 {
00873 DEBUG_EXIT( "ConvInitSolution()" );
00874
00875 return(1);
00876 }
00877
00878 pressure.PresTotlInit /= PCHNG;
00879 if( ConvPresTempEdenIoniz() )
00880 {
00881 DEBUG_EXIT( "ConvInitSolution()" );
00882
00883 return(1);
00884 }
00885
00886 # undef PCHNG
00887 }
00888
00889 DEBUG_EXIT( "ConvInitSolution()" );
00890
00891
00892 return( lgAbort );
00893 }
00894