cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
iter_startend.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 /*IterStart set and save values of many variables at start of iteration */
4 /*IterRestart restart iteration */
5 #include "cddefines.h"
6 #include "cddrive.h"
7 #include "physconst.h"
8 #include "iso.h"
9 #include "grainvar.h"
10 #include "taulines.h"
11 #include "hydrogenic.h"
12 #include "struc.h"
13 #include "dynamics.h"
14 #include "prt.h"
15 #include "hyperfine.h"
16 #include "dense.h"
17 #include "magnetic.h"
18 #include "continuum.h"
19 #include "geometry.h"
20 #include "h2.h"
21 #include "he.h"
22 #include "grains.h"
23 #include "atomfeii.h"
24 #include "pressure.h"
25 #include "stopcalc.h"
26 #include "conv.h"
27 #include "mean.h"
28 #include "ca.h"
29 #include "thermal.h"
30 #include "atoms.h"
31 #include "wind.h"
32 #include "opacity.h"
33 #include "timesc.h"
34 #include "trace.h"
35 #include "colden.h"
36 #include "secondaries.h"
37 #include "hmi.h"
38 #include "radius.h"
39 #include "phycon.h"
40 #include "called.h"
41 #include "mole.h"
42 #include "rfield.h"
43 #include "ionbal.h"
44 #include "atmdat.h"
45 #include "lines.h"
46 #include "molcol.h"
47 #include "input.h"
48 #include "rt.h"
49 #include "iterations.h"
50 /* these are the static saved variables, set in IterStart and reset in IterRestart */
51 
60 
61 /* arguments are atomic number, ionization stage*/
63 static double HeatSave[LIMELM][LIMELM];
65 dndtSave;
66 /* save values of dr and drNext */
67 static double drSave , drNextSave;
68 
69 /* also places to save them between iterations */
70 static long int IonLowSave[LIMELM],
72 
73 /* these are used to reset the level populations of the h and he model atoms */
74 /*static realnum hnsav[LIMELM][LMHLVL+1], */
75 static bool lgHNSAV = false;
76 
77 static realnum ***HOpacRatSav ,
79 
81 
82 static double ortho_save , para_save;
83 static double ***hnsav,
84  edsav;
85 
88 
89 void IterStart(void)
90 {
91  long int i,
92  ion,
93  ipISO,
94  n ,
95  nelem;
96  double fhe1nx,
97  ratio;
98 
99  DEBUG_ENTRY( "IterStart()" );
100 
101  /* allocate two matrices if first time in this core load
102  * these variables are malloced here because they are file static -
103  * used to save values at start of first zone, and reset to this value at
104  * start of new iteration */
105  if( !lgHNSAV )
106  {
107  /* set flag so we never do this again */
108  lgHNSAV = true;
109 
110  HOpacRatSav = (realnum ***)MALLOC(sizeof(realnum **)*NISO );
111 
112  hnsav = (double ***)MALLOC(sizeof(double **)*NISO );
113 
114  for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
115  {
116 
117  HOpacRatSav[ipISO] = (realnum **)MALLOC(sizeof(realnum *)*LIMELM );
118 
119  hnsav[ipISO] = (double **)MALLOC(sizeof(double *)*LIMELM );
120 
121  /* now do the second dimension */
122  for( nelem=ipISO; nelem<LIMELM; ++nelem )
123  {
124  HOpacRatSav[ipISO][nelem] = (realnum *)MALLOC(sizeof(realnum)*(unsigned)(iso.numLevels_max[ipISO][nelem]+1) );
125 
126  hnsav[ipISO][nelem] = (double *)MALLOC(sizeof(double)*(unsigned)(iso.numLevels_max[ipISO][nelem]+1) );
127  }
128  }
129  saveMoleSource =
130  (double**)MALLOC(sizeof(double*)*(unsigned)LIMELM );
131  saveMoleSink =
132  (double**)MALLOC(sizeof(double*)*(unsigned)LIMELM );
134  (realnum***)MALLOC(sizeof(realnum**)*(unsigned)LIMELM );
135  for(nelem=0; nelem<LIMELM; ++nelem )
136  {
137  /* chemistry source and sink terms for ionization ladders */
138  saveMoleSource[nelem] =
139  (double*)MALLOC(sizeof(double)*(unsigned)(nelem+2) );
140  saveMoleSink[nelem] =
141  (double*)MALLOC(sizeof(double)*(unsigned)(nelem+2) );
142  SaveMoleChTrRate[nelem] =
143  (realnum**)MALLOC(sizeof(realnum*)*(unsigned)(nelem+2) );
144  for( ion=0; ion<nelem+2; ++ion )
145  {
146  SaveMoleChTrRate[nelem][ion] =
147  (realnum*)MALLOC(sizeof(realnum)*(unsigned)(nelem+2) );
148  }
149  }
150  }
151 
152  /* IterStart is called at the start of EVERY iteration */
153  if( trace.lgTrace )
154  {
155  fprintf( ioQQQ, " IterStart called.\n" );
156  }
157 
158  /* check if this is the last iteration */
159  if( iteration > iterations.itermx )
160  {
161  iterations.lgLastIt = true;
162  }
163  else
164  {
165  iterations.lgLastIt = false;
166  }
167 
168  /* reset some vars dealing with H2
169  H2_Reset(); */
170 
171  /* flag set true in RT_line_one_tauinc if maser ever capped */
172  rt.lgMaserCapHit = false;
173 
174  /* zero out charge transfer heating and cooling fractions */
175  atmdat.HCharHeatMax = 0.;
176  atmdat.HCharCoolMax = 0.;
177 
178  /* the reason this calculation stops */
179  strncpy( StopCalc.chReasonStop, "reason not specified.",
180  sizeof(StopCalc.chReasonStop) );
181 
182  /* zero fractions of He0 destruction due to 23S */
183  he.nzone = 0;
184  he.frac_he0dest_23S = 0.;
186 
187  /* save expansion rate here */
189 
192  /* remains neg if not evaluated */
195 
196  for( i=0; i < mole.num_comole_calc; i++ )
197  {
198  timesc.AgeCOMoleDest[i] = 0.;
199  }
200  timesc.BigCOMoleForm = 0.;
201  timesc.TimeH21cm = 0.;
202  thermal.CoolHeatMax = 0.;
203  hydro.HCollIonMax = 0.;
204 
205  atmdat.HIonFracMax = 0.;
206 
207  /* will save highest n=2p/1s ratio for hydrogen atom*/
208  hydro.pop2mx = 0.;
209  hydro.lgHiPop2 = false;
210 
211  hydro.nLyaHot = 0;
212  hydro.TLyaMax = 0.;
213 
214  /* evaluate the gas and radiation pressure */
215  /* this sets values of pressure.PresTotlCurr */
216  pressure.PresInteg = 0.;
217  pressure.pinzon = 0.;
218  PresTotCurrent();
219 
220  dynamics.HeatMax = 0.;
221  dynamics.CoolMax = 0.;
222 
223  /* reset these since we now have a valid solution */
224  pressure.pbeta = 0.;
225  pressure.RadBetaMax = 0.;
226  pressure.lgPradCap = false;
227  pressure.lgPradDen = false;
228  /* this flag will say we hit the sonic point */
229  pressure.lgSonicPoint = false;
230  pressure.lgStrongDLimbo = false;
231 
232  /* define two timescales for equilibrium, Compton and thermal */
233  timesc.tcmptn = 1.e0/(rfield.qtot*6.65e-25*dense.eden);
235 
236  /* this will be total mass in computed structure */
237  dense.xMassTotal = 0.;
238 
239  for( nelem=0; nelem < LIMELM; nelem++ )
240  {
241 
242  if( dense.lgElmtOn[nelem] )
243  {
244 
245  /* now zero out the ionic fractions */
246  for( ion=0; ion < (nelem + 2); ion++ )
247  {
248  xIonFsave[nelem][ion] = dense.xIonDense[nelem][ion];
249  }
250 
251  for( ion=0; ion < LIMELM; ion++ )
252  {
253  HeatSave[nelem][ion] = thermal.heating[ion][nelem];
254  }
255  }
256  }
257 
258  /* >>chng 02 jan 28, add ipISO loop */
259  for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
260  {
261  /* >>chng 03 apr 11, from 0 to ipISO */
262  for( nelem=ipISO; nelem < LIMELM; nelem++ )
263  {
264  if( dense.lgElmtOn[nelem] )
265  {
266  for( n=0; n < iso.numLevels_max[ipISO][nelem]; n++ )
267  {
268  HOpacRatSav[ipISO][nelem][n] = iso.ConOpacRatio[ipISO][nelem][n];
269  hnsav[ipISO][nelem][n] = StatesElem[ipISO][nelem][n].Pop;
270  }
271  }
272  iso.TwoNu_induc_dn_max[ipISO][nelem] = 0.;
273  }
274  }
275 
277  rfield.EnergyBremsThin = 0.;
278 
279  atoms.nNegOI = 0;
280  for( i=0; i< N_OI_LEVELS; ++i )
281  {
282  atoms.popoi[i] = 0.;
283  }
284 
285  /* save molecular fractions and ionization range */
286  for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
287  {
288  /* save molecular densities */
289  xMolFracsSave[nelem] = dense.xMolecules[nelem];
290  gas_phase_save[nelem] = dense.gas_phase[nelem];
291  IonLowSave[nelem] = dense.IonLow[nelem];
292  IonHighSave[nelem] = dense.IonHigh[nelem];
293  }
294  /*fprintf(ioQQQ,"DEBUG IterStart save set to gas_phase hden %.4e \n",
295  dense.gas_phase[0]);*/
296 
297  edsav = dense.eden;
299 
300  /* save molecular densities */
301  for( i=0; i<N_H_MOLEC; i++)
302  {
303  hmsav[i] = hmi.Hmolec[i];
304  }
307 
308  for( i=0; i<mole.num_comole_calc; ++i )
309  {
310  COmole[i]->hevmol_save = COmole[i]->hevmol;
311  }
312 
313  for( nelem=ipHYDROGEN; nelem < LIMELM; ++nelem )
314  {
315  /* these have one more ion than above */
316  for( ion=0; ion<nelem+2; ++ion )
317  {
318  long int ion2;
319  /* zero out the source and sink arrays */
320  saveMoleSource[nelem][ion] = mole.source[nelem][ion];
321  saveMoleSink[nelem][ion] = mole.sink[nelem][ion];
322  /*>>chng 06 jun 27, move from here, where leak occurs, to
323  * above where xMoleChTrRate first created */
324  /*mole.xMoleChTrRate[nelem][ion] = (realnum*)MALLOC(sizeof(realnum)*(unsigned)(nelem+2));*/
325  for( ion2=0; ion2<nelem+2; ++ion2 )
326  {
327  SaveMoleChTrRate[nelem][ion][ion2] = mole.xMoleChTrRate[nelem][ion][ion2];
328  }
329  }
330  }
331 
334 
336 
339 
343 
345 
348 
352 
353  /* save zone thickness, and next zone thickness */
354  drSave = radius.drad;
356  /* remember largest and smallest dr over iteration */
359 
361 
362  colden.TotMassColl = 0.;
363  colden.tmas = 0.;
364  colden.wmas = 0.;
365  colden.rjnmin = FLT_MAX;
366  colden.ajmmin = FLT_MAX;
367 
368  thermal.nUnstable = 0;
369  thermal.lgUnstable = false;
370 
375 
376  /* plus 1 is to zero out point where unit integration occurs */
377  for( i=0; i < rfield.nupper+1; i++ )
378  {
379  /* diffuse fields continuum */
380  /* >>chng 03 nov 08, recover SummedDif */
382  rfield.ConEmitReflec[i] = 0.;
383  rfield.ConEmitOut[i] = 0.;
384  rfield.ConInterOut[i] = 0.;
385  /* save OTS continuum for next iteration */
386  rfield.otssav[i][0] = rfield.otscon[i];
387  rfield.otssav[i][1] = rfield.otslin[i];
388  rfield.outlin[i] = 0.;
389  rfield.outlin_noplot[i] = 0.;
392 
393  /* save initial absorption, scattering opacities for next iteration */
396 
397  /* will accumulate optical depths through cloud */
398  opac.TauAbsFace[i] = opac.taumin;
400  /* >>chng 99 dec 04, having exactly 1 zone thickness for first zone caused discontinuity
401  * for heating in very high T models in func_map.in. zone 1 and 2 were 20% different,
402  * since tau in is 1e-20, e2 is 0.9999, and so some H ots
403  * these were not here at all - changed to dr/2 */
404  /* attenuation of flux by optical depths IN THIS ZONE
405  * DirectionalCosin is 1/COS(theta), is usually 1, reset with illuminate command,
406  * option for illumination of slab at an angle */
407  /* >>chng 04 oct 09, from drad to radius.drad_x_fillfac - include fill fac, PvH */
409 
410  /* e(-tau) in inward direction, up to illuminated face */
411  opac.ExpmTau[i] = (realnum)opac.ExpZone[i];
412 
413  /* e2(tau) in inward direction, up to illuminated face, this is used to get the
414  * recombination escape probability */
416  }
417 
419 
420  /* this zeros out arrays to hold mean ionization fractions
421  * later entered by mean
422  * read out by setlim */
423  MeanZero();
424 
425  /* zero out the column densities */
426  for( i=0; i < NCOLD; i++ )
427  {
428  colden.colden[i] = 0.;
429  }
430  colden.He123S = 0.;
431  colden.H0_ov_Tspin = 0.;
432  colden.OH_ov_Tspin = 0.;
433  colden.coldenH2_ov_vel = 0.;
434 
435  /* upper and lower levels of H0 1s */
438 
439  h2.ortho_colden = 0.;
440  h2.para_colden = 0.;
441 
442  for( i=0; i < mole.num_comole_calc; i++ )
443  {
444  COmole[i]->hevcol = 0.;
445  }
446  for( i=0; i < 5; i++ )
447  {
448  colden.C2Pops[i] = 0.;
449  colden.C2Colden[i] = 0.;
450  /* pops and column density for SiII atom */
451  colden.Si2Pops[i] = 0.;
452  colden.Si2Colden[i] = 0.;
453  }
454  for( i=0; i < 3; i++ )
455  {
456  /* pops and column density for CI atom */
457  colden.C1Pops[i] = 0.;
458  colden.C1Colden[i] = 0.;
459  /* pops and column density for OI atom */
460  colden.O1Pops[i] = 0.;
461  colden.O1Colden[i] = 0.;
462  /* pops and column density for CIII atom */
463  colden.C3Pops[i] = 0.;
464  }
465  for( i=0; i < 4; i++ )
466  {
467  /* pops and column density for CIII atom */
468  colden.C3Colden[i] = 0.;
469  }
470 
471  /* these are some line of sight emission measures */
472  colden.dlnenp = 0.;
473  colden.dlnenHep = 0.;
474  colden.dlnenHepp = 0.;
475  colden.dlnenCp = 0.;
476  colden.H0_ov_Tspin = 0.;
477  colden.OH_ov_Tspin = 0.;
478 
479  /* now zero heavy element molecules */
480  molcol("ZERO",ioQQQ);
481 
482  /* this will be sum of all free free heating over model */
484 
485  thermal.thist = 0.;
486  thermal.tlowst = 1e20f;
487 
488  wind.AccelAver = 0.;
489  wind.acldr = 0.;
490  ionbal.ilt = 0;
491  ionbal.iltln = 0;
492  ionbal.ilthn = 0;
493  ionbal.ihthn = 0;
494  ionbal.ifail = 0;
495 
496  secondaries.SecHIonMax = 0.;
497  for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
498  {
499  for( ion=0; ion<nelem+1; ++ion )
500  {
501  supsav[nelem][ion] = secondaries.csupra[nelem][ion];
502  }
503  }
507  for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
508  {
509  for( ion=0; ion<nelem+1; ++ion )
510  {
511  ionbal.CompRecoilHeatRateSave[nelem][ion] = ionbal.CompRecoilHeatRate[nelem][ion];
512  ionbal.CompRecoilIonRateSave[nelem][ion] = ionbal.CompRecoilIonRate[nelem][ion];
513  }
514  }
515 
516  /* these will keep track of the number of convergence failures that occur */
517  conv.nTeFail = 0;
518  conv.nTotalFailures = 0;
519  conv.nPreFail = 0;
520  conv.failmx = 0.;
521  conv.nIonFail = 0;
522  conv.nPopFail = 0;
523  conv.nNeFail = 0;
524  conv.nGrainFail = 0;
525  /* abort flag */
526  lgAbort = false;
527 
528  /* zero out averaging routine */
529  aver("zero",0.,0.," ");
530 
531  GrainStartIter();
532 
533  rfield.comtot = 0.;
534 
535  co.codfrc = 0.;
536  co.codtot = 0.;
537 
538  co.COCoolBigFrac = 0.;
539  co.lgCOCoolCaped = false;
540 
541  hmi.HeatH2DexcMax = 0.;
542  hmi.CoolH2DexcMax = 0.;
543  hmi.h2dfrc = 0.;
544  hmi.h2line_cool_frac = 0.;
545  hmi.h2dtot = 0.;
546  thermal.HeatLineMax = 0.;
547  thermal.GBarMax = 0.;
548  hyperfine.cooling_max = 0.;
549  hydro.cintot = 0.;
550  geometry.lgZoneTrp = false;
551 
552  gv.lgNegGrnDrg = false;
553  gv.TotalDustHeat = 0.;
554  gv.dphmax = 0.;
555  hmi.h2pmax = 0.;
556  gv.dclmax = 0.;
557  gv.GrnElecDonateMax = 0.;
558  gv.GrnElecHoldMax = 0.;
559 
560  /************************************************************************
561  *
562  * allocate space for lines arrays
563  *
564  ************************************************************************/
565 
566  /* there are three types of calls to lines()
567  * ipass = -1, first call, only count number off lines
568  * ipass = 0, second pass, save labels and wavelengths
569  * ipass = 1, integrate intensity*/
570  LineSave.ipass = -1;
571  lines();
572  ASSERT( LineSave.nsum > 0 );
573 
574  /* in a grid this may not be the first time here, return old memory
575  * and grab new appropriate for this size, since number of lines to
576  * be stored can change
577  * as now coded this memory is freed and reallocated on every iteration.
578  * this allows some details of line emission to change for last iteration,
579  * but may not really be necessary */
580  if( LineSv!=NULL )
581  {
582  if( trace.lgTrace )
583  {
584  fprintf( ioQQQ, "IterStart has freed LindSv \n" );
585  }
586  free( LineSv );
587  LineSv = NULL;
588  }
589 
590  /* this is the large main line array */
591  LineSv = (LinSv*)MALLOC((unsigned)LineSave.nsum*sizeof( LinSv ) );
592 
593  for( i=0; i < LineSave.nsum; i++ )
594  {
595  /* this is array of wavelengths, only defined for some lines,
596  * negative is sentinel saying whether line should be added into total continuum */
597  LineSv[i].wavelength = -1;
598  }
599 
600  /* there are three calls to lines()
601  * first call with ipass = -1 was above, only counted number
602  * of lines and malloced space. This is second call and will do
603  * one-time creation of line labels */
604  LineSave.ipass = 0;
605  /* this will count number of lines */
606  LineSave.nsum = 0;
607  lines();
608  /* has to be positive */
609  ASSERT( LineSave.nsum > 0);
610  /* in the future calls to lines will result in integrations */
611  LineSave.ipass = 1;
612 
613  if( trace.lgTrace )
614  fprintf( ioQQQ, "%7ld lines printed in main line array\n",
615  LineSave.nsum );
616 
617  /* now zero out some variables set by last call to LINES */
618  hydro.cintot = 0.;
619  thermal.totcol = 0.;
620  rfield.comtot = 0.;
622  thermal.power = 0.;
623  thermal.HeatLineMax = 0.;
624  thermal.GBarMax = 0.;
625  hyperfine.cooling_max = 0.;
626 
627  gv.TotalDustHeat = 0.;
628  gv.dphmax = 0.;
629  hmi.h2pmax = 0.;
630  gv.dclmax = 0.;
631  gv.GrnElecDonateMax = 0.;
632  gv.GrnElecHoldMax = 0.;
633 
634  co.codfrc = 0.;
635  co.codtot = 0.;
636 
637  hmi.h2dfrc = 0.;
638  hmi.h2line_cool_frac = 0.;
639  hmi.h2dtot = 0.;
640  timesc.sound = 0.;
641 
642  if( LineSave.lgNormSet )
643  {
644  /* set normalization line index to zero from negative initial value so that
645  * we pass the assert in cdLine, in order to get the norm line index */
646  LineSave.ipNormWavL = 0;
647  LineSave.ipNormWavL = cdLine( LineSave.chNormLab , LineSave.WavLNorm , &fhe1nx , &ratio );
648  if( LineSave.ipNormWavL < 0 )
649  {
650  /* did not find the line if return is negative */
651  fprintf( ioQQQ, "PROBLEM could not find the normalisation line.\n");
652  fprintf( ioQQQ,
653  "IterStart could not find a line with a wavelength of %f and a label of %s\n",
655  fprintf( ioQQQ, "Please check the emission line output to find the correct line identification.\n");
656  fprintf( ioQQQ, "Sorry.\n");
657  cdEXIT(EXIT_FAILURE);
658  }
659  }
660  else
661  {
662  /* set normalization line index to zero from negative initial value so that
663  * we pass the assert in cdLine, in order to get the norm line index */
664  LineSave.ipNormWavL = 0;
665  LineSave.ipNormWavL = cdLine( "TOTL" , 4861. , &fhe1nx , &ratio );
666  if( LineSave.ipNormWavL < 0 )
667  {
668  /* did not find the line if return is negative */
669  fprintf( ioQQQ, "PROBLEM could not find the default normalisation line.\n");
670  fprintf( ioQQQ,
671  "IterStart could not find a line with a wavelength of 4861 and a label of TOTL\n" );
672  fprintf( ioQQQ, "Please check the emission line output to find the correct line identification.\n");
673  fprintf( ioQQQ, "Sorry.\n");
674  cdEXIT(EXIT_FAILURE);
675  }
676  }
677 
678  /* set up stop line command on first iteration
679  * find index for lines and save for future iterations
680  * StopCalc.nstpl is zero (false) if no stop line commands entered */
681  if( iteration == 1 && StopCalc.nstpl )
682  {
683  /* nstpl is number of stop line commands, 0 if none entered */
684  for( long int nStopLine=0; nStopLine < StopCalc.nstpl; nStopLine++ )
685  {
686  double relint, absint ;
687  /* returns array index for line in array stack if we found the line,
688  * return negative of total number of lines as debugging aid if line not found */
689  StopCalc.ipStopLin1[nStopLine] = cdLine( StopCalc.chStopLabel1[nStopLine],
690  /* wavelength of line in angstroms, not format printed by code */
691  StopCalc.StopLineWl1[nStopLine], &relint, &absint );
692 
693  if( StopCalc.ipStopLin1[nStopLine]<0 )
694  {
695  fprintf( ioQQQ,
696  " IterStart could not find first line in STOP LINE command, number %ld with label *%s* and wl %.1f\n",
697  StopCalc.ipStopLin1[nStopLine] ,
698  StopCalc.chStopLabel1[nStopLine],
699  StopCalc.StopLineWl1[nStopLine]);
700  cdEXIT(EXIT_FAILURE);
701  }
702 
703  StopCalc.ipStopLin2[nStopLine] = cdLine( StopCalc.chStopLabel2[nStopLine],
704  /* wavelength of line in angstroms, not format printed by code */
705  StopCalc.StopLineWl2[nStopLine], &relint, &absint );
706 
707  if( StopCalc.ipStopLin2[nStopLine] < 0 )
708  {
709  fprintf( ioQQQ,
710  " IterStart could not find second line in STOP LINE command, line number %ld with label *%s* and wl %.1f\n",
711  StopCalc.ipStopLin1[nStopLine] ,
712  StopCalc.chStopLabel1[nStopLine],
713  StopCalc.StopLineWl1[nStopLine]);
714  cdEXIT(EXIT_FAILURE);
715  }
716 
717  if( trace.lgTrace )
718  {
719  fprintf( ioQQQ,
720  " stop line 1 is number %5ld wavelength is %f label is %4.4s\n",
721  StopCalc.ipStopLin1[nStopLine],
722  LineSv[StopCalc.ipStopLin1[nStopLine]].wavelength,
723  LineSv[StopCalc.ipStopLin1[nStopLine]].chALab );
724  fprintf( ioQQQ,
725  " stop line 2 is number %5ld wavelength is %f label is %4.4s\n",
726  StopCalc.ipStopLin2[nStopLine],
727  LineSv[StopCalc.ipStopLin2[nStopLine]].wavelength,
728  LineSv[StopCalc.ipStopLin2[nStopLine]].chALab );
729  }
730  }
731  }
732 
733  /* option to only print last iteration */
734  if( prt.lgPrtLastIt )
735  {
736  if( iteration == 1 )
737  {
738  called.lgTalk = false;
739  }
740 
741  /* initial condition of TALK may be off if AMOEBA used
742  * sec part is for print last command
743  * lgTalkForcedOff is set true when cdTalk is called with false
744  * to turn off printing */
746  {
748  }
749  }
750 
751  if( trace.lgTrace && trace.lgOptcBug )
752  {
753  if( iteration > 1 )
754  {
755  fprintf( ioQQQ, " Outer optical depths defined.\n" );
756  }
757  else
758  {
759  fprintf( ioQQQ, " Outer optical depths NOT defined.\n" );
760  }
761 
762  /* heavy element lines */
763  for( i=0; i < nLevel1; i++ )
764  {
765  fprintf( ioQQQ, "%6f:%8.1e%8.1e%8.1e",
766  TauLines[i].WLAng,
767  TauLines[i].Emis->TauIn,
768  TauLines[i].Emis->TauTot,
769  TauLines[i].Emis->Pesc );
770  }
771 
772  /*Atomic and molecular lines-Humeshkar Nemala*/
773  for( i=0; i <linesAdded2; i++)
774  {
775  fprintf( ioQQQ, "%6f:%8.1e%8.1e%8.1e",
776  atmolEmis[i].tran->WLAng,
777  atmolEmis[i].TauIn,
778  atmolEmis[i].TauTot,
779  atmolEmis[i].Pesc);
780  }
781 
782  fprintf( ioQQQ, "\n" );
783  }
784 
785  if( opac.lgCaseB )
786  {
787  if( trace.lgTrace )
788  {
789  fprintf( ioQQQ, " IterStart does not change mid-zone optical depth since CASE B\n" );
790  }
791  }
792 
793  /* check if induced recombination can be ignored */
794  hydro.FracInd = 0.;
795  hydro.fbul = 0.;
796 
797  /* remember some things about effects of induced rec on H only
798  * don't do ground state since SPHERE turns it off */
799  for( i=ipH2p; i < (iso.numLevels_max[ipH_LIKE][ipHYDROGEN] - 1); i++ )
800  {
801  if( Transitions[ipH_LIKE][ipHYDROGEN][i+1][i].Emis->Aul <= iso.SmallA )
802  continue;
803 
808  if( ratio > hydro.FracInd )
809  {
810  hydro.FracInd = (realnum)ratio;
811  hydro.ndclev = i;
812  }
813 
814  ratio = Transitions[ipH_LIKE][ipHYDROGEN][i+1][i].Emis->pump/
815  (Transitions[ipH_LIKE][ipHYDROGEN][i+1][i].Emis->pump +
816  Transitions[ipH_LIKE][ipHYDROGEN][i+1][i].Emis->Aul);
817 
818  if( ratio > hydro.fbul )
819  {
820  hydro.fbul = (realnum)ratio;
821  hydro.nbul = i;
822  }
823  }
824 
825  /* this was set in call to lines above */
826  ASSERT( LineSave.nsum > 0);
827 
828  if( trace.lgTrace )
829  fprintf( ioQQQ, " IterStart returns.\n" );
830  return;
831 }
832 
833 /*IterRestart reset many variables at the start of a new iteration
834  * called by cloudy after calculation is completed, when more iterations
835  * are needed - the iteration has been incremented when this routine is
836  * called so iteration == 2 after first iteration, we are starting
837  * the second iteration */
838 void IterRestart(void)
839 {
840  long int i,
841  ion,
842  ipISO ,
843  n,
844  nelem;
845  double SumOTS;
846 
847  DEBUG_ENTRY( "IterRestart()" );
848 
850 
851  /* this is case where temperature floor has been set, if it was hit
852  * then we did a constant temperature calculation, and must go back to
853  * a thermal solution
854  * test on thermal.lgTemperatureConstantCommandParsed distinguishes
855  * from temperature floor option, so not reset if constant temperature
856  * was actually set */
858  {
860  thermal.ConstTemp = 0.;
861  }
862 
863  lgAbort = false;
864 
865  /* reset some parameters needed for magnetic field */
866  Magnetic_reinit();
867 
868  opac.stimax[0] = 0.;
869  opac.stimax[1] = 0.;
870 
871  H2_Reset();
872 
873  ca.oldcla = 0.;
874 
875  for( nelem=ipHYDROGEN; nelem < LIMELM; ++nelem )
876  {
877  /* these have one more ion than above */
878  for( ion=0; ion<nelem+2; ++ion )
879  {
880  long int ion2;
881  /* zero out the source and sink arrays */
882  mole.source[nelem][ion] = saveMoleSource[nelem][ion];
883  mole.sink[nelem][ion] = saveMoleSink[nelem][ion];
884  /*>>chng 06 jun 27, move from here, where leak occurs, to
885  * above where xMoleChTrRate first created */
886  /*mole.xMoleChTrRate[nelem][ion] = (realnum*)MALLOC(sizeof(realnum)*(unsigned)(nelem+2));*/
887  for( ion2=0; ion2<nelem+2; ++ion2 )
888  {
889  mole.xMoleChTrRate[nelem][ion][ion2] = SaveMoleChTrRate[nelem][ion][ion2];
890  }
891  }
892  }
893 
894  /* reset molecular abundances */
895  for( i=0; i<N_H_MOLEC; i++)
896  {
897  hmi.Hmolec[i] = hmsav[i];
898  }
900  /*fprintf(ioQQQ," IterRestar sets H2 total to %.2e\n",hmi.H2_total );*/
903 
904  /* zero out the column densities */
905  for( i=0; i < NCOLD; i++ )
906  {
907  colden.colden[i] = 0.;
908  }
909  colden.He123S = 0.;
910  colden.coldenH2_ov_vel = 0.;
911 
912  for( i=0; i < mole.num_comole_calc; i++ )
913  {
914  /* column densities */
915  COmole[i]->hevcol = 0.;
916  /* largest fraction of atoms in molecules */
917  COmole[i]->xMoleFracMax = 0.;
918  }
919 
920  for( i=0; i< mole.num_comole_calc; ++i )
921  {
922  COmole[i]->hevmol = COmole[i]->hevmol_save;
923 
924  /* check for NaN */
925  ASSERT( !isnan( COmole[i]->hevmol ) );
926  }
927 
928  /* close out this iteration if dynamics or time dependence is enabled */
929  if( dynamics.lgAdvection )
930  DynaEndIter();
931 
936 
939 
940  hmi.BiggestH2 = -1.f;
944 
951 
955 
957  rfield.EnergyBremsThin = 0.;
958  rfield.lgUSphON = false;
959 
960  radius.glbdst = 0.;
961  /* zone thickness, and next zone thickness */
962  radius.drad = drSave;
964 
965  /* was min dr ever used? */
966  radius.lgDrMinUsed = false;
967 
968  /* simple estimate of H-beta intensity */
971 
972  /* total luminosity */
974 
975  /* set debugger on now if NZONE desired is 0 */
976  if( (trace.nznbug == 0 && iteration >= trace.npsbug) && trace.lgTrOvrd )
977  {
978  if( trace.nTrConvg==0 )
979  {
980  trace.lgTrace = true;
981  }
982  else
983  /* trace convergence entered = but with negative flag = make positive,
984  * abs and not mult by -1 since may trigger more than one time */
985  trace.nTrConvg = abs( trace.nTrConvg );
986 
987  fprintf( ioQQQ, " IterRestart called.\n" );
988  }
989  else
990  {
991  trace.lgTrace = false;
992  }
993 
994  /* zero secondary suprathermals variables for first ionization attem */
995  for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
996  {
997  for( ion=0; ion<nelem+1; ++ion )
998  {
999  secondaries.csupra[nelem][ion] = supsav[nelem][ion];
1000  }
1001  }
1005  for( nelem=0; nelem<LIMELM; ++nelem)
1006  {
1007  for( ion=0; ion<nelem+1; ++ion )
1008  {
1009  ionbal.CompRecoilHeatRate[nelem][ion] = ionbal.CompRecoilHeatRateSave[nelem][ion];
1010  ionbal.CompRecoilIonRate[nelem][ion] = ionbal.CompRecoilIonRateSave[nelem][ion];
1011  }
1012  }
1013 
1014  wind.lgVelPos = true;
1015  wind.AccelMax = 0.;
1016  wind.AccelAver = 0.;
1017  wind.acldr = 0.;
1018  wind.windv = wind.windv0;
1019 
1020  thermal.nUnstable = 0;
1021  thermal.lgUnstable = false;
1022 
1023  pressure.pbeta = 0.;
1024  pressure.RadBetaMax = 0.;
1025  pressure.lgPradCap = false;
1026  pressure.lgPradDen = false;
1027  /* this flag will say we hit the sonic point */
1028  pressure.lgSonicPoint = false;
1029  pressure.lgStrongDLimbo = false;
1030  pressure.PresInteg = 0.;
1031  pressure.pinzon = 0.;
1032 
1033  dense.eden = edsav;
1034  dense.EdenHCorr = edsav;
1035  dense.EdenTrue = edsav;
1037 
1038  for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
1039  {
1040  /* reset molecular densities */
1041  dense.gas_phase[nelem] = gas_phase_save[nelem];
1042  dense.xMolecules[nelem] = xMolFracsSave[nelem];
1043  dense.IonLow[nelem] = IonLowSave[nelem];
1044  dense.IonHigh[nelem] = IonHighSave[nelem];
1045  }
1046  /*fprintf(ioQQQ,"DEBUG IterRestart gas_phase set to save hden %.4e\n",
1047  dense.gas_phase[0]);*/
1048 
1049  for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
1050  {
1051  for( nelem=ipISO; nelem < LIMELM; nelem++ )
1052  {
1053  if( dense.lgElmtOn[nelem] )
1054  {
1055  for( n=ipH1s; n < iso.numLevels_max[ipISO][nelem]; n++ )
1056  {
1057  iso.ConOpacRatio[ipISO][nelem][n] = HOpacRatSav[ipISO][nelem][n];
1058  StatesElem[ipISO][nelem][n].Pop = hnsav[ipISO][nelem][n];
1059  }
1060  }
1061  }
1062  }
1063 
1064  if( trace.lgTrace && trace.lgNeBug )
1065  {
1066  fprintf( ioQQQ, " EDEN set to%12.4e by IterRestart.\n",
1067  dense.eden );
1068  }
1069 
1070 # if 0
1071  for( nelem=ipHELIUM; nelem < LIMELM; nelem++ )
1072  {
1073  if( dense.lgElmtOn[nelem] )
1074  {
1075  /* reset total gas phase abundance to the original set */
1076  dense.gas_phase[nelem] = abund.solar[nelem]*dense.gas_phase[ipHYDROGEN];
1077  }
1078  else
1079  {
1080  /* >>chng 04 apr 20, set to zero if element is off */
1081  dense.gas_phase[nelem] = 0.;
1082  }
1083  }
1084 # endif
1085  for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
1086  {
1087  for( ion=0; ion < (nelem + 2); ion++ )
1088  {
1089  dense.xIonDense[nelem][ion] = xIonFsave[nelem][ion];
1090  }
1091  for( i=0; i < LIMELM; i++ )
1092  {
1093  thermal.heating[nelem][i] = HeatSave[nelem][i];
1094  }
1095  }
1096 
1097  GrainRestartIter();
1098 
1099  /* continuum was saved in flux_total_incident */
1100  for( i=0; i < rfield.nupper; i++ )
1101  {
1102  /* time-constant part of beamed continuum */
1104  /* continuum flux_time_beam_save has the initial value of the
1105  * time-dependent beamed continuum */
1108  /* the isotropic continuum source is time steady */
1112 
1114  rfield.SummedCon[i] = rfield.flux[i] + rfield.SummedDif[i];
1115  rfield.trans_coef_zone[i] = 1.f;
1116  rfield.trans_coef_total[i] = 1.f;
1117 
1119  rfield.otscon[i] = rfield.otssav[i][0];
1120  rfield.otslin[i] = rfield.otssav[i][1];
1121  /* >>chng 01 apr 10, add these two resets */
1122  /* >>chng 03 may 20, had set to otscon and otslin, which is whatever
1123  * was left from previous iteration - change to zero */
1124  rfield.otsconNew[i] = 0.;
1125  rfield.otslinNew[i] = 0.;
1126  rfield.ConInterOut[i] = 0.;
1127  rfield.OccNumbBremsCont[i] = 0.;
1128  rfield.OccNumbDiffCont[i] = 0.;
1129  rfield.OccNumbContEmitOut[i] = 0.;
1130  rfield.outlin[i] = 0.;
1131  rfield.outlin_noplot[i] = 0.;
1132  rfield.ConOTS_local_OTS_rate[i] = 0.;
1133  rfield.ConOTS_local_photons[i] = 0.;
1134 
1138  opac.albedo[i] =
1139  opac.opacity_sct[i]/MAX2(1e-30,opac.opacity_sct[i] + opac.opacity_abs[i]);
1140  opac.tmn[i] = 1.;
1141  /* >>chng 99 dec 04, having exactly 1 for first zone caused discontinuity
1142  * for heating in very high T models in func_map.in. zone 1 and 2 were 20% different,
1143  * since tau in is 1e-20, e2 is 0.9999, and so some H ots
1144  opac.ExpmTau[i] = 1.;
1145  opac.E2TauAbsFace[i] = 1.;*/
1146  /* >>chng 99 dec 04, having exactly 1 for first zone caused discontinuity
1147  * for heating in very high T models in func_map.in. zone 1 and 2 were 20% different,
1148  * since tau in is 1e-20, e2 is 0.9999, and so some H ots
1149  * these were not here at all*/
1150  /* attenuation of flux by optical depths IN THIS ZONE
1151  * DirectionalCosin is 1/COS(theta), is usually 1, reset with illuminate command,
1152  * option for illumination of slab at an angle */
1154 
1155  /* e(-tau) in inward direction, up to illuminated face */
1156  opac.ExpmTau[i] = (realnum)opac.ExpZone[i];
1157 
1158  /* e2(tau) in inward direction, up to illuminated face */
1160  rfield.reflin[i] = 0.;
1161  rfield.ConEmitReflec[i] = 0.;
1162  rfield.ConEmitOut[i] = 0.;
1163  rfield.ConRefIncid[i] = 0.;
1164 
1165  /* escape in the outward direction
1166  * on second and later iterations define outward E2 */
1167  if( iteration > 1 )
1168  {
1169  /* e2 from current position to outer edge of shell */
1171  opac.E2TauAbsOut[i] = (realnum)e2( tau );
1172  ASSERT( opac.E2TauAbsOut[i]>=0. && opac.E2TauAbsOut[i]<=1. );
1173  }
1174  else
1175  opac.E2TauAbsOut[i] = 1.;
1176 
1177  }
1178 
1179  /* update continuum */
1180  RT_OTS_Update(&SumOTS , 0.);
1181 
1182  thermal.FreeFreeTotHeat = 0.;
1183 
1184  if( called.lgTalk )
1185  {
1186  fprintf( ioQQQ, "\f\n Start Iteration Number%2ld %75.75s\n\n\n",
1187  iteration, input.chTitle );
1188  }
1189 
1190  /* reset variable to do with FeII atom, in FeIILevelPops */
1191  FeIIReset();
1192  return;
1193 }
1194 
1195 /* do some work with ending the iteration */
1196 void IterEnd(void)
1197 {
1198  long i;
1199  double fac,
1200  tau;
1201 
1202  DEBUG_ENTRY( "IterEnd()" );
1203 
1204  /* give indication of geometry */
1205  fac = radius.depth/radius.rinner;
1206  if( fac < 0.1 )
1207  {
1208  geometry.lgGeoPP = true;
1209  }
1210  else
1211  {
1212  geometry.lgGeoPP = false;
1213  }
1214 
1215  /* >>chng 06 apr 03, add last radius and drad */
1216  for( i=0; i<struc.nzone; ++i )
1217  {
1218  struc.depth_last[i] = struc.depth[i];
1219  struc.drad_last[i] = struc.drad[i];
1220  }
1222 
1223  /* all continue were attenuated in last zone in radius_increment to represent the intensity
1224  * in the middle of the next zone - this was too much since we now want
1225  * intensity emergent from outer edge of last zone */
1226  for( i=0; i < rfield.nflux; i++ )
1227  {
1228  {
1229  enum{DEBUG_LOC=false};
1230  if( DEBUG_LOC)
1231  {
1232  fprintf(ioQQQ,"i=%li opac %.2e \n", i,
1234  }
1235  }
1237  ASSERT( tau > 0. );
1238  fac = sexp( tau );
1239 
1240  /* >>chng 02 dec 14, add first test to see whether product of ratio is within
1241  * range of a float - ConInterOut can be finite and fac just above a small float,
1242  * so ratio exceeds largest size of a float */
1243  /*if( fac > SMALLFLOAT )*/
1244  if( (realnum)(fac/SDIV(rfield.ConInterOut[i]))>SMALLFLOAT && fac > SMALLFLOAT )
1245  {
1246  rfield.ConInterOut[i] /= (realnum)fac;
1247  rfield.outlin[i] /= (realnum)fac;
1248  rfield.outlin_noplot[i] /= (realnum)fac;
1249  }
1250  }
1251 
1252  /* remember thickness of previous iteration */
1254  return;
1255 }

Generated for cloudy by doxygen 1.8.3.1