cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
radius_next.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*radius_next use adaptive logic to find next zone thickness */
4 /*ContRate called by radius_next to find energy of maximum continuum-gas interaction */
5 /*GrainRateDr called by radius_next to find grain heating rate dr */
14 #include "cddefines.h"
15 #include "lines_service.h"
16 #include "iso.h"
17 #include "geometry.h"
18 #include "h2.h"
19 #include "mole.h"
20 #include "hyperfine.h"
21 #include "opacity.h"
22 #include "dense.h"
23 #include "heavy.h"
24 #include "grainvar.h"
25 #include "elementnames.h"
26 #include "conv.h"
27 #include "rfield.h"
28 #include "dynamics.h"
29 #include "thermal.h"
30 #include "hmi.h"
31 #include "coolheavy.h"
32 #include "timesc.h"
33 #include "doppvel.h"
34 #include "stopcalc.h"
35 #include "colden.h"
36 #include "phycon.h"
37 #include "rt.h"
38 #include "trace.h"
39 #include "wind.h"
40 #include "punch.h"
41 #include "taulines.h"
42 #include "pressure.h"
43 #include "iterations.h"
44 #include "struc.h"
45 #include "radius.h"
46 
47 #if 0
48 /*ChkRate called by radius_next to check how rates of destruction of various species changes */
49 STATIC void ChkRate(
50  /* element number on physical scale */
51  long int nelem,
52  /* change in destruction rate */
53  double *dDestRate,
54  /* old and new destruction rates */
55  double *DestRateOld,
56  double *DestRateNew,
57  /* stage of ionization on the physical scale */
58  long int *istage);
59 #endif
60 
61 /*ContRate called by radius_next to find energy of maximum continuum-gas interaction */
62 STATIC void ContRate(double *freqm,
63  double *opacm);
64 
65 /*GrainRateDr called by radius_next to find grain heating rate dr */
66 STATIC void GrainRateDr(double *grfreqm,
67  double *gropacm);
68 
69 /*radius_next use adaptive logic to find next zone thickness
70  * return 0 if ok, 1 for abort */
71 int radius_next(void)
72 {
73  char chLbl[11];
74  bool lgDoPun;
75 
76  long int level , ipStrong;
77 
78  double thickness_total , drThickness , DepthToGo , AV_to_go;
79  int mole_dr_change;
80 
81  double DrGrainHeat,
82  GlobDr,
83  SaveOHIIoHI,
84  SpecDr,
85  Strong,
86  TauDTau,
87  TauInwd,
88  drSolomon_BigH2 ,
89  dEfrac,
90  dHdStep,
91  dRTaue,
92  dTdStep,
93  dnew,
94  drConPres,
95  drH2_heat_cool = 0. ,
96  dH2_heat_cool = 0.,
97  drH2_abund = 0. ,
98  dr_mole_abund = 0.,
99  dH2_abund=0.,
100  dCO_abund=0.,
101  drDepth,
102  drDest,
103  drEfrac,
104  drFail,
105  drFluc,
106  drHMase,
107  drHe1ovHe2,
108  drHion,
109  drInter,
110  drLeiden_hack ,
111  drLineHeat,
112  drSphere,
113  drTab,
114  drdHdStep,
115  drdTdStep,
116  drThermalFront ,
117  drmax,
118  dVelRelative,
119  fac2,
120  factor,
121  freqm,
122  grfreqm=0.,
123  gropacm=0.,
124  hdnew,
125  opacm,
126  recom_dr_last_iter ,
127  OldDR ,
128  winddr,
129  WindAccelDR,
130  x;
131 
132  double change_heavy_frac=-1. , change_heavy_frac_big , dr_change_heavy ,
133  frac_change_heavy_frac_big, Efrac_old , Efrac_new;
134  long int ichange_heavy_nelem=-1 , nelem , ion , ichange_heavy_ion=-1;
135 
136  static double OHIIoHI,
137  OldHeat = 0.,
138  OldTe = 0.,
139  OlddTdStep = 0.,
140  OldH2Abund=0.,
141  OldWindVelocity=0.,
142  Old_He_atom_ov_ion = 0,
143  Old_H2_heat_cool;
144  static long int iteration_last=-1;
145 
146  static double BigRadius = 1e30;
147  bool lgFirstCall;
148 
149  DEBUG_ENTRY( "radius_next()" );
150 
151 
152  /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
153  *
154  * changes in logic
155  * 95 oct 19, drSphere now 3% of radius, was 2%, fewer zone
156  * 95 oct 19, subtracted grain opacity from total opacity used
157  * to get thickness in routine ContRate
158  *
159  *>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
160  *
161  * free statement labels >= 13
162  *
163  *-----------------------------------------------------------------------
164  *
165  * this sub determines the thickness of the next zone
166  * if is called one time for each zone
167  * flag lgNxtDROff is true if this is initialization of radius_next,
168  * is false if we are to use logic to find dr
169  *
170  *----------------------------------------------------------------------- */
171 
172  /* >>chng 03 sep 21 - decide whether this is the first call */
173  if( (iteration != iteration_last) && (nzone==0) )
174  {
175  /* this is the first call in this iteration */
176  iteration_last = iteration;
177  lgFirstCall = true;
178  }
179  else
180  lgFirstCall = false;
181 
182  if( trace.lgTrace )
183  {
184  fprintf( ioQQQ, " radius_next called\n" );
185  }
186 
187  /* save current dr */
188  OldDR = radius.drad;
189 
190  /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
191  /* 1 '' radius_next keys from change in H ionization'',e11.3)')
192  * check on change in hydrogen ionizaiton */
193  if( lgFirstCall )
194  {
195  if( dense.xIonDense[ipHYDROGEN][1] > 0. && dense.xIonDense[ipHYDROGEN][0] > 0. )
196  {
198  }
199  else
200  {
201  OHIIoHI = 0.;
202  }
203  SaveOHIIoHI = OHIIoHI;
204  drHion = BigRadius;
205  /* else if(hii.gt.0. .and. hi.gt.0. .and. OHIIoHI.gt.0. ) then
206  * >>chng 97 jul 9, for deep in PDR vastly now ionz H slowed down works */
207  }
208 
209  else if( (dense.xIonDense[ipHYDROGEN][1] > 0. && dense.xIonDense[ipHYDROGEN][0] > 0.) && OHIIoHI > 1e-15 )
210  {
211  double atomic_frac = (dense.xIonDense[ipHYDROGEN][1]/dense.xIonDense[ipHYDROGEN][0]);
212  /* ratio of current HII/HI to old value - < 1 when becoming more neutral */
213  /* this is now change in atomic fraction */
214  x = 1. - atomic_frac /OHIIoHI;
215  if( atomic_frac > 0.05 && atomic_frac < 0.9 )
216  {
217  /* >>chng 96 feb 23 from 0.7 to 0.3 because punched through H i-front
218  * >>chng 97 may 5, from 0.3 to 0.2 since punched through H i-front */
219  /* >>chng 02 jun 05 from 0.2 to 0.05 poorly resolved i-front, also added two-branch logic*/
220  drHion = radius.drad*MAX2( 0.2 , 0.05/MAX2(1e-10,x) );
221  }
222  else if( x > 0. )
223  {
224  /* >>chng 96 feb 23 from 0.7 to 0.3 because punched through H i-front
225  * >>chng 97 may 5, from 0.3 to 0.2 since punched through H i-front */
226  drHion = radius.drad*MAX2( 0.2 , 0.2/MAX2(1e-10,x) );
227  }
228  else
229  {
230  drHion = BigRadius;
231  }
232  SaveOHIIoHI = OHIIoHI;
234  }
235 
236  else
237  {
238  SaveOHIIoHI = OHIIoHI;
239  if( dense.xIonDense[ipHYDROGEN][1] > 0. && dense.xIonDense[ipHYDROGEN][0] > 0. )
240  {
242  }
243  else
244  {
245  OHIIoHI = 0.;
246  }
247  drHion = BigRadius;
248  }
249 
250  /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
251  /* "radius_next keys from H maser, dt, ij=" possible hydrogen maser action */
252  if( rt.dTauMase < -0.01 )
253  {
254  /* maser so powerful, must limit inc in tay
255  * >>chng 96 aug 08, from 0.3 to 0.1 due to large maser crash */
256  drHMase = radius.drad*MAX2(0.1,-0.2/rt.dTauMase);
257  }
258  else
259  {
260  drHMase = BigRadius;
261  }
262 
263  /* >>chng 03 nov 09, try doing he in following, not above */
264  Old_He_atom_ov_ion = 0.;
265  drHe1ovHe2 = BigRadius;
266 
267  /* check on N0 - 1 ionization changes,
268  * >>chng 03 jun 06, add this test due to smashing into H ifront in blr89.in
269  */
270  dr_change_heavy = BigRadius;
271  change_heavy_frac_big = -1.;
272  frac_change_heavy_frac_big = -1.;
273  /* >chng 03 jun 09, from 0.05 to 0.1, initial tests with zoning */
274 # define CHANGE_ION_HEAV 0.2f
275 # define CHANGE_ION_HHE 0.15f
276  if( nzone > 4 )
277  {
278  for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
279  {
280  if( dense.lgElmtOn[nelem] )
281  {
282  realnum change;
283  /* this is the limit to the ionization we will check -
284  * also used in prt_comment to check on whether oscillations occurred */
285  realnum frac_limit;
286  if( nelem<=ipHELIUM || nelem==ipCARBON )
287  {
288  frac_limit = 1e-4f;
289  change = CHANGE_ION_HHE;
290  }
291  else
292  {
293  /* this var is used to print warnings,
294  * make sure we converge below it */
295  frac_limit = struc.dr_ionfrac_limit/2.f;
296  change = CHANGE_ION_HEAV;
297  }
298  /* >>chng 03 dec 09, go up to full range of ion, not just =2 */
299  for( ion=0; ion<=nelem+1; ++ion )
300  {
301  realnum abundnew = dense.xIonDense[nelem][ion]/SDIV( dense.gas_phase[nelem]);
302  realnum abundold = struc.xIonDense[nelem][ion][nzone-3]/SDIV( struc.gas_phase[nelem][nzone-3]);
303  if( abundnew > frac_limit && abundold > frac_limit )
304  {
305  /* NB must make sure this test is not done when nzone-x <0 */
306  /* -2 because previous zone, and nzone is off by one (on physical, not C, scale) */
307  /*realnum abundold = struc.xIonDense[nelem][ion][nzone-3]/SDIV( struc.gas_phase[nelem][nzone-3]);*/
308  realnum abundolder = struc.xIonDense[nelem][ion][nzone-4]/SDIV( struc.gas_phase[nelem][nzone-4]);
309  realnum abundoldest = struc.xIonDense[nelem][ion][nzone-5]/SDIV( struc.gas_phase[nelem][nzone-5]);
310  /* this is fractional change */
311  /* >>chng 04 feb 28, take min of old and new abund, to try to pick up
312  * rapid changing Ca+ sooner */
313  change_heavy_frac = fabs(abundnew-abundold)/MIN2(abundold,abundnew);
314  /* want fractional change to be less than this factor */
315  if( (change_heavy_frac > change) && (change_heavy_frac > change_heavy_frac_big) &&
316  /* >>chng 03 dec 07, add test that abund is not oscillating */
317  /* also test that abundance is increasing - we are headed into a front */
318  ((abundnew-abundold)>0.) &&
319  ((abundold-abundolder)>0.) &&
320  ((abundolder-abundoldest)>0.) )
321  {
322  ichange_heavy_nelem = nelem;
323  ichange_heavy_ion = ion;
324  change_heavy_frac_big = change_heavy_frac;
325  frac_change_heavy_frac_big = abundnew;
326  /* >>chng 03 dec 06, from min of 0.5 to min of 0.25, crash into He i-front
327  * in hizqso.in */
328  /* >>chng 04 mar 03, min had become 0.1, forced oscillations in nova.in
329  * in silicon, zoning changed greatly, causing change in diffuse lin
330  * pumping. put back to 0.25 */
331  dr_change_heavy = radius.drad * MAX2(0.25, change / change_heavy_frac );
332  }
333  }
334  }
335  }
336  }
337  }
338 
339  /* if Leiden hacks are on then use increase in dust extinction as control
340  * >>chng 05 aug 13, add this */
341  if(!co.lgUMISTrates)
342  {
343  /* Draine field is only defined over narrow range in FUV - must not let change
344  * in extinction become too large -
345  * prefactor is change in optical depth */
346  drLeiden_hack = MAX2( 0.02 , 0.05*rfield.extin_mag_V_point) / SDIV( geometry.FillFac * rfield.opac_mag_V_point );
347  }
348  else
349  {
350  drLeiden_hack = BigRadius;
351  }
352  /* >>chng 04 feb 15, kill this block - not used */
353 
354  /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
355  /* check how heating is changing
356  * '' radius_next keys from change in heating; current, delta='', */
357  if( nzone <= 1 || thermal.lgTemperatureConstant )
358  {
359  drdHdStep = BigRadius;
360  dHdStep = FLT_MAX;
361  }
362  else
363  {
364  dHdStep = fabs(thermal.htot-OldHeat)/thermal.htot;
365  if( dHdStep > 0. )
366  {
367  if( dense.gas_phase[ipHYDROGEN] >= 1e13 )
368  {
369  /* drdHdStep = drad * MAX( 0.8 , 0.05/dHdStep ) */
370  drdHdStep = radius.drad*MAX2(0.8,0.075/dHdStep);
371  }
372  else if( dense.gas_phase[ipHYDROGEN] >= 1e11 )
373  {
374  /* drdHdStep = drad * MAX( 0.8 , 0.075/dHdStep ) */
375  drdHdStep = radius.drad*MAX2(0.8,0.100/dHdStep);
376  }
377  else
378  {
379  /* changed from .15 to .12 for outer edge of coolhii too steep dT
380  * changed to .10 for symmetry, big change in some rates, 95nov14
381  * changed from .10 to .125 since parispn seemed to waste zones
382  * >>chng 96 may 21, from .125 to .15 since pn's still waste zones */
383  drdHdStep = radius.drad*MAX2(0.8,0.15/dHdStep);
384  }
385  }
386  else
387  {
388  drdHdStep = BigRadius;
389  }
390  }
391  OldHeat = thermal.htot;
392 
393  /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
394  /* pressure due to incident continuum if in eos */
395  if( strcmp(dense.chDenseLaw,"CPRE") == 0 && pressure.lgContRadPresOn )
396  {
397  if( nzone > 2 && pressure.pinzon > 0. )
398  {
399  /* pinzon is pressrue from acceleration onto previos zone
400  * want this less than some fraction of total pressure */
401  /* >>chng 06 feb 01, change from init pres to current total pressure
402  * in const press high U ulirgs current pressure may be quite larger
403  * than init pressure due to continuum absorption */
404  drConPres = 0.05*pressure.PresTotlCurr/(wind.AccelTot*
406  }
407  else
408  {
409  drConPres = BigRadius;
410  }
411  }
412  else
413  {
414  drConPres = BigRadius;
415  }
416 
417  /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
418  /* check how temperature is changing
419  * 1 '' radius_next keys from change in temperature; current, delta='', */
420  if( nzone <= 1 )
421  {
422  drdTdStep = BigRadius;
423  dTdStep = FLT_MAX;
424  OlddTdStep = dTdStep;
425  }
426  else
427  {
428  /* change in temperature; current= */
429  dTdStep = (phycon.te-OldTe)/phycon.te;
430  /* >>chng 02 dec 08, desired change in temperature per zone must not
431  * be smaller than allower error in temperature. For now use relative error
432  * in heating - - cooling balance. Better would be to also use c-h deriv wrt t
433  * to get slope */
435  x = MAX2( 0.01 , x );
436  x = MIN2( 0.05 , x );
437  /* >>chng 02 dec 11 rjrw change back to 0.03 -- improve dynamics.dRad criterion instead */
438  x = 0.03;
439  /* >>chng 02 dec 07, do not do this if there is mild te jitter,
440  * so that dT is changing sign - this happens
441  * in ism.ini, very stable temperature with slight noise up and down */
442  if( dTdStep*OlddTdStep > 0. )
443  {
444  /* >>chng 96 may 30, variable depending on temp
445  * >>chng 96 may 31, allow 0.7 smaller, was 0.8
446  * >>chng 97 may 05, from 0.7 to 0.5 stop from punching through thermal front
447  * >>chng 04 mar 30, from 0.7 to 0.5 stop from punching through thermal front,
448  * for some reason factor was 0.7, not 0.5 as per previous change */
449  double absdt = fabs(dTdStep);
450  drdTdStep = radius.drad*MAX2(0.5,x/absdt);
451  }
452  else
453  {
454  drdTdStep = BigRadius;
455  }
456  }
457  OlddTdStep = dTdStep;
458  OldTe = phycon.te;
459 
460  /* >>chng 02 oct 06, only check on opacity - interaction if not
461  * constant temperature - there were constant temperature models that
462  * extended to infinity but hung with last few photons and this test.
463  * better to ignore this test which is really for thermal models */
464  /* >>chng 06 mar 20, do not call if recombination logic in place */
466  {
467  /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
468  /* find freq where opacity largest and interaction rate is fastest
469  * "cont inter nu=%10.2e opac=%10.2e\n" */
470  ContRate(&freqm,&opacm);
471 
472  /* use optical depth at max interaction energy
473  * >>chng 96 jun 06 was drChange=0.15 changed to 0.3 for high Z models
474  * taking too many zones
475  * drInter = drChange / MAX(1e-30,opacm*FillFac) */
476 
477  drInter = 0.3/MAX2(1e-30,opacm*geometry.FillFac*geometry.DirectionalCosin);
478  }
479  else
480  {
481  drInter = BigRadius;
482  freqm = 0.;
483  opacm = 1.;
484  }
485 
486 # if 0
487  /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
488  /* check on changes in destruction rates for various atoms */
489  ChkRate(ipCARBON,&dDRCarbon,&DestOldCarb, &DestNewCarb,&icarstag);
490  ChkRate(ipNITROGEN,&dDRNitrogen,&DestOldNit, &DestNewNit,&initstag);
491  ChkRate(ipOXYGEN,&dDROxygen,&DestOldOxy, &DestNewOxy,&ioxystag);
492  ChkRate(ipIRON,&dDRIron,&DestOldIron, &DestNewIron,&iironstag);
493 
494  dDestRate = MAX4(dDROxygen,dDRIron,dDRCarbon,dDRNitrogen);
495 
496  if( fp_equal( dDRCarbon, dDestRate ) )
497  {
498  dDestRate = dDRCarbon;
499  DestOld = DestOldCarb;
500  DestNew = DestNewCarb;
501  istage = icarstag;
502  strcpy( chDestAtom, "Carbon " );
503  }
504 
505  else if( fp_equal( dDRNitrogen, dDestRate ) )
506  {
507  dDestRate = dDRNitrogen;
508  DestOld = DestOldNit;
509  DestNew = DestNewNit;
510  istage = initstag;
511  strcpy( chDestAtom, "Nitrogen" );
512  }
513 
514  else if( fp_equal( dDROxygen, dDestRate ) )
515  {
516  dDestRate = dDROxygen;
517  DestOld = DestOldOxy;
518  DestNew = DestNewOxy;
519  strcpy( chDestAtom, "Oxygen " );
520  istage = ioxystag;
521  }
522 
523  else if( fp_equal( dDRIron, dDestRate ) )
524  {
525  dDestRate = dDRIron;
526  DestOld = DestOldIron;
527  DestNew = DestNewIron;
528  istage = iironstag;
529  strcpy( chDestAtom, "Iron " );
530  }
531 
532  else
533  {
534  fprintf( ioQQQ, " insanity following ChkRate\n" );
535  ShowMe();
536  cdEXIT(EXIT_FAILURE);
537  }
538 
539  /* radius_next keys from change in dest rates, atom= */
540  if( dDestRate > 0. )
541  {
542  /* if( te.gt.40 000. ) then
543  * added different accuracy for hot gas since tend to jump over
544  * big te range for small change in heat (intrinsically unstable)
545  * drDest = drad * MAX( 0.5 , 0.10/dDestRate )
546  * else
547  * was drChange, changed to .15 for parishii go through HeII=HeI I front
548  * drDest = drad * MAX( 0.5 , 0.15/dDestRate )
549  * >>chng 95 dec 18 from above to below to stop oscillations
550  * >>chng 95 dec 27 from min of .5 to .75 to stop zone size from changing rapidly
551  * drDest = drad * MAX( 0.8 , 0.20 /dDestRate )
552  * >>chng 96 jan 14 from .2 to .25 to use less zones
553  * >>chng 96 may 30 from min of 0.8 to 0.5 to prevent crashing into He i-front */
554  drDest = radius.drad*MAX2(0.5,0.20/dDestRate);
555  /* endif */
556  }
557  else
558  {
559  drDest = BigRadius;
560  }
561 # endif
562  drDest = BigRadius;/*03 dec 12 is this needed? */
563 
564  /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
565  /* check whether change in wind velocity constrains DRAD
566  * WJH 22 May 2004: disable when we are near the sonic point since
567  * the velocity may be jumping all over the place but we just want
568  * to push through it as quickly as we can */
570  {
571  double v = fabs(wind.windv);
572  /* this is relative change in velocity */
573  dVelRelative = fabs(wind.windv-OldWindVelocity)/
575 
576 # define WIND_CHNG_VELOCITY_RELATIVE 0.01
577  /*fprintf(ioQQQ,"DEBUG rad %.3f vel %.2e dVelRelative %.2e",
578  log10(radius.Radius) , wind.windv , dVelRelative );*/
579 
580  /* do not use this on first zone since do not have old velocity */
581  if( dVelRelative > WIND_CHNG_VELOCITY_RELATIVE && nzone>1 )
582  {
583  /* factor less than one if change larger than WIND_CHNG_VELOCITY_RELATIVE */
584  double factor = WIND_CHNG_VELOCITY_RELATIVE / dVelRelative;
585  /*fprintf(ioQQQ," DEBUG factor %.2e", factor );*/
586  factor = MAX2(0.8 , factor );
587  winddr = radius.drad * factor;
588  }
589  else
590  {
591  winddr = BigRadius;
592  }
593  /*fprintf(ioQQQ," DEBUG \n" );*/
594 
595  WindAccelDR = BigRadius;
596 
597  /* set dr from advective term,
598  * the 1/20 came from looking at one set of structure plots */
600  {
601  /* >>chng 02 dec 11, from 5 to 20 */
602  winddr = MIN2( winddr , dynamics.dRad / 20. );
603  /*>>chng 04 oct 06, set dVelRelative to dynamics.dRad since dVelRelative is printed as part
604  * of reason for choosing this criteria, want to reflect valid reason. */
605  dVelRelative = dynamics.dRad;
606  }
607  }
608  else
609  {
610  winddr = BigRadius;
611  WindAccelDR = BigRadius;
612  dVelRelative = 0.;
613  }
614  /* remember old velocity */
615  OldWindVelocity = wind.windv;
616 
617  /* inside out globule */
618  if( strcmp(dense.chDenseLaw,"GLOB") == 0 )
619  {
620 # define DNGLOB 0.10
621  if( radius.glbdst < 0. )
622  {
623  fprintf( ioQQQ, " Globule distance is negative, internal overflow has occured, sorry.\n" );
624  fprintf( ioQQQ, " This is routine radius_next, GLBDST is%10.2e\n",
625  radius.glbdst );
626  cdEXIT(EXIT_FAILURE);
627  }
630  if( fac2/factor > 1. + DNGLOB )
631  {
632  /* DNGLOB is relative change in density allowed this zone, 0.10 */
633  GlobDr = radius.drad*DNGLOB/(fac2/factor - 1.);
634  }
635  else
636  {
637  GlobDr = BigRadius;
638  }
639  /* GlobDr = GLBDST * DNGLOB * (GLBRAD/GLBDST)**(-GLBPOW) / GLBPOW */
640  GlobDr = MIN2(GlobDr,radius.glbdst/20.);
641  }
642  else
643  {
644  GlobDr = BigRadius;
645  }
646 
647  /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
648  hdnew = 0.;
649  if( strncmp( dense.chDenseLaw , "DLW" , 3) == 0 )
650  {
651  /* one of the special density laws, first get density at possible next dr */
652  if( strcmp(dense.chDenseLaw,"DLW1") == 0 )
653  {
655  radius.drad);
656  }
657  else if( strcmp(dense.chDenseLaw,"DLW2") == 0 )
658  {
660  radius.drad);
661  }
662  else
663  {
664  fprintf( ioQQQ, " dlw insanity in radius_next\n" );
665  cdEXIT(EXIT_FAILURE);
666  }
667  drTab = fabs(hdnew-dense.gas_phase[ipHYDROGEN])/MAX2(hdnew,dense.gas_phase[ipHYDROGEN]);
668  drTab = radius.drad*MAX2(0.2,0.10/MAX2(0.01,drTab));
669  }
670  else
671  {
672  drTab = BigRadius;
673  }
674 
675  /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
676  /* special density law */
677  if( strcmp(dense.chDenseLaw,"DLW1") == 0 )
678  {
681  /* DNGLOB is relative change in density allowed this zone, 0.10 */
682  if( dnew == 0. )
683  {
684  SpecDr = radius.drad*3.;
685  }
686  else if( dnew/DNGLOB > 1.0 )
687  {
688  SpecDr = radius.drad/(dnew/DNGLOB);
689  }
690  else
691  {
692  SpecDr = BigRadius;
693  }
694  }
695  else
696  {
697  SpecDr = BigRadius;
698  }
699 
700  /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
701  /* check grain line heating dominates
702  * this is important in PDR and HII region calculations
703  * >>chng 97 jul 03, added following check */
704  if( thermal.heating[0][13]/thermal.htot > 0.2 )
705  {
706  /* >>chng 01 jan 03, following returns 0 when NO light at all,
707  * had failed with botched assert */
708  GrainRateDr(&grfreqm,&gropacm);
709  DrGrainHeat = 1.0/MAX2(1e-30,gropacm*geometry.FillFac*geometry.DirectionalCosin);
710  }
711  else
712  {
713  gropacm = 0.;
714  grfreqm = 0.;
715  DrGrainHeat = BigRadius;
716  }
717 
718  /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
719  /* check if line heating dominates
720  * this is important in high metallicity models */
721  if( thermal.heating[0][22]/thermal.htot > 0.2 )
722  {
723  FndLineHt(&level,&ipStrong,&Strong);
724  if( Strong/thermal.htot > 0.1 )
725  {
726  if( level == 1 )
727  {
728  /* a level1 line was the heat source (usual case)
729  * drLineHeat = MAX(1.0,TauLines(ipLnTauIn,ipStrong)*0.2) /
730  * 1 TauLines(ipLnDTau,ipStrong) */
731  TauInwd = TauLines[ipStrong].Emis->TauIn;
732  TauDTau = TauLines[ipStrong].Emis->PopOpc * TauLines[ipStrong].Emis->opacity /
733  DoppVel.doppler[TauLines[ipStrong].Hi->nelem-1];
734  }
735  else if( level == 2 )
736  {
737  /* a atom_level2 line was the heat source
738  * (bad case since not as well treated)
739  * drLineHeat = MAX(1.0,WindLine(ipLnTauIn,ipStrong)*0.2) /
740  * 1 WindLine(ipLnDTau,ipStrong) */
741  TauInwd = TauLine2[ipStrong].Emis->TauIn;
742  TauDTau = TauLine2[ipStrong].Emis->PopOpc * TauLine2[ipStrong].Emis->opacity /
743  DoppVel.doppler[TauLine2[ipStrong].Hi->nelem-1];
744  }
745  else if( level == 3 )
746  {
747  /* a 12CO line was the heat source
748  * (bad case since not as well treated)
749  * drLineHeat = MAX(1.0,WindLine(ipLnTauIn,ipStrong)*0.2) /
750  * 1 WindLine(ipLnDTau,ipStrong) */
751  TauInwd = C12O16Rotate[ipStrong].Emis->TauIn;
752  TauDTau = C12O16Rotate[ipStrong].Emis->PopOpc * C12O16Rotate[ipStrong].Emis->opacity /
753  DoppVel.doppler[C12O16Rotate[ipStrong].Hi->nelem-1];
754  }
755  else if( level == 4 )
756  {
757  /* a 13CO line was the heat source
758  * (bad case since not as well treated)
759  * drLineHeat = MAX(1.0,WindLine(ipLnTauIn,ipStrong)*0.2) /
760  * 1 WindLine(ipLnDTau,ipStrong) */
761  TauInwd = C13O16Rotate[ipStrong].Emis->TauIn;
762  TauDTau = C13O16Rotate[ipStrong].Emis->PopOpc * C13O16Rotate[ipStrong].Emis->opacity /
763  DoppVel.doppler[C13O16Rotate[ipStrong].Hi->nelem-1];
764  }
765  else if( level == 5 )
766  {
767  /* >>chng 03 dec 07, this branch had been left off, caught by Hiroaki Oyaizu */
768  /* a hyperfine transition */
769  TauInwd = HFLines[ipStrong].Emis->TauIn;
770  TauDTau = HFLines[ipStrong].Emis->PopOpc * HFLines[ipStrong].Emis->opacity /
771  DoppVel.doppler[HFLines[ipStrong].Hi->nelem-1];
772  }
773  else
774  {
775  /* this is insane, since Strong was set, but not level */
776  fprintf( ioQQQ, " PROBLEM radius_next Strong line heating set, but not level.\n" );
777  TotalInsanity();
778  }
779  /* in following logic cannot use usual inverse opacity,
780  * since line heating competes with escape probability,
781  * so is effective at surprising optical depths */
782  if( TauDTau > 0. )
783  {
784  drLineHeat = MAX2(1.,TauInwd)*0.4/TauDTau;
785  }
786  else
787  {
788  drLineHeat = BigRadius;
789  }
790  }
791  else
792  {
793  TauInwd = 0.;
794  drLineHeat = BigRadius;
795  ipStrong = 0;
796  Strong = 0.;
797  }
798  }
799  else
800  {
801  TauInwd = 0.;
802  drLineHeat = BigRadius;
803  ipStrong = 0;
804  level = 0;
805  Strong = 0.;
806  }
807 
808  /* >>chng 03 mar 03, add this logic */
809  /* do not let change in cooling/heating due to H2 become too large */
810  drH2_heat_cool = BigRadius;
811  if( lgFirstCall )
812  {
813  Old_H2_heat_cool = 0.;
814  }
815  else if( !thermal.lgTemperatureConstant )
816  {
817  /* this is case where temperature has not been set */
818  /* compare total heating - cooling due to h2 with total due to everything */
819  double H2_heat_cool = (fabs(hmi.HeatH2Dexc_used)+fabs(hmi.HeatH2Dish_used)) / thermal.htot;
820  if( H2_heat_cool > 0.1 )
821  {
822  dH2_heat_cool = fabs( H2_heat_cool - Old_H2_heat_cool );
823  dH2_heat_cool = SDIV(dH2_heat_cool);
824  /* >>chng 05 oct 27, had been taking 20% of original radius - this caused zoning
825  * to become very fine and may not have been needed - change from 0.2 to 0.5 */
826  /*drH2_heat_cool = radius.drad*MAX2(0.2,0.05/dH2_heat_cool);*/
827  drH2_heat_cool = radius.drad*MAX2(0.5,0.05/dH2_heat_cool);
828  }
829  else
830  {
831  drH2_heat_cool = BigRadius;
832  }
833  }
834  Old_H2_heat_cool = (fabs(hmi.HeatH2Dexc_used)+fabs(hmi.HeatH2Dish_used)) / thermal.htot;
835 
836  /* >>chng 03 mar 04, add this logic */
837  /* do not let change in H2 and heavy element molecular abundances become too large
838  * "change in heav ele mole abundance, d(mole)/elem" */
839  drH2_abund = BigRadius;
840  dr_mole_abund = BigRadius;
841  mole_dr_change = -1;
842  if( nzone>=4 )
843  {
844  /* first do H2 abundance */
845  double Old_H2_abund;
846  /* >>chng 04 dec 15, do not use special logic when large h2 is turned on */
847  int i;
848  Old_H2_abund = struc.H2_molec[ipMH2g][nzone-3] + struc.H2_molec[ipMH2s][nzone-3];
849  /* >>chng 03 jun 16, limit from 0.01 to 0.001, some models fell over H2 front due to
850  * large zoning, when large H2 was just this caused oscillations in solomon process */
851  /* >>chng 03 nov 22, from > 0.001 to > 3e-4, models that start almost in H2 have
852  * rapid increase in H2 at shallow depths, start sensing this sooner */
853  /* >>chng 03 dec 10, from 3e-4 to 1e-4, to get smaller chagned in leiden1.in test */
854  /* radius_next keys from change in H2 abundance, d(H2)/H */
855  /* >>chng 04 mar 11, start sensing H2 earlier since otherwise step size
856  * needs to become way too small tooo quickly. limit changed from 1e-4 to 1e-6 */
857  /* >>chng 04 jun 29, fromo > 1e-6 to >1e-8, some pdr's had too large a change in H2 */
858 
859  if( 2.*hmi.H2_total/dense.gas_phase[ipHYDROGEN] > 1e-8 )
860  {
861  double fac = 0.2;
862  /* this is percentage change in H2 density - "change in H2 abundance" */
863  dH2_abund = 2.*fabs( hmi.H2_total - Old_H2_abund ) / hmi.H2_total;
864  /* in testing th85ism the dH2_abund did come out zero */
865  /* >>chng 03 jun 16, change d(H2) from 0.05 to 0.1, fine resolution of H2/H exposed
866  * small numerical oscillations in Solomon process */
867  /* >>chng 03 nov 22, smallest possible ratio of dr(next)/dr changed from
868  * 0.2 to 0.05, models that started almost in H2 front need to be able to sense it */
869  /*drH2_abund = radius.drad*MAX2(0.05,fac/SDIV(dH2_abund) );*/
870  /* >>chng 04 mar 15, with such small possible changes in dr, only 0.05, a thermal front
871  * can easily cause large changes in H2 and T that are not due to zoning, but to the
872  * discontinuity. make smallest change larger to prevent hald due to dr */
873  /* >>chng 05 aug 23, thermal front allowed dr to become much too small
874  * chng from 0.02 to .6 */
875  dH2_abund = SDIV(dH2_abund);
876  drH2_abund = radius.drad*MAX2(0.2,fac/dH2_abund );
877  }
878  else
879  {
880  drH2_abund = BigRadius;
881  }
882 
883  /* check how molecular fractions of all heavy elements are chaning relative
884  * to total gas phase abundance */
885  dr_mole_abund = BigRadius;
886  /* >>chng 04 jun 02, upper limit had been all species but now limit to real
887  * molecules since do not want this logic to work with the ions */
888  for(i=0; i<mole.num_comole_calc; ++i)
889  {
890  realnum abund,
891  abund_limit;
892  if( !dense.lgElmtOn[COmole[i]->nelem_hevmol] || COmole[i]->n_nuclei <= 1 )
893  continue;
894  /* >>chng 03 sep 21, add CO logic */
895  /* >>chng 04 mar 30, generalize to any molecule at all */
896  /* >>chng 04 mar 31 lower limit to abund had been 0.01, change
897  * to 0.001 to pick up approach to molecular fronts */
898  abund = COmole[i]->hevmol/dense.gas_phase[COmole[i]->nelem_hevmol];
899  /* is this an ice? need special logic for ice since density increases
900  * exponentially when grain temp falls below sublimation temperature
901  * >>chng 05 dec 06 - detect changes for smaller abundances for ices
902  * due to large changes in ice abundances */
903  if( COmole[i]->lgGas_Phase )
904  {
905  abund_limit = 1e-3f;
906  }
907  else
908  {
909  /* this is an ice - track its abundance at lower values so that
910  * we resolve the sublimation transition region */
911  abund_limit = 1e-5f;
912  }
913 
914  if( abund > abund_limit )
915  {
916  double drmole_one, relative_change, relative_change_denom;
917  /* >>chng 05 dec 08, use smaller abundance for the denominator since just taking
918  * current abundance will overlook case where current density is vastly large
919  * than old density */
920  if( struc.CO_molec[i][nzone-3]>SMALLFLOAT )
921  {
922  relative_change_denom = MIN2( COmole[i]->hevmol , struc.CO_molec[i][nzone-3] );
923  }
924  else
925  {
926  relative_change_denom = COmole[i]->hevmol;
927  }
928  /* the relative change in the abundance */
929  relative_change =
930  /*fabs( COmole[i]->hevmol - struc.CO_molec[i][nzone-3] ) / COmole[i]->hevmol;*/
931  /* >>chng 05 dec 08, use smaller abundance */
932  fabs( COmole[i]->hevmol - struc.CO_molec[i][nzone-3] ) / relative_change_denom;
933  /* >>chng 03 jun 16, change from 0.05 to 0.1, fine resolution of H2/H exposed
934  * small numerical oscillations in Solomon process */
935  /* >>chng 04 jun 02, from 0.1 back to 0.05, more extensive CO etc network
936  * caused oscillations in SiO abundance and Si Si+ density. */
937  /* >>chng 04 aug 03, from 0.05 to 0.035, leiden pdr model v2 had major
938  * jump in eden */
939  /* >>chng 04 oct 18, from 0.035 back to 0.05, leiden pdr v2 actually due to having
940  * PAHs in fully molecular limit (??), this caused cool flow pdr grid to trip on
941  * too small dr */
942  relative_change = SDIV(relative_change);
943  /*>>chng 05 dec 08, relative_change must be less than one - with
944  * revised logic above can be bigger than one */
945  if( relative_change > 1. )
946  relative_change = 1./relative_change;
947  /*drmole_one = radius.drad*MAX2(0.2,0.035/relative_change );*/
948  /* >>chng 05 aug 23, thermal front allowed dr to become much too small
949  * chng from 0.02 to .6 */
950  /*drmole_one = radius.drad*MAX2(0.2,0.05/relative_change );*/
951  drmole_one = radius.drad*MAX2(0.6,0.05/relative_change );
952  /* final dr will be the smallest we encounter */
953  if( drmole_one < dr_mole_abund )
954  {
955  /* this is the dr used to set next dr - keep track of which moe was changing */
956  dr_mole_abund = drmole_one;
957  mole_dr_change = i;
958  dCO_abund = relative_change;
959  }
960  }
961  }
962  }
963 
964  /* some consideration due to big H2 molecule */
965  drSolomon_BigH2 = H2_DR();
966 
967  /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
968  /* can't make drmax large deep in model since causes feedback
969  * oscillations with changes in heating or destruction rates
970  * >>chng 96 oct 15, change from 5 to 10 */
971  if( nzone < 5 )
972  {
973  /* >>chng 96 nov 11, had been 4 * drad up to 11, change to following
974  * to be similar to c90.01, max of 1.3 between 5 and 11
975  * >>chng 04 oct 29 geometry.DirectionalCosin was ioncorrect applied to this factor */
976  drmax = 4.*radius.drad;
977  }
978  else
979  {
980  drmax = 1.3*radius.drad;
981  }
982 
983  /* >>chng 05 apr 05, do not sense temp oscillation, so that we can move past this
984  * point if it occurs */
985 # if 0
986  /* look for oscillations in electron density or tempeature - freeze dr if these occur */
987  dr_ne_oscil = BigRadius;
988  dr_te_oscil = BigRadius;
989  if( nzone >= 11 )
990  {
991  /* >>chng 96 oct 15, do not let zones increase if oscillations present */
992  /* >>chng 96 oct 31, error to declare oscillation propto toler, the
993  *heating cooling tolerance */
996  limit = nzone -2;
997  ASSERT( limit < struc.nzlim );
998  for( k=nzone - 10; k < limit; k++ )
999  {
1000  /* check that square of change both chng sign and is
1001  * greater than square of heat-tool error */
1002  if( (struc.testr[k-1] - struc.testr[k])/struc.testr[k]*
1003  (struc.testr[k] - struc.testr[k+1])/struc.testr[k] <
1004  -errorHC )
1005  {
1006  dr_te_oscil = radius.drad;
1007  }
1008  /* small residiual is to allow 0.01 rel error */
1009  if( (struc.ednstr[k-1] - struc.ednstr[k])/struc.ednstr[k]*
1010  (struc.ednstr[k] - struc.ednstr[k+1])/struc.ednstr[k] <
1011  -errorNe )
1012  {
1013  /* >>chng 96 oct 15, do not let zones increase if oscillations present */
1014  /* radius_next keys from electron density oscillation*/
1015  dr_ne_oscil = radius.drad;
1016  }
1017  }
1018  }
1019  /* >>chng 05 apr 05, do not sense temp oscillation, so that we can move past this
1020  * point if it occurs */
1021  dr_ne_oscil = BigRadius;
1022  dr_te_oscil = BigRadius;
1023 # endif
1024 
1025  /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
1026  /* check on several convergence criteria */
1027 
1028  /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
1029  if( !conv.lgConvTemp )
1030  {
1031  drFail = radius.drad/2.;
1032  }
1033  else
1034  {
1035  drFail = BigRadius;
1036  }
1037 
1038  /* time dependent recombination - conditions become very homogeneous and
1039  * crash into I fronts - use last iteration's radius as guess of current case */
1040  recom_dr_last_iter = BigRadius;
1042  {
1043  /* first time through nzone == 1 */
1044  static long int nzone_recom = -1;
1045  if( nzone<=1 )
1046  {
1047  /* initialization case */
1048  nzone_recom = 0;
1049  }
1050  else if( radius.depth < struc.depth_last[struc.nzone_last-1] &&
1051  nzone_recom < struc.nzone_last )
1052  {
1053  ASSERT( nzone_recom>=0 && nzone_recom<struc.nzone_last );
1054  /* case where we are within previous computed structure
1055  * first possibly increase nzone_recom so that it points deeper
1056  * than this zone */
1057  while( struc.depth_last[nzone_recom] < radius.depth &&
1058  nzone_recom < struc.nzone_last-1 )
1059  ++nzone_recom;
1060  recom_dr_last_iter = struc.drad_last[nzone_recom]*3.;
1061  }
1062  else
1063  {
1064  /* case where we have overrun the previous iteration's saved structure */
1065  recom_dr_last_iter = BigRadius;
1066  }
1067  }
1068 
1069  /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
1070  /* change in electron density - radius_next keys from change in elec den,
1071  * remember old electron density during first call */
1072  /* this is low ionization solution */
1073  if( nzone > 2 )
1074  {
1075  /* next is-2 since nzone is on physics not c scale, and we want previous zone */
1076  Efrac_old = struc.ednstr[nzone-3]/struc.gas_phase[ipHYDROGEN][nzone-3];
1077  Efrac_new = dense.eden/dense.gas_phase[ipHYDROGEN];
1078  dEfrac = fabs(Efrac_old-Efrac_new) / Efrac_new;
1079 
1080  if( dEfrac > SMALLFLOAT )
1081  {
1082  double fac = 0.04;
1083  /* >>chng 03 dec 25, use smaller rel change in elec frac when most elec in ipMole or grains */
1084  /* >>chng 04 sep 14, change to from from metals but comment out */
1085  /* >>chng 04 sep 17, change to from from metals - uncomment */
1086  if( dense.eden_from_metals > 0.5 )
1087  {
1088  /* >>chng 04 sep 18, change 0.02 from 0.01 */
1089  /* >>chng 04 sep 18, change 0.02 from 0.04 */
1090  fac = 0.04;
1091  }
1092  /* >>chng 04 feb 23, add test on hydrogen being predomintly
1093  * recombined due to three-body recom, which is very sensitive
1094  * to the electron density - but only do this in partially ionized medium */
1095  else if( iso.RecomCollisFrac[ipH_LIKE][ipHYDROGEN] > 0.8 &&
1098 
1099  {
1100  fac = 0.02;
1101  }
1102  /* >>chng 04 sep 17, change to 0.1 from 0.2 */
1103  drEfrac = radius.drad*MAX2(0.1,fac/dEfrac);
1104  }
1105  else
1106  {
1107  drEfrac = BigRadius;
1108  }
1109  }
1110  else
1111  {
1112  dEfrac = 0.;
1113  drEfrac = BigRadius;
1114  Efrac_old = 0.;
1115  Efrac_new = 0.;
1116  }
1117 
1118  /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
1119  /* do not let thickness get too large
1120  * 1 '' radius_next keys from relative depth'',e11.3)') */
1121  if( nzone > 20 )
1122  {
1123  /*drDepth = radius.depth/20.;*/
1124  /* >>chng 02 nov 05, change from 1/20 to 1/10 wasted zones early on */
1125  drDepth = radius.depth/10.;
1126  }
1127  else
1128  {
1129  drDepth = BigRadius;
1130  }
1131 
1132  /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
1133  /* case where stopping thickness or edge specified, need to approach slowly */
1134  thickness_total = BigRadius;
1135  DepthToGo = BigRadius;
1136  if( StopCalc.HColStop < 5e29 )
1137  {
1138  double coleff = SDIV( dense.gas_phase[ipHYDROGEN]*geometry.FillFac);
1139  DepthToGo = MIN2(DepthToGo ,
1140  (StopCalc.HColStop-colden.colden[ipCOL_HTOT]) / coleff );
1141  /* >>chng 03 oct 28, forgot to div col den above by eff density */
1142  thickness_total = MIN2(thickness_total , StopCalc.HColStop / coleff );
1143  }
1144 
1145  if( StopCalc.colpls < 5e29 )
1146  {
1147  double coleff = (double)SDIV( (dense.xIonDense[ipHYDROGEN][1])*geometry.FillFac);
1148  DepthToGo = MIN2(DepthToGo ,
1149  (StopCalc.colpls-colden.colden[ipCOL_Hp]) / coleff );
1150  thickness_total = MIN2(thickness_total , StopCalc.colpls / coleff );
1151  }
1152 
1153  if( StopCalc.col_h2 < 5e29 )
1154  {
1155  /* >>chng 03 apr 15, add this molecular hydrogen */
1156  double coleff = (double)SDIV( hmi.H2_total*geometry.FillFac);
1157  DepthToGo = MIN2(DepthToGo ,
1159  thickness_total = MIN2(thickness_total , StopCalc.col_h2 / coleff );
1160  }
1161 
1162  if( StopCalc.col_h2_nut < 5e29 )
1163  {
1164  /* >>chng 03 apr 15, add this molecular hydrogen */
1165  double coleff = (double)SDIV( (2*hmi.H2_total+dense.xIonDense[ipHYDROGEN][1])*geometry.FillFac);
1166  DepthToGo = MIN2(DepthToGo ,
1168  thickness_total = MIN2(thickness_total , StopCalc.col_h2_nut / coleff );
1169  }
1170 
1171  if( StopCalc.col_H0_ov_Tspin < 5e29 )
1172  {
1173  /* >>chng 05 jan 09, add n(H0)/Tspin */
1174  double coleff = (double)SDIV( dense.xIonDense[ipHYDROGEN][0] / hyperfine.Tspin21cm*geometry.FillFac );
1175  DepthToGo = MIN2(DepthToGo ,
1176  (StopCalc.col_H0_ov_Tspin - colden.H0_ov_Tspin) / coleff );
1177  thickness_total = MIN2(thickness_total , StopCalc.col_H0_ov_Tspin / coleff );
1178  }
1179 
1180  if( StopCalc.col_monoxco < 5e29 )
1181  {
1182  /* >>chng 03 apr 15, add this, CO */
1183  double coleff = (double)SDIV( (findspecies("CO")->hevmol)*geometry.FillFac);
1184  DepthToGo = MIN2(DepthToGo ,
1185  (StopCalc.col_monoxco-findspecies("CO")->hevcol) / coleff );
1186  thickness_total = MIN2(thickness_total , StopCalc.col_monoxco / coleff );
1187  }
1188 
1189  if( StopCalc.colnut < 5e29 )
1190  {
1191  double coleff = (double)SDIV( (dense.xIonDense[ipHYDROGEN][0])*geometry.FillFac);
1192  DepthToGo = MIN2(DepthToGo ,
1193  (StopCalc.colnut - colden.colden[ipCOL_H0]) / coleff );
1194  thickness_total = MIN2(thickness_total , StopCalc.colnut / coleff );
1195  }
1196 
1197  /* this is case where outer radius is set */
1198  if( radius.router[iteration-1] < 5e29 )
1199  {
1200  thickness_total = MIN2(thickness_total , radius.router[iteration-1] );
1201  DepthToGo = MIN2(DepthToGo ,
1203  }
1204 
1205  /* this is case where stopping optical depth was specified */
1206  if( StopCalc.iptnu != rfield.nupper )
1207  {
1208  /* end optical depth has been specified */
1210  DepthToGo = MIN2(DepthToGo ,
1212  }
1213  /* stop AV - usually this is dust, but we consider all opacity sources,
1214  * so always include this */
1215  /* compute some average grain properties */
1216  AV_to_go = BigRadius;
1218  {
1219  /* by default stop av is very large, and opacity can be very small, so ratio
1220  * goes to inf - work with logs to see how big the number is */
1221  double ave = log10(StopCalc.AV_extended - rfield.extin_mag_V_extended) -
1222  log10(rfield.opac_mag_V_extended);
1223  double avp = log10(StopCalc.AV_point - rfield.extin_mag_V_point) -
1224  log10(rfield.opac_mag_V_point);
1225  AV_to_go = MIN2( ave , avp );
1226  if( AV_to_go < 37. )
1227  {
1228  AV_to_go = pow(10., AV_to_go );
1229  /* this is to make sure that we go slightly deeper than AV so that
1230  * we ttrigger this stop */
1231  AV_to_go *= 1.0001;
1232  }
1233  else
1234  AV_to_go = BigRadius;
1235  /*fprintf(ioQQQ,"DEBUG next dr %.3e %.3e %.3e\n", AV_to_go , ave, avp );*/
1236  }
1237 
1238  /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
1239  /* set dr if one of above tests have triggered */
1240  if( DepthToGo <= 0. )
1241  {
1242  TotalInsanity();
1243  }
1244  else if( DepthToGo < BigRadius )
1245  {
1246  /* want to approach the outer edge slowly,
1247  * the need for this logic is most evident in brems.in -
1248  * HI fraction varies across coronal model */
1249  drThickness = MIN2( thickness_total/10. , DepthToGo );
1250  }
1251  else
1252  {
1253  drThickness = BigRadius;
1254  }
1255 
1256  /*fprintf(ioQQQ,
1257  "DEBUG depth2go z%li drThickness2 = %e %e %e\n",
1258  nzone , drThickness , thickness_total , DepthToGo );*/
1259 
1260  /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
1261  /* spherical models, do not want delta R/R big */
1262  drSphere = radius.Radius*0.04;
1263 
1264  /* optical depth to electron scattering */
1265  /* >>chng 04 jun 16, add filling factor, was missing */
1266  dRTaue = radius.drChange/(dense.eden*6.65e-25*geometry.FillFac);
1267  /* >>chng 02 oct 06, increase dr when constant temperature,
1268  * to prevent some ct models from taking too many cells */
1270  dRTaue *= 3.;
1271 
1272  /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
1273  if( dense.flong != 0. )
1274  {
1275  drFluc = 0.628/2./dense.flong;
1276  /* >>chng 04 sep 18, caused cautions that ionization jumps occurred.
1277  * set to have the value */
1278  drFluc /= 2.;
1279  }
1280  else
1281  {
1282  drFluc = BigRadius;
1283  }
1284 
1285  /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
1286  /* if density fluctuations in place then override change in heat
1287  * for dr set */
1288  if( strcmp(dense.chDenseLaw,"SINE") == 0 && dense.lgDenFlucOn )
1289  {
1290  drdHdStep = BigRadius;
1291  }
1292 
1293  /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
1294  /*active dr sets */
1295  /* we are deep into model, use logic since we have several zones
1296  * of old data */
1297  radius.drNext = MIN2( MIN4( drmax, drInter, drLineHeat, winddr ),
1298  MIN4( drFluc, GlobDr, DrGrainHeat, dr_change_heavy ) );
1299  radius.drNext = MIN2( MIN4( radius.drNext, SpecDr, drFail, WindAccelDR ),
1300  MIN3( drSphere, radius.sdrmax, dRTaue ) );
1301  radius.drNext = MIN3( MIN3( radius.drNext, drDest, drdTdStep ),
1302  MIN3( drdHdStep, drConPres, drTab ),
1303  MIN3( drSolomon_BigH2, drLeiden_hack, recom_dr_last_iter ) );
1304  radius.drNext = MIN3( MIN4( radius.drNext, drHion, drDepth, dr_mole_abund ),
1305  MIN4( AV_to_go, drEfrac, drHMase, drThickness ),
1306  MIN3( drHe1ovHe2, drH2_heat_cool, drH2_abund ) );
1307 
1308  /*fprintf(ioQQQ,
1309  "DEBUG depth2go drNext = %e \n", radius.drNext );*/
1310 
1311  /* keep dr constant in first two zones, if it wants to increase,
1312  * but always allow it to decrease.
1313  * to guard against large increases in efrac for balmer cont photo dominated models,
1314  */
1315  if( nzone <= 1 && radius.drNext > OldDR )
1316  {
1317  radius.drNext = OldDR;
1318  }
1319 
1320  /* option to force min drad */
1321  if( radius.drNext < radius.sdrmin )
1322  {
1324  }
1325 
1328  /* a pressure failure has occurred - keep zone the same time, hoping to pass through
1329  * troubled region */
1330  if( !conv.lgConvPres || !conv.lgConvTemp )
1331  {
1333  }
1334 
1335  /* >>chng 05 aug 05, in case of thermal front, where temp is falling quickly and
1336  * conditions change very fast, the zone thickness does not really affect the change
1337  * in conditions and can cause zoning to become very very thin, which causes
1338  * an abort. this occurs between 200 and 1000K. if we are doing temp soln,
1339  * temp is between these values, and temp is changing rapidly, do not make sone
1340  * thickness much smaller */
1341  drThermalFront = BigRadius;
1342  /* do not use thermal front logic in case of recombination time dep cloud */
1343  if( nzone >=5 && !dynamics.lgStatic )
1344  {
1345  /* >>chng 05 aug 23, upper bound of thermal from from 1000K to 4000K */
1346  /*if( phycon.te > 200. && phycon.te < 1000. && */
1347  if( phycon.te > 200. && phycon.te < 3000. &&
1348  /* >>chng 05 aug 23, from > 10% in zone to to two zones > 5%,
1349  * to fix leiden v3 with large H2 */
1350  (struc.testr[nzone-3] - phycon.te)/phycon.te > 0.02 &&
1351  (struc.testr[nzone-4] - phycon.te)/phycon.te > 0.02 &&
1352  (struc.testr[nzone-5] - phycon.te)/phycon.te > 0.02 )
1353  {
1354  /* the 0.91 is to make dr unique, so that print statement that
1355  * follows will identify this reason */
1356  drThermalFront = radius.drad * 0.91;
1357  radius.drNext = drThermalFront;
1358  }
1359  }
1360 
1361  /* dr = zero is a logical mistake */
1362  if( radius.drNext <= 0. )
1363  {
1364  fprintf( ioQQQ, " radius_next finds insane drNext:%10.2e\n",
1365  radius.drNext );
1366  fprintf( ioQQQ, " all drs follow:%10.2e%10.2e%10.2e%10.2e\n all drs follow:%10.2e%10.2e%10.2e%10.2e%10.2e\n all drs follow:%10.2e%10.2e%10.2e%10.2e\n",
1367  drmax, drInter, drLineHeat, winddr, drFluc, GlobDr, SpecDr,
1368  drFail, drSphere, radius.sdrmax, dRTaue,
1369  OldH2Abund, drDepth );
1370  cdEXIT(EXIT_FAILURE);
1371  }
1372 
1373  /* set flag if dr set by maser */
1374  if( fp_equal( radius.drNext, drHMase ) )
1375  {
1376  rt.lgMaserSetDR = true;
1377  }
1378 
1379  /* all this is to only punch on last iteration
1380  * the punch dr command is not really a punch command, making this necessary
1381  * lgDRon is set true if "punch dr" entered */
1382  if( punch.lgDROn )
1383  {
1384  if( punch.lgDRPLst )
1385  {
1386  /* lgDRPLst was set true if "punch" had "last" on it */
1387  if( iterations.lgLastIt )
1388  {
1389  lgDoPun = true;
1390  }
1391  else
1392  {
1393  lgDoPun = false;
1394  }
1395  }
1396  else
1397  {
1398  lgDoPun = true;
1399  }
1400  }
1401  else
1402  {
1403  lgDoPun = false;
1404  }
1405  if( (trace.lgTrace && trace.lgDrBug) || lgDoPun )
1406  {
1407  if( !conv.lgConvTemp && nzone > 0 )
1408  {
1409  fprintf( punch.ipDRout, "#>>>> A temperature failure occured.\n" );
1410  }
1411  if( !conv.lgConvPres && nzone > 0 )
1412  {
1413  fprintf( punch.ipDRout, "#>>>> A pressure failure occured.\n" );
1414  }
1415 
1416  /* this is common part of each line, the zone count, depth, chosen dr, and depth2go */
1417  fprintf( punch.ipDRout , "%ld\t%.5e\t%.3e\t%.3e\t", nzone, radius.depth, radius.drNext, radius.Depth2Go );
1418 
1419  /*=======begin active dr sets */
1420  if( fp_equal( radius.drNext, drLineHeat ) )
1421  {
1422  if( level == 1 )
1423  {
1424  strcpy( chLbl, chLineLbl(&TauLines[ipStrong]) );
1425  fprintf( punch.ipDRout, "level 1 line heating,%10.10s TauIn%10.2e%10.2e%10.2e\n",
1426  chLbl, TauInwd, TauLines[ipStrong].Emis->pump,
1427  TauLines[ipStrong].Emis->Pesc );
1428  }
1429  else if( level == 2 )
1430  {
1431  strcpy( chLbl, chLineLbl(&TauLine2[ipStrong]));
1432  fprintf( punch.ipDRout, "level 2 line heating,%10.10s TauIn%10.2e%10.2e%10.2e\n",
1433  chLbl, TauInwd, TauLine2[ipStrong].Emis->pump,
1434  TauLine2[ipStrong].Emis->Pesc );
1435  }
1436  else
1437  {
1438  fprintf( ioQQQ, " insanity pr line heat\n" );
1439  cdEXIT(EXIT_FAILURE);
1440  }
1441  }
1442 
1443  else if( fp_equal( radius.drNext, drDepth ) )
1444  {
1445  fprintf( punch.ipDRout, "relative depth\n");
1446  }
1447 
1448  else if( fp_equal( radius.drNext, drThermalFront ) )
1449  {
1450  fprintf( punch.ipDRout, "thermal front logic\n");
1451  }
1452 
1453  else if( fp_equal( radius.drNext, dr_change_heavy ) )
1454  {
1455  ASSERT( ichange_heavy_nelem < LIMELM+3 );
1456  fprintf( punch.ipDRout,
1457  "change in ionization, element %s ion (C scale) %li rel change %.2e ion frac %.2e\n",
1458  elementnames.chElementName[ichange_heavy_nelem],
1459  ichange_heavy_ion ,
1460  change_heavy_frac_big ,
1461  frac_change_heavy_frac_big);
1462  }
1463 
1464  else if( fp_equal( radius.drNext, drThickness ) )
1465  {
1466  fprintf( punch.ipDRout, "depth to go\n");
1467  }
1468 
1469  else if( fp_equal( radius.drNext, AV_to_go ) )
1470  {
1471  fprintf( punch.ipDRout, "A_V to go\n");
1472  }
1473 
1474  else if( fp_equal( radius.drNext, drTab ) )
1475  {
1476  fprintf( punch.ipDRout, "spec den law, new old den%10.2e%10.2e\n",
1477  hdnew, dense.gas_phase[ipHYDROGEN] );
1478  }
1479 
1480  else if( fp_equal( radius.drNext, drHMase ) )
1481  {
1482  fprintf( punch.ipDRout, "H maser dTauMase=%10.2e %li %li %li %li\n",
1483  rt.dTauMase,
1484  rt.mas_species,
1485  rt.mas_ion,
1486  rt.mas_hi,
1487  rt.mas_lo );
1488  }
1489 
1490  else if( fp_equal( radius.drNext, drHe1ovHe2 ) )
1491  {
1492  /* radius_next keys from change in He2/He1 ionization, old=%11.3e sv=%11.3e FrLya*/
1493  fprintf( punch.ipDRout, "change in He0/He+ ionization, ratio %.2e\n",
1494  Old_He_atom_ov_ion );
1495  }
1496 
1497  else if( fp_equal( radius.drNext, drH2_heat_cool ) )
1498  {
1499  /* radius_next keys from change in H2 heating, old=%11.3e sv=%11.3e FrLya*/
1500  fprintf( punch.ipDRout, "change in H2 heating/cooling, d(c,h)/H %.2e\n",
1501  dH2_heat_cool );
1502  }
1503 
1504  else if( fp_equal( radius.drNext, drH2_abund ) )
1505  {
1506  /* radius_next keys from change in H2 abundance, old=%11.3e sv=%11.3e FrLya*/
1507  fprintf( punch.ipDRout, "change in H2 abundance, d(H2)/H %.2e\n",
1508  dH2_abund );
1509  }
1510 
1511  else if( fp_equal( radius.drNext, dr_mole_abund ) )
1512  {
1513  /* radius_next keys from change in CO mole abundance */
1514  fprintf( punch.ipDRout, "change in heav ele mole abundance, d(mole)/elem %.2e mole=%i=%s\n",
1515  dCO_abund , mole_dr_change , COmole[mole_dr_change]->label);
1516  }
1517 
1518  else if( fp_equal( radius.drNext, drSolomon_BigH2 ) )
1519  {
1520  /* radius_next keys from change in H2 abundance, old=%11.3e sv=%11.3e FrLya*/
1521  fprintf( punch.ipDRout, "change in big H2 Solomon rate line opt depth\n");
1522  }
1523 
1524  else if( fp_equal( radius.drNext, recom_dr_last_iter ) )
1525  {
1526  /* radius_next keys from previous iteration recom logic */
1527  fprintf( punch.ipDRout, "previous iteration recom logic\n");
1528  }
1529 
1530  else if( fp_equal( radius.drNext, drHion ) )
1531  {
1532  fprintf( punch.ipDRout, "change in H ionization fm to%11.3e%11.3e\n",
1533  SaveOHIIoHI, OHIIoHI );
1534  }
1535 
1536  else if( fp_equal( radius.drNext, drEfrac ) )
1537  {
1538  fprintf( punch.ipDRout,
1539  "change in elec den, rel chng:%11.3e, cur %11.3e old=%11.3e\n",
1540  dEfrac , Efrac_old , Efrac_new );
1541  }
1542 
1543  else if( fp_equal( radius.drNext, drdHdStep ) )
1544  {
1545  fprintf( punch.ipDRout,
1546  "change in heating; current %10.3e delta=%10.3e\n",
1547  thermal.htot, dHdStep );
1548  }
1549 
1550  else if( fp_equal( radius.drNext, drLeiden_hack ) )
1551  {
1552  fprintf( punch.ipDRout,
1553  "Leiden hack\n" );
1554  }
1555 
1556  else if( fp_equal( radius.drNext, drConPres ) )
1557  {
1558  fprintf( punch.ipDRout, "change in con accel\n" );
1559  }
1560 
1561  else if( fp_equal( radius.drNext, drdTdStep ) )
1562  {
1563  fprintf( punch.ipDRout,
1564  "change in temperature; current= %10.3e, dT/T= %.3f\n",
1565  phycon.te, dTdStep );
1566  }
1567 
1568  else if( fp_equal( radius.drNext, radius.sdrmin ) )
1569  {
1570  fprintf( punch.ipDRout, "sdrmin\n" );
1571  }
1572 
1573  else if( fp_equal( radius.drNext, radius.sdrmax ) )
1574  {
1575  fprintf( punch.ipDRout, "sdrmax\n" );
1576  }
1577 
1578  else if( fp_equal( radius.drNext, drSphere ) )
1579  {
1580  fprintf( punch.ipDRout, "sphericity\n" );
1581  }
1582 
1583  else if( fp_equal( radius.drNext, dRTaue ) )
1584  {
1585  fprintf( punch.ipDRout,
1586  "optical depth to electron scattering\n" );
1587  }
1588 
1589  else if( fp_equal( radius.drNext, drFail ) )
1590  {
1591  fprintf( punch.ipDRout,
1592  "temperature failure.\n" );
1593  }
1594 
1595  else if( fp_equal( radius.drNext, drmax ) )
1596  {
1597  fprintf( punch.ipDRout,
1598  "DRMAX; nu opc dr=%10.2e%10.2e%10.2e\n",
1599  freqm, opacm, radius.drChange/
1600  SDIV(opacm) );
1601  }
1602 
1603  else if( fp_equal( radius.drNext, drInter ) )
1604  {
1605  fprintf( punch.ipDRout,
1606  "cont inter nu=%10.2e opac=%10.2e\n",
1607  freqm, opacm );
1608  }
1609 
1610  else if( fp_equal( radius.drNext, DrGrainHeat ) )
1611  {
1612 
1613  fprintf( punch.ipDRout,
1614  "grain heating nu=%10.2e opac=%10.2e\n",
1615  grfreqm, gropacm );
1616  }
1617 
1618  else if( fp_equal( radius.drNext, winddr ) )
1619  {
1620  fprintf( punch.ipDRout,
1621  "Wind, dVelRelative=%.3e windv=%.3e\n",
1622  dVelRelative ,
1623  wind.windv);
1624  }
1625 
1626  else if( fp_equal( radius.drNext, drFluc ) )
1627  {
1628  fprintf( punch.ipDRout,
1629  "density fluctuations\n" );
1630  }
1631 
1632  else if( fp_equal( radius.drNext, GlobDr ) )
1633  {
1634  fprintf( punch.ipDRout,
1635  "GLOB law new dr=%10.2e HDEN=%10.2e\n",
1636  GlobDr,
1638  }
1639 
1640  else if( fp_equal( radius.drNext, SpecDr ) )
1641  {
1642  fprintf( punch.ipDRout,
1643  "special law new dr=%10.2e HDEN=%10.2e\n",
1644  SpecDr,
1646  }
1647 
1648  else if( fp_equal( radius.drNext, OldDR ) )
1649  {
1650  fprintf( punch.ipDRout, "old DR.\n" );
1651  }
1652 
1653  else
1654  {
1655  fprintf( punch.ipDRout,
1656  " %4ld radius_next keys from insanity %10.2e\n",
1657  nzone, radius.drNext );
1658 
1659  fprintf( ioQQQ,
1660  " %4ld radius_next keys from insanity %10.2e\n",
1661  nzone, radius.drNext );
1662  cdEXIT(EXIT_FAILURE);
1663  }
1664 
1665  /*======end active dr sets */
1666  }
1667 
1668  /* this is general code that prevents zone thickness drNext from
1669  * becoming too thin, something that can happen for various bad reasons
1670  * HOWEVER we do not want to do this test for certain density laws,
1671  * for which very small zone thicknesses are unavoidable
1672  * the special cases are:
1673  * special density law,
1674  * globule density law,
1675  * carbon +-0 i front
1676  * the fluctuations command
1677  * drMinimum was set in radius_first to either sdrmin (set drmin) or
1678  * some fraction of the initial radius - it is always set
1679  * to something non-trivial.
1680  * sdrmin is set with the "set dr" command - is SMALLFLOAT by default */
1681  if( ((strcmp(dense.chDenseLaw,"DLW1") != 0 &&
1682  strcmp(dense.chDenseLaw ,"GLOB") != 0) )&&
1683  /* >>chng 04 feb 19, do not use this test - errors can still happen
1684  * when all C is atomic! */
1685  /*dense.xIonDense[ipCARBON][0]/dense.gas_phase[ipCARBON] < 0.05) && */
1686  (dense.flong == 0.) &&
1687  /* >>chng 01 aug 11, add check against stopping on depth to go */
1688  radius.drNext != DepthToGo )
1689  {
1690  /* don't let dr get smaller than drMinimum, if this resets drNext
1691  * then code stops with warning that zones got too thin
1692  * this can happen due to numerical oscillations, which the code
1693  * tries to damp out by making the zones thinner.
1694  * scale by radius.Radius/radius.rinner to stop very spherical
1695  * simulations from false trigger on dr too small */
1697  /* drMinimum is drad * hden, to make it proportional to optical depth.
1698  * This avoids false trigger across thermal fronts. */
1700  radius.drMinimum )
1701  {
1703  /* leaving this at true will cause the model to stop with a warning
1704  * that the zone thickness is too small */
1705  radius.lgDrMinUsed = true;
1706  /* set abort handler */
1707  lgAbort = true;
1708  /* must decrement nzone, since we will not complete this zone, and will not have
1709  * valid structure data for it */
1710  --nzone;
1711  fprintf( ioQQQ,
1712  "\n DISASTER PROBLEM radius_next finds dr too small and aborts. "
1713  "This is zone %ld iteration %ld. The thickness was %.2e\n",
1714  nzone,
1715  iteration,
1716  radius.drNext);
1717  fprintf( ioQQQ,
1718  " If this simulation seems reasonable you can override this limit "
1719  "with the command SET DRMIN %.2e\n\n",
1720  radius.drNext*2);
1721  /* set abort flag */
1722  lgAbort = true;
1723  return(1);
1724  }
1725  }
1726 
1727  /* factor to allow for slop in floating numbers */
1728 # define Z 1.0001
1729 
1730  /* following is to make thickness of model exact
1731  * n.b., on last zone, drNext can be NEGATIVE!!
1732  * DEPTH was incremented at start of zone calc in newrad,
1733  * has been outer edge of zone all throughout */
1735  radius.depth)*Z);
1736 
1737  /* this means outer limit exceeded */
1738  if( radius.drNext < 0. )
1739  {
1740  radius.lgDrNeg = true;
1741  }
1742  else
1743  {
1744  radius.lgDrNeg = false;
1745  }
1746 
1747  ASSERT( radius.drNext > 0. );
1748 
1749  if( trace.lgTrace )
1750  {
1751  fprintf( ioQQQ, " radius_next chooses next drad drNext=%.4e; this drad was%12.4e\n",
1752  radius.drNext, radius.drad );
1753  }
1754 
1755  return( 0 );
1756 }
1757 
1758 /*ContRate called by radius_next to find energy of maximum continuum-gas interaction */
1759 STATIC void ContRate(double *freqm,
1760  double *opacm)
1761 {
1762  long int i,
1763  ipHi,
1764  iplow,
1765  limit;
1766  double FreqH,
1767  Freq_nonIonizing,
1768  Opac_Hion,
1769  Opac_nonIonizing,
1770  Rate_max_Hion,
1771  Rate_max_nonIonizing;
1772 
1773  DEBUG_ENTRY( "ContRate()" );
1774 
1775  /*
1776  * find maximum continuum interaction rate,
1777  * these should be reset in following logic without exception,
1778  * if they are still zero at the end we have a logical error
1779  */
1780  Rate_max_nonIonizing = 0.;
1781  Freq_nonIonizing = 0.;
1782  Opac_nonIonizing = 0.;
1783 
1784  /* this must be reset to val >= 0 */
1785  *opacm = -1.;
1786  *freqm = -1.;
1787 
1788  /* do up to carbon photo edge if carbon is turned on */
1789  /* >>>chng 00 apr 07, add test for whether element is turned on */
1790  if( dense.lgElmtOn[ipCARBON] )
1791  {
1792  /* carbon is turned on, use carbon 1 edge */
1793  ipHi = Heavy.ipHeavy[ipCARBON][0] - 1;
1794  }
1795  else
1796  {
1797  /* carbon turned off, use hydrogen Balmer continuum */
1799  }
1800 
1801  /* start at plasma frequency */
1802  for( i=rfield.ipPlasma; i < ipHi; i++ )
1803  {
1804  /* this does not have grain opacity since grains totally passive
1805  * at energies smaller than CI edge */
1806  if( rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*(opac.opacity_abs[i] -
1807  gv.dstab[i]*dense.gas_phase[ipHYDROGEN]) > Rate_max_nonIonizing )
1808  {
1809  Rate_max_nonIonizing = rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*
1811  Freq_nonIonizing = rfield.anu[i];
1812  Opac_nonIonizing = opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN];
1813  }
1814  }
1815 
1816  /* not every continuum extends beyond C edge-this whole loop can add to zero
1817  * use total opacity here
1818  * test is to put in fir continuum if free free heating is important */
1819  if( CoolHeavy.brems_heat_total/thermal.htot > 0.05 )
1820  {
1821  /* this is index for energy where cloud free free optical depth is unity,
1822  * is zero if no freq are opt thin */
1823  iplow = MAX2(1 , rfield.ipEnergyBremsThin );
1824  }
1825  else
1826  {
1827  /* >>>chng 00 apr 07, from Heavy.ipHeavy[0][5] to ipHi defined above, since
1828  * would crash if element not defined */
1829  iplow = ipHi;
1830  }
1831  /* do not go below plasma frequency */
1832  iplow = MAX2( rfield.ipPlasma , iplow );
1833 
1834  /* this energy range from carbon edge to hydrogen edge */
1836  for( i=iplow; i < limit; i++ )
1837  {
1838  if( rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*(opac.opacity_abs[i] -
1839  gv.dstab[i]*dense.gas_phase[ipHYDROGEN]) > Rate_max_nonIonizing )
1840  {
1841  Rate_max_nonIonizing = rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*
1843  Freq_nonIonizing = rfield.anu[i];
1844  Opac_nonIonizing = opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN];
1845  }
1846  }
1847 
1848  /* variables to check continuum interactions over Lyman continuum */
1849  Rate_max_Hion = 0.;
1850  Opac_Hion = 0.;
1851  FreqH = 0.;
1852 
1853  /* not every continuum extends beyond 1 Ryd-this whole loop can add to zero */
1854  for( i=iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1; i < rfield.nflux; i++ )
1855  {
1856  if( rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*(opac.opacity_abs[i] -
1857  gv.dstab[i]*dense.gas_phase[ipHYDROGEN]) > Rate_max_Hion )
1858  {
1859  /* Rate_max_Hion = anu(i)*flux(i)/widflx(i)*opac(i) */
1860  Rate_max_Hion = rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*
1862  FreqH = rfield.anu[i];
1863  Opac_Hion = opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN];
1864  }
1865  }
1866 
1867  /* use Lyman continuum if its opacity is larger than non-h ion */
1868  if( Rate_max_nonIonizing < 1e-30 && Opac_Hion > SMALLFLOAT )
1869  {
1870  /* this happens for laser source - use Lyman continuum */
1871  *opacm = Opac_Hion;
1872  *freqm = FreqH;
1873  }
1874  /* >>chng 05 aug 03, add last test on Opac_Hion for case where we go very
1875  * deep and very little radiation is left */
1876  else if( Opac_Hion > Opac_nonIonizing && Rate_max_Hion/Rate_max_nonIonizing > 1e-10 && Opac_Hion > SMALLFLOAT )
1877  {
1878  /* use Lyman continuum */
1879  *opacm = Opac_Hion;
1880  *freqm = FreqH;
1881  }
1882  else
1883  {
1884  /* not much rate in the Lyman continuum, stick with low energy */
1885  *opacm = Opac_nonIonizing;
1886  *freqm = Freq_nonIonizing;
1887  }
1888 
1889 # if 0
1890  /*>>chng 06 aug 12, do not let zones become optically thick to
1891  * peak of dust emission - one of NA's vastly optically thick dust clouds
1892  * did not conserve total energy - very opticallly thick to ir dust emission
1893  * so ir was absorbed and reemitted many times
1894  * this makes sure the cells remain optically thin to dust emission */
1895  if( gv.lgDustOn && gv.lgGrainPhysicsOn )
1896  {
1897  double fluxGrainPeak = -1.;
1898  long int ipGrainPeak = -1;
1899  long int ipGrainPeak2 = -1;
1900 
1901  i = 0;
1902  while( rfield.anu[i] < 0.0912 && i < (rfield.nflux-2) )
1903  {
1904  /* >>chng 06 aug 23. Only want to do this test for the IR dust continuum, therefore only
1905  * check on optical depth for wavelengths greater than 1 micron */
1906  if( gv.GrainEmission[i]*rfield.anu[i]*opac.opacity_abs[i] > fluxGrainPeak )
1907  {
1908  ipGrainPeak = i;
1909  fluxGrainPeak = gv.GrainEmission[i]*rfield.anu[i]*opac.opacity_abs[i];
1910  }
1911  ++i;
1912  }
1913  ASSERT( fluxGrainPeak>=0. && ipGrainPeak >=0 );
1914 
1915  /* >>chng 06 aug 23. We have just found the wavelength and flux at the peak of the IR emission.
1916  * Now we want to find the wavelength, short of the peak, which corresponds to 5% of the
1917  * peak emission. This wavelength will be where the code checks to make sure the zone has
1918  * not become optically thick. Since dust opacity generally decreases with wavelength, this assures that
1919  * the optical depths are small for all wavelengths where the flux is appreciable. Tests show that
1920  * this allows for flux/luminosity conservation to within ~5% for an AV of 1e4 mag and at least 2 iterations*/
1921  i = ipGrainPeak;
1922  /* >>chng 06 aug 23. Only want to do this test for the IR dust continuum, therefore only
1923  * check on opacities for wavelengths shortward of the peak and greater than 1 micron
1924  * this routine can be called with the dust emission totally zero - it is only evaluated
1925  * late to save time. don't do the following tests if peak dust emission is zero */
1926  if( fluxGrainPeak > 0. )
1927  {
1928  while( rfield.anu[i] < 0.0912 && i < (rfield.nflux-2) )
1929  {
1930  /* find wavelength where flux = 5% of peak flux, shortward of the peak */
1931  if( gv.GrainEmission[i]*rfield.anu[i]*opac.opacity_abs[i] > 0.05*fluxGrainPeak )
1932  {
1933  /* This will be the array number and flux at 10% of the peak */
1934  ipGrainPeak2 = i;
1935  }
1936  ++i;
1937  }
1938  ASSERT( ipGrainPeak2>=0 );
1939 
1940  /* use this as a limit to the dr if opacity is larger */
1941  if( opac.opacity_abs[ipGrainPeak2] > *opacm )
1942  {
1943  /* scale opacity by factor of 10, which further decreases the zone thickness*/
1944  *opacm = opac.opacity_abs[ipGrainPeak2]*10.;
1945  *freqm = rfield.anu[ipGrainPeak2];
1946  /*fprintf(ioQQQ,"DEBUT opac peak %.2e %.2e \n",
1947  *opacm , *freqm );*/
1948  }
1949  }
1950  }
1951 # endif
1952 
1953  {
1954  /* following should be set true to print contributors */
1955  enum {DEBUG_LOC=false};
1956  if( DEBUG_LOC )
1957  {
1958  fprintf(ioQQQ,"conratedebug \t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
1959  Rate_max_nonIonizing,Freq_nonIonizing,Opac_nonIonizing,
1960  Rate_max_Hion,FreqH ,Opac_Hion,*freqm,*opacm
1961  );
1962 
1963  }
1964  }
1965 
1966  /* these were set to -1 at start, and must have been reset in one of the
1967  * two loops. Logic error if still <0. */
1968  /* >>chng 05 aug 03, change logic to -1 on entry and check at least zero
1969  * here - will be zero if NO radiation field exists, perhaps since totally
1970  * attenuated */
1971  ASSERT( *opacm >= 0. && *freqm >= 0. );
1972  return;
1973 }
1974 
1975 /*GrainRateDr called by radius_next to find grain heating rate dr */
1976 STATIC void GrainRateDr(double *grfreqm,
1977  double *gropacm)
1978 {
1979  long int i,
1980  iplow;
1981  double xMax;
1982 
1983  DEBUG_ENTRY( "GrainRateDr()" );
1984 
1985  /* in all following changed from anu2 to anu july 25 95
1986  *
1987  * find maximum continuum interaction rate */
1988 
1989  /* not every continuum extends beyond C edge-this whole loop can add to zero
1990  * use total opacity here
1991  * test is to put in fir continuum if free free heating is important */
1992  if( CoolHeavy.brems_heat_total/thermal.htot > 0.05 )
1993  {
1994  /* this is pointer to energy where cloud free free optical depth is unity,
1995  * is zero if no freq are opt thin */
1996  iplow = MAX2(1 , rfield.ipEnergyBremsThin );
1997  }
1998  else
1999  {
2000  /* do up to carbon photo edge if carbon is turned on */
2001  /* >>>chng 00 apr 07, add test for whether element is turned on */
2002  if( dense.lgElmtOn[ipCARBON] )
2003  {
2004  /* carbon is turned on, use carbon 1 edge */
2005  iplow = Heavy.ipHeavy[ipCARBON][0];
2006  }
2007  else
2008  {
2009  /* carbon truned off, use hydrogen balmer continuum */
2011  }
2012 
2013  }
2014 
2015  xMax = -1.;
2016  /* integrate up to H edge */
2017  for( i=iplow-1; i < Heavy.ipHeavy[ipHYDROGEN][0]; i++ )
2018  {
2019  if( rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*opac.opacity_abs[i] > xMax )
2020  {
2021  xMax = rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*
2022  opac.opacity_abs[i];
2023  *grfreqm = rfield.anu[i];
2024  *gropacm = opac.opacity_abs[i];
2025  }
2026  }
2027  /* integrate up to heii edge if he is turned on,
2028  * this logic will not make sense if grains on but he off, which in itself makes no sense*/
2029  if( dense.lgElmtOn[ipHELIUM] )
2030  {
2031  for( i=Heavy.ipHeavy[ipHYDROGEN][0]-1; i < Heavy.ipHeavy[ipHELIUM][1]; i++ )
2032  {
2033  if( rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*opac.opacity_abs[i] > xMax )
2034  {
2035  xMax = rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*
2036  opac.opacity_abs[i];
2037  *grfreqm = rfield.anu[i];
2038  *gropacm = opac.opacity_abs[i];
2039  }
2040  }
2041  }
2042 
2043  /* possible that there is NO ionizing radiation, in extreme cases,
2044  * if so then xMax will still be negative */
2045  if( xMax <= 0. )
2046  {
2047  *gropacm = 0.;
2048  *grfreqm = 0.;
2049  }
2050  return;
2051 }
2052 
2053 #undef WIND_CHNG_VELOCITY_RELATIVE
2054 
2055 #if 0
2056 /*ChkRate called by radius_next to check how rates of destruction of various species changes */
2057 STATIC void ChkRate(
2058  /* element number on C scale */
2059  long int nelem,
2060  /* change in destruction rate */
2061  double *dDestRate,
2062  /* old and new destruction rates */
2063  double *DestRateOld,
2064  double *DestRateNew,
2065  /* stage of ionization on the physical scale */
2066  long int *istage)
2067 {
2068  long int i;
2069 
2070  double average,
2071  dDest;
2072 
2073  static double OldDest[LIMELM][LIMELM];
2074 
2075  DEBUG_ENTRY( "ChkRate()" );
2076 
2077  /* return if this element is not turned on */
2078  if( !dense.lgElmtOn[nelem] )
2079  {
2080  *dDestRate = 1e-3;
2081  *DestRateOld = 1e-3;
2082  *DestRateNew = 1e-3;
2083  *istage = 0;
2084  return;
2085  }
2086 
2087  /* for first zone, and during search, we will do nothing but
2088  * still must return finite numbers */
2089  *istage = 1;
2090  *dDestRate = 0.;
2091  *DestRateOld = 0.;
2092  *DestRateNew = 0.;
2093  *dDestRate = 0.;
2094 
2095  if( nzone <= 1 )
2096  {
2097  for( i=0; i < nelem+1; i++ )
2098  {
2099  OldDest[nelem][i] = ionbal.RateIonizTot[nelem][i];
2100  }
2101  }
2102 
2103  else if( dense.xIonDense[nelem][0]/dense.gas_phase[nelem] < 0.9 )
2104  {
2105  /* do not use this method if everything is atomic */
2106  for( i=0; i < (nelem); i++ )
2107  {
2108  /* last check below, .5 chosen so that do not key off
2109  * predominantly neutral species where self-opacity
2110  * could cause oscillations */
2111  if( ((dense.xIonDense[nelem][i]/dense.gas_phase[nelem] > 0.01 &&
2112  dense.xIonDense[nelem][i]/dense.gas_phase[nelem] < 0.9) &&
2113  dense.xIonDense[nelem][i+1]/dense.gas_phase[nelem] > .05) &&
2114  OldDest[nelem][i] > 0. )
2115  {
2116  /* last check on old dest in case just lowered ionization
2117  * stage, so no history */
2118  /* check that rate is positive */
2119  if( ionbal.RateIonizTot[nelem][i] <= 0. )
2120  {
2121  fprintf( ioQQQ, " ChkRate gets insane destruction rate for ion%4ld%4ld%10.2e\n",
2122  nelem+1, i, ionbal.RateIonizTot[nelem][i] );
2123  cdEXIT(EXIT_FAILURE);
2124  }
2125 
2126  /* do not consider unless of middling ionization, and
2127  * rate is going down (to prevent dr osciallating)
2128  * no absolute value in following since do not want to
2129  * consider cases where ionization rate increases */
2130  average = (OldDest[nelem][i] + ionbal.RateIonizTot[nelem][i])* 0.5;
2131 
2132  dDest = (OldDest[nelem][i] - ionbal.RateIonizTot[nelem][i])/ average;
2133  /* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ term + if rate going down */
2134 
2135  if( *dDestRate < dDest )
2136  {
2137  /* biggest rate so far, remember change in rates and ionization stage */
2138  *dDestRate = dDest;
2139  *istage = i+1;
2140  *DestRateOld = OldDest[nelem][i];
2141  *DestRateNew = ionbal.RateIonizTot[nelem][i];
2142  }
2143 
2144  }
2145  OldDest[nelem][i] = ionbal.RateIonizTot[nelem][i];
2146  }
2147  }
2148  return;
2149 }
2150 #endif

Generated for cloudy by doxygen 1.8.1.2