27 bool lgHi_Down =
false;
32 bool lgLo_Down =
false;
33 long int ion_lo_save =
dense.IonLow[nelem],
34 ion_hi_save =
dense.IonHigh[nelem];
82 trimhi =
MIN2( trimhi , 1e-20f );
94 trimhi =
MIN2( trimhi , 1e-17f );
99 trimhi =
MIN2( trimhi , 1e-12f );
109 if(
dense.IonHigh[nelem] ==2)
111 trimhi =
MIN2(trimhi, 1e-20f);
125 (
dense.xIonDense[nelem][
dense.IonHigh[nelem]]/
dense.gas_phase[nelem] <
dense.density_low_limit ||
126 dense.xIonDense[nelem][
dense.IonHigh[nelem]] <
dense.density_low_limit ) &&
132 dense.xIonDense[nelem][
dense.IonHigh[nelem]-1] +=
134 dense.xIonDense[nelem][
dense.IonHigh[nelem]] = 0.;
136 if(
dense.IonHigh[nelem] == nelem+1-
NISO )
138 long int ipISO = nelem -
dense.IonHigh[nelem];
140 for(
long level = 0; level <
iso_sp[ipISO][nelem].numLevels_max; level++ )
141 iso_sp[ipISO][nelem].st[level].Pop() = 0.;
145 --
dense.IonHigh[nelem];
153 (
dense.xIonDense[nelem][
dense.IonLow[nelem]]/
dense.gas_phase[nelem] <
dense.density_low_limit ||
154 dense.xIonDense[nelem][
dense.IonLow[nelem]] <
dense.density_low_limit ) &&
155 dense.IonLow[nelem] <
dense.IonHigh[nelem] - 1 )
159 dense.xIonDense[nelem][
dense.IonLow[nelem]+1] +=
161 dense.xIonDense[nelem][
dense.IonLow[nelem]] = 0.;
165 if(
dense.IonLow[nelem] == nelem+1-
NISO )
167 long int ipISO = nelem -
dense.IonLow[nelem];
169 for(
long level = 0; level <
iso_sp[ipISO][nelem].numLevels_max; level++ )
170 iso_sp[ipISO][nelem].st[level].Pop() = 0.;
174 ++
dense.IonLow[nelem];
185 if(
dense.IonHigh[nelem] > 1 &&
186 (
dense.IonLow[nelem]==
dense.IonHigh[nelem]-1) &&
187 (
dense.xIonDense[nelem][
dense.IonLow[nelem]] <
dense.density_low_limit) )
191 "PROBLEM DISASTER the density of ion %li of element %s is too small to be computed on this cpu.\n",
192 dense.IonLow[nelem]+1,
195 "Turn off the element with the command ELEMENT XXX OFF or do not consider "
196 "gas with low density, the current hydrogen density is %.2e cm-3.\n",
199 "The outer radius of the cloud is %.2e cm - should a stopping "
203 "The abort flag is being set - the calculation is stopping.\n");
213 ( (
dense.xIonDense[nelem][
dense.IonHigh[nelem]]/
dense.gas_phase[nelem] <=
215 (
dense.xIonDense[nelem][
dense.IonHigh[nelem]] <=
dense.density_low_limit ) )
220 if(
dense.IonHigh[nelem]>1 ||
221 (
dense.IonHigh[nelem]==1&&
dense.xIonDense[nelem][1]<100.*
dense.density_low_limit) )
223 dense.xIonDense[nelem][
dense.IonHigh[nelem]-1] +=
225 dense.xIonDense[nelem][
dense.IonHigh[nelem]] = 0.;
227 if(
dense.IonHigh[nelem] == nelem+1-
NISO )
229 long int ipISO = nelem -
dense.IonHigh[nelem];
231 for(
long level = 0; level <
iso_sp[ipISO][nelem].numLevels_max; level++ )
232 iso_sp[ipISO][nelem].st[level].Pop() = 0.;
234 --
dense.IonHigh[nelem];
240 lgHi_Up_enable =
true;
243 lgHi_Up_enable =
false;
245 if(
dense.IonHigh[nelem]>=nelem+1 )
246 lgHi_Up_enable =
false;
249 lgHi_Up_enable =
false;
252 lgHi_Up_enable =
false;
265 dense.xIonDense[nelem][
dense.IonHigh[nelem]]/
dense.gas_phase[nelem] > 1e-4f &&
267 abundnew > abundold*1.01 )
271 ++
dense.IonHigh[nelem];
275 dense.xIonDense[nelem][
dense.IonHigh[nelem]-1] -=
277 long int ipISO = nelem -
dense.IonHigh[nelem];
278 if (ipISO >= 0 && ipISO <
NISO)
279 iso_sp[ipISO][nelem].st[0].Pop() =
dense.xIonDense[nelem][
dense.IonHigh[nelem]];
287 if(
dense.xIonDense[nelem][
dense.IonLow[nelem]]/
dense.gas_phase[nelem] > 1e-3f &&
288 dense.IonLow[nelem] > 0 )
291 --
dense.IonLow[nelem];
295 dense.xIonDense[nelem][
dense.IonLow[nelem]+1] -=
297 long int ipISO = nelem -
dense.IonLow[nelem];
298 if (ipISO >= 0 && ipISO <
NISO)
299 iso_sp[ipISO][nelem].st[0].Pop() =
dense.xIonDense[nelem][
dense.IonLow[nelem]];
303 else if(
nzone < 10 &&
305 (
dense.IonLow[nelem] < (
dense.IonHigh[nelem] - 2) ) )
308 dense.xIonDense[nelem][
dense.IonLow[nelem]+1] +=
310 dense.xIonDense[nelem][
dense.IonLow[nelem]] = 0.;
313 if(
dense.IonLow[nelem] == nelem+1-
NISO )
315 long int ipISO = nelem -
dense.IonLow[nelem];
317 for(
long level = 0; level <
iso_sp[ipISO][nelem].numLevels_max; level++ )
318 iso_sp[ipISO][nelem].st[level].Pop() = 0.;
320 ++
dense.IonLow[nelem];
324 else if( (
dense.xIonDense[nelem][
dense.IonLow[nelem]] <
dense.density_low_limit) &&
328 (
dense.IonLow[nelem] <
dense.IonHigh[nelem]-1) )
330 while(
dense.xIonDense[nelem][
dense.IonLow[nelem]] <
dense.density_low_limit &&
333 (
dense.IonLow[nelem] <
dense.IonHigh[nelem]-1) )
336 dense.xIonDense[nelem][
dense.IonLow[nelem]+1] +=
338 dense.xIonDense[nelem][
dense.IonLow[nelem]] = 0.;
341 if(
dense.IonLow[nelem] == nelem+1-
NISO )
343 long int ipISO = nelem -
dense.IonLow[nelem];
345 for(
long level = 0; level <
iso_sp[ipISO][nelem].numLevels_max; level++ )
346 iso_sp[ipISO][nelem].st[level].Pop() = 0.;
348 ++
dense.IonLow[nelem];
362 dense.xIonDense[nelem][
dense.IonLow[nelem]-1] == 0. );
365 dense.xIonDense[nelem][
dense.IonLow[nelem]] >=
dense.density_low_limit ||
366 dense.xIonDense[nelem][
dense.IonLow[nelem]]/
dense.gas_phase[nelem] >=
dense.density_low_limit ||
369 (
dense.IonLow[nelem]==
dense.IonHigh[nelem]-1 ));
373 enum {DEBUG_LOC=
false};
374 if( DEBUG_LOC && nelem ==
ipHELIUM )
376 if( lgHi_Down ||lgHi_Up ||lgLo_Up ||lgLo_Down )
378 fprintf(
ioQQQ,
"DEBUG TrimZone\t%li\t",
nzone );
381 fprintf(
ioQQQ,
"high dn %li to %li",
383 dense.IonHigh[nelem] );
387 fprintf(
ioQQQ,
"high up %li to %li",
389 dense.IonHigh[nelem] );
393 fprintf(
ioQQQ,
"low up %li to %li",
395 dense.IonLow[nelem] );
399 fprintf(
ioQQQ,
"low dn %li to %li",
401 dense.IonLow[nelem] );
403 fprintf(
ioQQQ,
"\n" );
410 if(
dense.IonHigh[nelem] < nelem+1 )
414 if( lgHi_Down || lgHi_Up || lgLo_Up || lgLo_Down )
416 conv.lgIonStageTrimed =
true;
419 enum {DEBUG_LOC=
false};
422 fprintf(
ioQQQ,
"DEBUG ion_trim zone\t%.2f \t%li\t",
fnzone, nelem);
424 fprintf(
ioQQQ,
"\tHi_Down");
426 fprintf(
ioQQQ,
"\tHi_Up");
428 fprintf(
ioQQQ,
"\tLo_Up");
430 fprintf(
ioQQQ,
"\tLo_Down");
438 for( ion=0; ion<
dense.IonLow[nelem]; ++ion )
442 for( ion=
dense.IonHigh[nelem]+1; ion<nelem+1; ++ion )
447 for( ion=
dense.IonLow[nelem]; ion<=
dense.IonHigh[nelem]; ++ion )