cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
parse_save.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*ParseSave parse the save command */
4 /*SaveFilesInit initialize save file pointers, called from cdInit */
5 /*CloseSaveFiles close all save files */
6 /*ChkUnits check for keyword UNITS on line, then scan wavelength or energy units if present */
7 #include "cddefines.h"
8 #include "cddrive.h"
9 #include "physconst.h"
10 #include "elementnames.h"
11 #include "input.h"
12 #include "geometry.h"
13 #include "prt.h"
14 #include "optimize.h"
15 #include "rfield.h"
16 #include "hcmap.h"
17 #include "atomfeii.h"
18 #include "h2.h"
19 #include "mole.h"
20 #include "hmi.h"
21 #include "version.h"
22 #include "grainvar.h"
23 #include "parse.h"
24 #include "grid.h"
25 #include "save.h"
26 #include "mpi_utilities.h"
27 #include "parser.h"
28 
29 /* check for keyword UNITS on line, then scan wavelength or energy units if present */
30 STATIC void ChkUnits(Parser &p);
31 
32 /* NB NB NB NB NB NB NB NB NB NB
33  *
34  * if any "special" save commands are added (commands that copy the file pointer
35  * into a separate variable, e.g. like SAVE _DR_), be sure to add that file pointer
36  * to SaveFilesInit and CloseSaveFiles !!!
37  *
38  * SAVE FILE POINTERS SHOULD NEVER BE ALTERED BY ROUTINES OUTSIDE THIS MODULE !!!
39  *
40  * hence initializations of save file pointers should never be included in zero() !!
41  * the pointer might be lost without the file being closed
42  *
43  * NB NB NB NB NB NB NB NB NB NB */
44 
45 /* save file header headers - these are written into the string save.chHeader[save.nsave] when
46  * the command is parsed
47  * save.lgPunHeader[] determines whether header is saved
48  */
49 
50 
51 void ParseSave(Parser& p)
52 {
53  char chLabel[INPUT_LINE_LENGTH] ,
54  chFilename[INPUT_LINE_LENGTH] ,
55  chSecondFilename[INPUT_LINE_LENGTH];
56  bool lgSecondFilename;
57  /* pointer to column across line image for free format number reader*/
58  long int i,
59  nelem;
60 
61  char chTemp[MAX_HEADER_SIZE];
62 
63  DEBUG_ENTRY( "ParseSave()" );
64 
65  /* check that limit not exceeded */
66  if( save.nsave >= LIMPUN )
67  {
68  fprintf( ioQQQ,
69  "The limit to the number of SAVE options is %ld. Increase "
70  "LIMPUN in save.h if more are needed.\nSorry.\n",
71  LIMPUN );
73  }
74 
75  /* initialize this flag, forced true for special cases below (e.g. for FITS files) */
77 
78  /* LAST keyword is an option to save only on last iteration */
79  save.lgPunLstIter[save.nsave] = p.nMatch("LAST");
80 
81  /* get file name for this save output.
82  * GetQuote does the following -
83  * first copy original version of file name into chLabel,
84  * string does include null termination.
85  * set filename in OrgCard and second parameter to spaces so
86  * that not picked up below as keyword
87  * last parameter says whether to abort if no quote found */
88  if( p.GetQuote( chLabel , true ) )
89  /* this can't happen since routine would not return at all if no double quotes found */
90  TotalInsanity();
91 
92  /* check that name is not same as opacity.opc, a special file */
93  if( strcmp(chLabel , "opacity.opc") == 0 )
94  {
95  fprintf( ioQQQ, "ParseSave will not allow save file name %s, please choose another.\nSorry.\n",
96  chLabel);
98  }
99  else if( chLabel[0]=='\0' )
100  {
101  fprintf( ioQQQ, "ParseSave found a null file name between double quotes, please check command line.\nSorry.\n");
103  }
104 
105  /* now copy to chFilename, with optional grid prefix first */
106  strcpy( chFilename , save.chGridPrefix.c_str() );
107  /* this is optional prefix, normally a null string, set with set save prefix command */
108  strcat( chFilename , save.chFilenamePrefix.c_str() );
109  strcat( chFilename , chLabel );
110 
111  /* there may be a second file name, and we need to get it off the line
112  * before we parse options, last false parameter says not to abort if
113  * missing - this is not a problem at this stage */
114  if( p.GetQuote( chSecondFilename , false ) )
115  lgSecondFilename = false;
116  else
117  lgSecondFilename = true;
118 
119  /* CLOBBER clobber keyword is an option to overwrite rather than
120  * append to a given file */
121  if( p.nMatch("CLOB") )
122  {
123  if( p.nMatch(" NO ") )
124  {
125  /* do not clobber files */
126  save.lgNoClobber[save.nsave] = true;
127  }
128  else
129  {
130  /* clobber files */
131  save.lgNoClobber[save.nsave] = false;
132  }
133  }
134 
135  /* put version number and title of model on output file, but only if
136  * this is requested with a "title" on the line*/
137  /* >>chng 02 may 10, invert logic from before - default had been title */
138  /* put version number and title of model on output file, but only if
139  * there is not a "no title" on the line*/
140  if( !p.nMatch("NO TI") && p.nMatch("TITL"))
141  {
142  sprintf( save.chHeader[save.nsave],
143  "# %s %s\n",
145  }
146 
147  /* usually results for each iteration are followed by a series of
148  * hash marks, ####, which fool excel. This is an option to not print
149  * the line. If the keyword NO HASH no hash appears the hash marks
150  * will not occur */
151  if( p.nMatch("NO HA") )
152  save.lgHashEndIter[save.nsave] = false;
153 
154  /* save opacity must come first since elements appear as sub-keywords */
155  if( p.nMatch("OPAC") )
156  {
157  /* check for keyword UNITS on line, then scan wavelength or energy units if present,
158  * units are copied into save.chConPunEnr */
159  ChkUnits(p);
160 
161  strcpy( save.chSave[save.nsave], "OPAC" );
162 
163  /* "every" option to save this on every zone -
164  * not present then only last zone is saved */
165  if( p.nMatch("EVER" ) )
166  {
167  /* save every zone */
168  save.lgSaveEveryZone[save.nsave] = true;
170  }
171  else
172  {
173  /* only save last zone */
174  save.lgSaveEveryZone[save.nsave] = false;
176  }
177 
178  if( p.nMatch("TOTA") )
179  {
180  /* DoPunch will call save_opacity to parse the subcommands
181  * save total opacity */
182  strcpy( save.chOpcTyp[save.nsave], "TOTL" );
183  sprintf( save.chHeader[save.nsave],
184  "#nu/%s\tTot opac\tAbs opac\tScat opac\tAlbedo\telem\n",
186  }
187 
188  else if( p.nMatch("FIGU") )
189  {
190  /* do figure for hazy */
191  strcpy( save.chOpcTyp[save.nsave], "FIGU" );
192  sprintf( save.chHeader[save.nsave],
193  "#nu/%s\tH\tHe\ttot opac\n",
195  }
196 
197  else if( p.nMatch("FINE") )
198  {
199  /* save the fine opacity array */
200  rfield.lgSaveOpacityFine = true;
201  strcpy( save.chOpcTyp[save.nsave], "FINE" );
202  /* check for keyword UNITS on line, then scan wavelength or energy units if present,
203  * units are copied into save.chConPunEnr */
204  ChkUnits(p);
205 
206  sprintf( save.chHeader[save.nsave],
207  "#nu/%s\topac\n",
209 
210  /* range option - important since so much data - usually want to
211  * only give portion of the continuum */
212  if( p.nMatch("RANGE") )
213  {
214  /* get lower and upper range, eventually must be in Ryd */
215  double Energy1 = p.FFmtRead();
216  double Energy2 = p.FFmtRead();
217  if( p.lgEOL() )
218  {
219  fprintf(ioQQQ,"There must be two numbers, the lower and upper energy range in Ryd.\nSorry.\n");
221  }
222  if( p.nMatch("UNIT" ) )
223  {
224  // apply units to range option
225  const char *energyUnits = p.StandardEnergyUnit();
226  Energy unitChange;
227  unitChange.set(Energy1, energyUnits );
228  Energy1 = unitChange.Ryd();
229  unitChange.set(Energy2, energyUnits );
230  Energy2 = unitChange.Ryd();
231  }
232  /* get lower and upper rang in Ryd */
233  save.punarg[save.nsave][0] = (realnum)MIN2( Energy1 , Energy2 );
234  save.punarg[save.nsave][1] = (realnum)MAX2( Energy1 , Energy2 );
235  //fprintf(ioQQQ , "DEBUG units change fine %.3e %.3e\n" , save.punarg[save.nsave][0] ,
236  // save.punarg[save.nsave][1] );
237  //cdEXIT(EXIT_FAILURE);
238  }
239  else
240  {
241  /* these mean full energy range */
242  save.punarg[save.nsave][0] = 0.;
243  save.punarg[save.nsave][1] = 0.;
244  }
245  /* optional last parameter - how many points to bring together */
246  save.punarg[save.nsave][2] = (realnum)p.FFmtRead();
247 
248  /* default is to average together ten */
249  if( p.lgEOL() )
250  save.punarg[save.nsave][2] = 10;
251 
252  if( save.punarg[save.nsave][2] < 1 )
253  {
254  fprintf(ioQQQ,"The number of fine opacities to skip must be > 0 \nSorry.\n");
256  }
257  }
258 
259  else if( p.nMatch("GRAI") )
260  {
261  /* save grain opacity command, give optical properties of gains in calculation */
262  strcpy( save.chSave[save.nsave], "DUSO" );
263  /* save grain opacity command in twice, here and above in opacity */
264  sprintf( save.chHeader[save.nsave],
265  "#grain\tnu\tabs+scat*(1-g)\tabs\tscat*(1-g)\tscat\tscat*(1-g)/[abs+scat*(1-g)]\n" );
266  }
267 
268  else if( p.nMatch("BREM") )
269  {
270  /* save bremsstrahlung opacity */
271  strcpy( save.chOpcTyp[save.nsave], "BREM" );
272  sprintf( save.chHeader[save.nsave],
273  "#nu\tbrem opac\n" );
274  }
275 
276  else if( p.nMatch("SHEL") )
277  {
278  /* save shells, a form of the save opacity command for showing subshell crossections*/
279  strcpy( save.chSave[save.nsave], "OPAC" );
280 
281  /* save subshell cross sections */
282  strcpy( save.chOpcTyp[save.nsave], "SHEL" );
283 
284  /* this is element */
285  save.punarg[save.nsave][0] = (realnum)p.FFmtRead();
286 
287  /* this is ion */
288  save.punarg[save.nsave][1] = (realnum)p.FFmtRead();
289 
290  /* this is shell */
291  save.punarg[save.nsave][2] = (realnum)p.FFmtRead();
292 
293  if( p.lgEOL() )
294  {
295  fprintf( ioQQQ, "There must be atom number, ion, shell\nSorry.\n" );
297  }
298  sprintf( save.chHeader[save.nsave],
299  "#sub shell cross section\n" );
300  }
301 
302  else if( p.nMatch("ELEM") )
303  {
304  /* save element opacity, produces n name.n files, one for each stage of
305  * ionization. the name is the 4-char version of the element's name, and
306  * n is the stage of ionization. the file name on the card is ignored.
307  * The code stops in save_opacity after these files are produced. */
308 
309  /* this will be used as check that we did find an element on the command lines */
310  /* nelem is -1 if an element was not found */
311  if( (nelem = p.GetElem() ) < 0 )
312  {
313  fprintf( ioQQQ, "I did not find an element name on the opacity element command. Sorry.\n" );
315  }
316 
317  /* copy string over */
319  }
320  else
321  {
322  fprintf( ioQQQ, " I did not recognize a keyword on this save opacity command.\n" );
323  fprintf( ioQQQ, " Sorry.\n" );
325  }
326  }
327 
328  /* save H2 has to come early since it has many suboptions */
329  else if( p.nMatchErase(" H2 ") )
330  {
331  /* this is in mole_h2_io.c */
333  }
334 
335  /* save HD has to come early since it has many suboptions */
336  else if( p.nMatchErase(" HD ") )
337  {
338  /* this is in mole_h2_io.c */
340  }
341 
342  /* save grain abundance will be handled later */
343  else if( p.nMatch("ABUN") && !p.nMatch("GRAI") )
344  {
345  /* save abundances */
346  strcpy( save.chSave[save.nsave], "ABUN" );
347  sprintf( save.chHeader[save.nsave],
348  "#abund H" );
349  for(nelem=ipHELIUM;nelem<LIMELM; ++nelem )
350  {
351  sprintf( chTemp, "\t%s",
353  strcat( save.chHeader[save.nsave], chTemp );
354  }
355  strcat( save.chHeader[save.nsave], "\n");
356  }
357 
358  else if( p.nMatch(" AGE") )
359  {
360  /* save ages */
361  strcpy( save.chSave[save.nsave], "AGES" );
362  sprintf( save.chHeader[save.nsave],
363  "#ages depth\tt(cool)\tt(H2 dest)\tt(CO dest)\tt(OH dest)\tt(H rec)\n" );
364  }
365 
366  else if( p.nMatch(" AGN") )
367  {
368  /* save tables needed for AGN3 */
369  strcpy( save.chSave[save.nsave], " AGN" );
370  /* this is the AGN option, to produce a table for AGN */
371 
372  /* charge exchange rate coefficients */
373  if( p.nMatch("CHAR") )
374  {
375  strcpy( save.chSave[save.nsave], "CHAG" );
376  sprintf( save.chHeader[save.nsave],
377  "#charge exchange rate coefficnt\n" );
378  }
379 
380  else if( p.nMatch("RECO") )
381  {
382  /* save recombination rates for AGN3 table */
383  strcpy( save.chSave[save.nsave], "RECA" );
384  sprintf( save.chHeader[save.nsave],
385  "#Recom rates for AGN3 table\n" );
386  }
387 
388  else if( p.nMatch("OPAC") )
389  {
390  /* create table for appendix in AGN */
391  strcpy( save.chOpcTyp[save.nsave], " AGN" );
392  strcpy( save.chSave[save.nsave], "OPAC" );
393  }
394 
395  else if( p.nMatch("HECS") )
396  {
397  /* create table for appendix in AGN */
398  strcpy( save.chSaveArgs[save.nsave], "HECS" );
399  sprintf( save.chHeader[save.nsave],
400  "#AGN3 he cs \n" );
401  }
402 
403  else if( p.nMatch("HEMI") )
404  {
405  /* HEMIS - continuum emission needed for chap 4 of AGN3 */
406  strcpy( save.chSaveArgs[save.nsave], "HEMI" );
407 
408  /* check for keyword UNITS on line, then scan wavelength or energy units if present,
409  * units are copied into save.chConPunEnr */
410  ChkUnits(p);
411  }
412  else if( p.nMatch("RECC") )
413  {
414  /* recombination cooling, for AGN */
415  strcpy( save.chSave[save.nsave], "HYDr" );
416  sprintf( save.chHeader[save.nsave],
417  "#T\tbAS\tb1\tbB\n" );
418  }
419  else
420  {
421  fprintf( ioQQQ, " I did not recognize this option on the SAVE HYDROGEN command.\n" );
422  fprintf( ioQQQ, " Sorry.\n" );
424  }
425  }
426 
427  else if( p.nMatch("AVER") )
428  {
429  /* save averages */
430  strcpy( save.chSave[save.nsave], "AVER" );
431  /* no need to print this standard line of explanation*/
432  /*sprintf( save.chHeader[save.nsave], " asserts\n" );*/
433 
434  /* actually get the averages from the input stream, and malloc the
435  * space in the arrays
436  * save io unit not used in read */
438  }
439 
440  /* save charge transfer */
441  else if( p.nMatch("CHAR") && p.nMatch("TRAN") )
442  {
443  /* NB in SaveDo only the first three characters are compared to find this option,
444  * search for "CHA" */
445  /* save charge transfer */
446  strcpy( save.chSave[save.nsave], "CHAR" );
447  sprintf( save.chHeader[save.nsave],
448  "#charge exchange rate coefficient\n" );
449  }
450 
451  // save chianti collision strengths in physical units
452  else if( p.nMatch("CHIA"))
453  {
454  strcpy( save.chSave[save.nsave], "CHIA" );
455 
456  }
457 
458  else if( p.nMatch("CHEM") )
459  {
460 
461  if( p.nMatch( "RATE" ) )
462  {
463  /* >>chng 06 May 30, NPA. Save reaction rates for selected species */
464  if( lgSecondFilename )
465  {
466  if( p.nMatch( "DEST" ) )
467  strcpy( save.chSaveArgs[save.nsave], "DEST" );
468  else if( p.nMatch( "CREA" ) )
469  strcpy( save.chSaveArgs[save.nsave], "CREA" );
470  else if( p.nMatch( "CATA" ) )
471  strcpy( save.chSaveArgs[save.nsave], "CATA" );
472  else if( p.nMatch( "ALL" ) )
473  strcpy( save.chSaveArgs[save.nsave], "ALL " );
474  else
475  strcpy( save.chSaveArgs[save.nsave], "DFLT" );
476 
477  strcpy( save.chSave[save.nsave], "CHRT" );
478  save.optname[save.nsave] = chSecondFilename;
479  // Haven't read chemistry database yet, so put off setting up header
480  //sprintf( save.chHeader[save.nsave], "#");
481  }
482 
483  else
484  {
485  fprintf(ioQQQ," A species label must appear within a second set of quotes (following the output filename).\n" );
486  fprintf( ioQQQ, " Sorry.\n" );
488  }
489 
490  }
491  }
492 
493  else if( p.nMatch("COMP") )
494  {
495  /* save Compton, details of the energy exchange problem */
496  strcpy( save.chSave[save.nsave], "COMP" );
497  sprintf( save.chHeader[save.nsave],
498  "#nu, comup, comdn\n" );
499  }
500 
501  else if( p.nMatch("COOL") )
502  {
503  /* save cooling, actually done by routine cool_save */
504  if( p.nMatch("EACH") )
505  {
506  strcpy( save.chSave[save.nsave], "EACH");
507  sprintf( save.chHeader[save.nsave],
508  "#depth(cm)\tTemp(K)\tCtot(erg/cm3/s)\t" );
509  for( int i = 0 ; i < LIMELM ; i++ )
510  {
512  strcat(save.chHeader[save.nsave], "\t" );
513  }
514  strcat(save.chHeader[save.nsave], "molecule\tdust\tH2cX\tCT C\tH-fb\tH2ln\tHDro\tH2+ \tFFcm\thvFB\teeff\tComp\tExtr\tExpn\tCycl\tHvin\tdima\n" );
515  }
516  else
517  {
518  strcpy( save.chSave[save.nsave], "COOL");
519  /*>>chng 06 jun 06, revise to be same as save cooling */
520  sprintf( save.chHeader[save.nsave],
521  "#depth cm\tTemp K\tHtot erg/cm3/s\tCtot erg/cm3/s\tcool fracs\n" );
522  }
523  }
524 
525  // punch the dominant rates for a given species
526  else if( p.nMatch("DOMI") && p.nMatch("RATE"))
527  {
528  if( !lgSecondFilename )
529  {
530  fprintf( ioQQQ,"This command requires two items in quotes (a filename and a species label). Only one set of quotes was found.\nSorry.\n");
532  }
533  /* in this case the "second filename" is really the species label. */
534  strncpy( save.chSpeciesDominantRates[save.nsave], chSecondFilename, CHARS_SPECIES );
535 
536  /* save dominant rates "species" */
537  strcpy( save.chSave[save.nsave], "DOMI" );
538  sprintf( save.chHeader[save.nsave],
539  "#depth cm\t%s col cm-2\tsrc s-1\tsnk s-1\n",
541  }
542 
543  else if( p.nMatch("DYNA") )
544  {
545  /* save something dealing with dynamics
546  * in SaveDo the DYN part of key is used to call DynaPunch,
547  * with the 4th char as its second argument. DynaSave uses that
548  * 4th letter to decide the job */
549  if( p.nMatch( "ADVE" ) )
550  {
551  /* save information relating to advection */
552  strcpy( save.chSave[save.nsave], "DYNa");
553  sprintf( save.chHeader[save.nsave],
554  "#advection depth\tHtot\tadCool\tadHeat\tdCoolHeatdT\t"
555  "Source[hyd][hyd]\tRate\tEnthalph\tadSpecEnthal\n" );
556  }
557  else
558  {
559  fprintf( ioQQQ, " I did not recognize a sub keyword on this SAVE DYNAMICS command.\n" );
560  fprintf( ioQQQ, " Sorry.\n" );
562  }
563  }
564 
565  else if( p.nMatch("ENTH") )
566  {
567  /* contributors to the total enthalpy */
568  strcpy( save.chSave[save.nsave], "ENTH" );
569  sprintf( save.chHeader[save.nsave],
570  "#depth\tTotal\tExcit\tIoniz\tBind\tKE\tther+PdV\tmag \n" );
571  }
572 
573  else if( p.nMatch("EXEC") && p.nMatch("TIME") )
574  {
575  /* output the execution time per zone */
576  strcpy( save.chSave[save.nsave], "XTIM" );
577  sprintf( save.chHeader[save.nsave],
578  "#zone\tdTime\tElapsed t\n" );
579  }
580 
581  else if( p.nMatch("FEII") || p.nMatch("FE II") )
582  {
583  /* something to do with FeII atom - several options
584  * FeII column densities */
585  if( p.nMatch("COLU") && p.nMatch("DENS") )
586  {
587  /* save FeII column density */
588  strcpy( save.chSave[save.nsave], "FENl" );
589 
590  /* file will give energy, statistical weight, and column density [cm-2] */
591  sprintf( save.chHeader[save.nsave],
592  "#FeII: energy\tstat wght\tcol den\n" );
593  }
594 
595  /* FeII continuum, only valid if large atom is set */
596  else if( p.nMatch("CONT") )
597  {
598  // save FeII continuum, options are total (default), inward,
599  // and outward
600  if( p.nMatch("INWA") )
601  {
602  // inward continuum
603  strcpy( save.chSave[save.nsave], "FEcI" );
604  //sprintf( save.chHeader[save.nsave],
605  // "#FeII inward: Wl(A)\tInt[erg cm-2 s-1]\n" );
606  }
607  else if( p.nMatch(" OUT") )
608  {
609  // outward continuum
610  strcpy( save.chSave[save.nsave], "FEcO" );
611  //sprintf( save.chHeader[save.nsave],
612  // "#FeII outward: Wl(A)\tInt[erg cm-2 s-1]\n" );
613  }
614  else
615  {
616  // total continuum
617  strcpy( save.chSave[save.nsave], "FEcT" );
618  //sprintf( save.chHeader[save.nsave],
619  // "#FeII total: Wl(A)\tInt[erg cm-2 s-1]\n" );
620  }
621 
622  // default units of spectral energy are Ryd, this can change to
623  // many other units
624  ChkUnits(p);
625 
626  // by default give numbers in two columns, row keyword says to
627  // write the numbers across as one long row
628  if( p.nMatch(" ROW") )
629  save.punarg[save.nsave][0] = 1;
630  else
631  // the default, two columns
632  save.punarg[save.nsave][0] = 2;
633  }
634 
635  else if( p.nMatch("DEPA") )
636  {
637  /* save out departure coefficient for the large FeII atom */
638  sprintf( save.chHeader[save.nsave],
639  "#FeII departure coefficient \n" );
640  /* optional keyword all means do all levels, if not then just do subset */
641  if( p.nMatch(" ALL") )
642  {
643  /* save all levels, calls routine FeIIPunDepart */
644  strcpy( save.chSave[save.nsave], "FE2D" );
645  }
646  else
647  {
648  /* save a very few selected levels, calls routine FeIIPunDepart */
649  strcpy( save.chSave[save.nsave], "FE2d" );
650  }
651  }
652 
653  else if( p.nMatch("LEVE") )
654  {
655  /* save level energies and statistical weights for large FeII atom */
656  sprintf( save.chHeader[save.nsave],
657  "#FeII energy(wn)\tstat weight\n" );
658  strcpy( save.chSave[save.nsave], "FE2l" );
659  }
660 
661  else if( p.nMatch("LINE") )
662  {
663  /* save FeII lines command
664  * three optional parameters, threshold for faintest
665  * line to print, lower and upper energy bounds */
666 
667  /* save intensities from large FeII atom */
668  strcpy( save.chSave[save.nsave], "FEli" );
669 
670  /* short keyword makes save half as big */
671  if( p.nMatch("SHOR") )
672  {
673  FeII.lgShortFe2 = true;
674  }
675  else
676  {
677  FeII.lgShortFe2 = false;
678  }
679 
680  /* first optional number changes the threshold of weakest line to print*/
681  /* fe2thresh is intensity relative to normalization line,
682  * normally Hbeta, and is set to zero in zero.c */
684  if( p.lgEOL() )
685  {
686  FeII.fe2thresh = 0.;
687  }
688 
689  /* it is a log if negative */
690  if( FeII.fe2thresh < 0. )
691  {
692  FeII.fe2thresh = (realnum)pow((realnum)10.f,FeII.fe2thresh);
693  }
694 
695  /* check for energy range (Rydberg) of lines to be saved,
696  * this is to limit size of output file */
697  FeII.fe2ener[0] = (realnum)p.FFmtRead();
698  if( p.lgEOL() )
699  {
700  FeII.fe2ener[0] = 0.;
701  }
702 
703  FeII.fe2ener[1] = (realnum)p.FFmtRead();
704  if( p.lgEOL() )
705  {
706  FeII.fe2ener[1] = 1e8;
707  }
708  /* if either is negative then both are logs */
709  if( FeII.fe2ener[0] < 0. || FeII.fe2ener[1] < 0. )
710  {
711  FeII.fe2ener[0] = (realnum)pow((realnum)10.f,FeII.fe2ener[0]);
712  FeII.fe2ener[1] = (realnum)pow((realnum)10.f,FeII.fe2ener[1]);
713  }
714 
715  /* entered in Ryd in above, convert to wavenumbers */
716  FeII.fe2ener[0] /= (realnum)WAVNRYD;
717  FeII.fe2ener[1] /= (realnum)WAVNRYD;
718 
719  /* these results are actually created by the FeIISaveLines routine
720  * that lives in the FeIILevelPops file */
721  sprintf( save.chHeader[save.nsave],
722  "#FeII ipHi\tipLo\tWL(A)\tlog I\tI/Inorm\t\tTau\n" );
723  }
724 
725  else if( p.nMatch("OPTI") && p.nMatch("DEPT") )
726  {
727  /* save optical depths for large FeII atom */
728  sprintf( save.chHeader[save.nsave],
729  "#FeII hi\tlow\twl(A)\ttau\n" );
730  strcpy( save.chSave[save.nsave], "FE2o" );
731  }
732 
733  else if( p.nMatch("POPU") )
734  {
735  /* save out level populations for the large FeII atom */
736  sprintf( save.chHeader[save.nsave],
737  "#FeII level populations [cm^-3]\n" );
738 
739  /* this is keyword RELATIVE that says to save relative to total Fe+,
740  * default is actual density (cm-3) */
741  if( p.nMatch("RELA") )
742  {
743  save.punarg[save.nsave][2] = 0.;
744  }
745  else
746  {
747  /* default is to save density (cm-3) */
748  save.punarg[save.nsave][2] = 1.;
749  }
750 
751  /* optional keyword all means do all levels, if not then just do subset */
752  if( p.nMatch(" ALL") )
753  {
754  /* save all levels, calls routine FeIIPunPop */
755  strcpy( save.chSave[save.nsave], "FE2P" );
756  save.punarg[save.nsave][0] = 0.;
758  }
759 
760  /* optional keyword range means read lower and upper bound, do these */
761  else if( p.nMatch("RANG") )
762  {
763  /* save range of levels, calls routine FeIIPunPop */
764  strcpy( save.chSave[save.nsave], "FE2P" );
765  save.punarg[save.nsave][0] = (realnum)p.FFmtRead();
766  save.punarg[save.nsave][1] = (realnum)p.FFmtRead();
767  if( p.lgEOL() || save.punarg[save.nsave][0] <0 ||
768  save.punarg[save.nsave][0]>= save.punarg[save.nsave][1] )
769  {
770  fprintf( ioQQQ, "There must be two numbers on this save "
771  "FeII populations range command.\n" );
772  fprintf( ioQQQ, "These give the lower and upper levels "
773  "for the range of FeII levels.\n" );
774  fprintf( ioQQQ, "The first, %g, must be less than the second, %g.\n",
775  save.punarg[save.nsave][0],
776  save.punarg[save.nsave][1]);
777  fprintf( ioQQQ, "Sorry.\n" );
779  }
780  }
781 
782  else
783  {
784  /* save a very few selected levels, calls routine FeIIPunPop */
785  strcpy( save.chSave[save.nsave], "FE2p" );
786  }
787  }
788  else
789  {
790  fprintf( ioQQQ, "There must be a second keyword on this SAVE FEII command.\n" );
791  fprintf( ioQQQ, "The ones I know about are COLUmn, CONTinuum, "
792  "DEPArture, LEVEls, LINE, OPTIcal DEPTh, and POPUlations.\n" );
793  fprintf( ioQQQ, "Sorry.\n" );
795  }
796  }
797 
798  /* the save continuum command, with many options,
799  * the first 3 char of the chSave flag will always be "CON"
800  * with the last indicating which one */
801  else if( p.nMatch("CONT") && !p.nMatch("XSPE") )
802  {
803  /* this flag is checked in PrtComment to generate a caution
804  * if continuum is saved but iterations not performed */
805  save.lgPunContinuum = true;
806 
807  /* check for keyword UNITS on line, then scan wavelength or energy units if present,
808  * units are copied into save.chConPunEnr */
809  ChkUnits(p);
810 
811  if( p.nMatch("BINS") )
812  {
813  /* continuum binning */
814  strcpy( save.chSave[save.nsave], "CONB" );
815 
816  sprintf( save.chHeader[save.nsave],
817  "#Continuum binning enrOrg/%s\tEnergy\twidth of cells\n",
819  }
820 
821  else if( p.nMatch("DIFF") )
822  {
823  /* diffuse continuum, the locally emitted lines and continuum */
824  strcpy( save.chSave[save.nsave], "COND" );
825 
826  /* by default gives lines and continuum separately only for
827  * last zone. The keyword ZONE says to give the total for every
828  * zone in one very low row */
829  if( p.nMatch("ZONE") )
830  {
831  sprintf( save.chHeader[save.nsave],
832  "#energy/%s then emission per zone\n",
834  save.punarg[save.nsave][0] = 2.;
835 
836  }
837  else
838  {
839  sprintf( save.chHeader[save.nsave],
840  "#energy/%s\tConEmitLocal\tDiffuseLineEmission\tTotal\n",
842  save.punarg[save.nsave][0] = 1.;
843  }
844  }
845 
846  else if( p.nMatch("EMIS") )
847  {
848  /* continuum volume emissivity and opacity as a function of radius */
849  strcpy( save.chSave[save.nsave], "CONS" );
850 
851  double num = p.FFmtRead();
852  if( p.lgEOL() )
853  p.NoNumb( "continuum emissivity frequency" );
855  if( save.emisfreq[save.nsave].Ryd() < rfield.emm ||
857  {
858  fprintf( ioQQQ, " The frequency is outside the Cloudy range\n Sorry.\n" );
860  }
861 
862  sprintf( save.chHeader[save.nsave],
863  "#Radius\tdepth\tnujnu\tkappa_abs\tkappa_sct @ %e Ryd\n",
864  save.emisfreq[save.nsave].Ryd() );
865  }
866 
867  else if( p.nMatch("EMIT") )
868  {
869  /* continuum emitted by cloud */
870  strcpy( save.chSave[save.nsave], "CONE" );
871 
872  sprintf( save.chHeader[save.nsave],
873  "#Energy/%s\treflec\toutward\ttotal\tline\tcont\n",
875  }
876 
877  else if( p.nMatch("FINE" ) )
878  {
879  rfield.lgSaveOpacityFine = true;
880  /* fine transmitted continuum cloud */
881  strcpy( save.chSave[save.nsave], "CONf" );
882 
883  sprintf( save.chHeader[save.nsave],
884  "#Energy/%s\tTransmitted\n",
886 
887  /* range option - important since so much data */
888  if( p.nMatch("RANGE") )
889  {
890  /* get lower and upper range, eventually must be in Ryd */
891  double Energy1 = p.FFmtRead();
892  double Energy2 = p.FFmtRead();
893  if( p.lgEOL() )
894  {
895  fprintf(ioQQQ,"There must be two numbers, the lower and upper energies in Ryd.\nSorry.\n");
897  }
898  if( p.nMatch("UNIT" ) )
899  {
900  // apply units to range option
901  const char *energyUnits = p.StandardEnergyUnit();
902  Energy unitChange;
903  unitChange.set(Energy1, energyUnits );
904  Energy1 = unitChange.Ryd();
905  unitChange.set(Energy2, energyUnits );
906  Energy2 = unitChange.Ryd();
907  }
908  /* get lower and upper rang in Ryd */
909  save.punarg[save.nsave][0] = (realnum)MIN2( Energy1 , Energy2 );
910  save.punarg[save.nsave][1] = (realnum)MAX2( Energy1 , Energy2 );
911  //fprintf(ioQQQ , "DEBUG units change fine %.3e %.3e\n" , save.punarg[save.nsave][0] ,
912  // save.punarg[save.nsave][1] );
913  //cdEXIT(EXIT_FAILURE);
914  }
915  else
916  {
917  /* these mean full energy range */
918  save.punarg[save.nsave][0] = 0.;
919  save.punarg[save.nsave][1] = 0.;
920  }
921  /* optional last parameter - how many points to bring together */
922  save.punarg[save.nsave][2] = (realnum)p.FFmtRead();
923 
924  /* default is to bring together ten */
925  if( p.lgEOL() )
926  save.punarg[save.nsave][2] = 10;
927 
928  if( save.punarg[save.nsave][2] < 1 )
929  {
930  fprintf(ioQQQ,"The number of fine opacities to skip must be > 0 \nSorry.\n");
932  }
933  }
934 
935  else if( p.nMatch("GRAI") )
936  {
937  /* save grain continuum in optically thin limit */
938  strcpy( save.chSave[save.nsave], "CONG" );
939 
940  sprintf( save.chHeader[save.nsave],
941  "#energy\tgraphite\trest\ttotal\n" );
942  }
943 
944  else if( p.nMatch("INCI") )
945  {
946  /* incident continuum */
947  strcpy( save.chSave[save.nsave], "CONC" );
948 
949  sprintf( save.chHeader[save.nsave],
950  "#Incident Continuum, Enr\tnFn \n" );
951  }
952 
953  else if( p.nMatch("INTE") )
954  {
955  /* continuum interactions */
956  strcpy( save.chSave[save.nsave], "CONi" );
957 
958  sprintf( save.chHeader[save.nsave],
959  "#Continuum interactions, inc, otslin. otscon, ConInterOut, outlin \n" );
960  /* this is option for lowest energy, if nothing then zero */
961  save.punarg[save.nsave][0] = (realnum)p.FFmtRead();
962  }
963 
964  else if( p.nMatch("IONI") )
965  {
966  /* save ionizing continuum*/
967  strcpy( save.chSave[save.nsave], "CONI" );
968 
969  /* this is option for lowest energy, if nothing then zero */
970  save.punarg[save.nsave][0] = (realnum)p.FFmtRead();
971 
972  /* this is option for smallest interaction to save, def is 1 percent */
973  save.punarg[save.nsave][1] = (realnum)p.FFmtRead();
974  if( p.lgEOL() )
975  save.punarg[save.nsave][1] = 0.01f;
976 
977  /* "every" option to save this on every zone -
978  * not present then only last zone is saved */
979  if( p.nMatch("EVER" ) )
980  {
981  /* save every zone */
982  save.lgSaveEveryZone[save.nsave] = true;
984  }
985  else
986  {
987  /* only save last zone */
988  save.lgSaveEveryZone[save.nsave] = false;
990  }
991 
992  /* put the header at the top of the file */
993  sprintf( save.chHeader[save.nsave],
994  "#cell(on C scale)\tnu\tflux\tflx*cs\tFinc\totsl\totsc\toutlin\toutcon\trate/tot\tintegral\tline\tcont\n" );
995  }
996 #ifdef USE_NLTE7
997  else if( p.nMatch("NLTE") )
998  {
999  /* continuum emitted by cloud */
1000  strcpy( save.chSave[save.nsave], "CONl" );
1001 
1002  sprintf( save.chHeader[save.nsave],
1003  " spectrum1 NeXY6 XNUMX\n");
1004  }
1005 #endif
1006 
1007  else if( p.nMatch("OUTW") )
1008  {
1009  /* outward only continuum */
1010  if( p.nMatch("LOCA") )
1011  {
1012  strcpy( save.chSave[save.nsave], "CONo" );
1013  sprintf( save.chHeader[save.nsave],
1014  "#Local Out ConInterOut+line SvOt*opc pass*opc\n" );
1015  }
1016  else
1017  {
1018  strcpy( save.chSave[save.nsave], "CONO" );
1019  sprintf( save.chHeader[save.nsave],
1020  "#Out Con OutIncid OutConD OutLinD OutConS\n" );
1021  }
1022  }
1023 
1024  else if( p.nMatch("TRAN") )
1025  {
1026  /* transmitted continuum */
1027  strcpy( save.chSave[save.nsave], "CONT" );
1028 
1029  sprintf( save.chHeader[save.nsave],
1030  "#ener\tTran Contin\ttrn coef\n" );
1031  }
1032 
1033  else if( p.nMatch(" TWO") )
1034  {
1035  /* total two photon continua rfield.TotDiff2Pht */
1036  strcpy( save.chSave[save.nsave], "CON2" );
1037 
1038  sprintf( save.chHeader[save.nsave],
1039  "#energy\t n_nu\tnuF_nu \n" );
1040  }
1041 
1042  else if( p.nMatch(" RAW") )
1043  {
1044  /* "raw" continua */
1045  strcpy( save.chSave[save.nsave], "CORA" );
1046 
1047  sprintf( save.chHeader[save.nsave],
1048  "#Raw Con anu\tflux\totslin\totscon\tConRefIncid\tConEmitReflec\tConInterOut\toutlin\tConEmitOut\tline\tcont\tnLines\n" );
1049  }
1050 
1051  else if( p.nMatch("REFL") )
1052  {
1053  /* reflected continuum */
1054  strcpy( save.chSave[save.nsave], "CONR" );
1055 
1056  sprintf( save.chHeader[save.nsave],
1057  "#Reflected\tcont\tline\ttotal\talbedo\tConID\n" );
1058  }
1059 
1060  else
1061  {
1062  /* this is the usual save continuum command,
1063  * ipType is index for continuum array to set either
1064  * iteration or cumulative output */
1065  int ipType = 0;
1066  if( p.nMatch( "CUMU" ) )
1067  ipType = 1;
1068  save.punarg[save.nsave][0] = (realnum)ipType;
1069  strcpy( save.chSave[save.nsave], "CON " );
1070  char chHold[100];
1071  strcpy( chHold, "#Cont " );
1072  if( ipType > 0 )
1073  strcpy( chHold , "#Cumul " );
1074  sprintf( save.chHeader[save.nsave],
1075  "%s nu\tincident\ttrans\tDiffOut\tnet trans\treflc\ttotal\treflin\toutlin\tlineID\tcont\tnLine\n" ,
1076  chHold );
1077 
1078  /* >>chng 06 apr 03, add "every" option to save this on every zone -
1079  * if every is not present then only last zone is saved */
1080  if( p.nMatch("EVER" ) )
1081  {
1082  /* save every zone */
1083  save.lgSaveEveryZone[save.nsave] = true;
1084  /* option to say how many to skip */
1085  save.nSaveEveryZone[save.nsave] = (long)p.FFmtRead();
1086  if( p.lgEOL() )
1088  }
1089  else
1090  {
1091  /* only save last zone */
1092  save.lgSaveEveryZone[save.nsave] = false;
1094  }
1095  }
1096  }
1097 
1098  /* save information about convergence of this model
1099  * reason - why it did not converge an iteration
1100  * error - zone by zone display of various convergence errors */
1101  else if( p.nMatch("CONV") )
1102  {
1103  if( p.nMatch("REAS") )
1104  {
1105  /* this does not count as a save option (really) */
1106  save.lgPunConv = true;
1107  /* this is done below */
1108  strcpy( save.chSave[save.nsave], "" );
1109  save.lgRealSave[save.nsave] = false;
1110  }
1111  else if( p.nMatch("ERRO") )
1112  {
1113  /* save zone by zone errors in pressure, electron density, and heating-cooling */
1114  /* convergence error */
1115  strcpy( save.chSave[save.nsave], "CNVE" );
1116  sprintf( save.chHeader[save.nsave],
1117  "#depth\tnPres2Ioniz\tP(cur)\tP%%error\tNE(cor)\tNE(cur)\tNE%%error\tHeat\tCool\tHC%%error\n" );
1118  }
1119  else if( p.nMatch("BASE") )
1120  {
1121  /* save converged quantities in Converge base for each pass through
1122  * solvers - not one pass per zone */
1123  strcpy( save.chSave[save.nsave], "CNVB" );
1124  strcpy( save.chSave[save.nsave], "" );
1125  save.lgRealSave[save.nsave] = false;
1126  }
1127  else
1128  {
1129  fprintf( ioQQQ, "There must be a second keyword on this command.\n" );
1130  fprintf( ioQQQ, "The ones I know about are REASON, ERROR, and BASE.\n" );
1131  fprintf( ioQQQ, "Sorry.\n" );
1133  }
1134  }
1135 
1136  else if( p.nMatch(" DR ") )
1137  {
1138  /* first occurrence of save dr to follow choice in change of zoning */
1139  save.lgDROn = true;
1140  strcpy( save.chSave[save.nsave], "" );
1141  save.lgRealSave[save.nsave] = false;
1142  }
1143 
1144  else if( p.nMatch("ELEM") && !p.nMatch("GAMMA") && !p.nMatch("COOL") ) // do not trip on SAVE COOLING EACH ELEMENT
1145  {
1146  /* option to save ionization structure of some element
1147  * will give each stage of ionization, vs depth */
1148  strcpy( save.chSave[save.nsave], "ELEM" );
1149 
1150  /* this returns element number on c scale */
1151  /* >>chng 04 nov 23, had converted to f scale, leave on c */
1152  nelem = p.GetElem();
1153  if( nelem < 0 || nelem >= LIMELM )
1154  {
1155  fprintf( ioQQQ, " I could not recognize a valid element name on this line.\n" );
1156  fprintf( ioQQQ, " Please check your input script. Bailing out...\n" );
1158  }
1159 
1160  /* this is the atomic number on the c scale */
1161  save.punarg[save.nsave][0] = (realnum)nelem;
1162 
1163  /* >>chng 04 nov 24, add DENSE option to print density rather than fraction */
1164  save.punarg[save.nsave][1] = 0;
1165  if( p.nMatch("DENS") )
1166  save.punarg[save.nsave][1] = 1.;
1167 
1168  /* start printing header line - first will be the depth in cm */
1169  sprintf( save.chHeader[save.nsave], "#depth");
1170 
1171  /* next come the nelem+1 ion stages */
1172  for(i=0; i<=nelem+1;++i )
1173  {
1174  sprintf( chTemp,
1175  "\t%2s%2li", elementnames.chElementSym[nelem],i+1);
1176  strcat( save.chHeader[save.nsave], chTemp );
1177  }
1178 
1179  /* finally some fine structure or molecular populations */
1180  /* >>chng 04 nov 23, add fs pops of C, O
1181  * >>chng 04 nov 25, add molecules */
1182  if( nelem==ipHYDROGEN )
1183  {
1184  sprintf( chTemp, "\tH2");
1185  strcat( save.chHeader[save.nsave], chTemp );
1186  }
1187  if( nelem==ipCARBON )
1188  {
1189  sprintf( chTemp, "\tC1\tC1*\tC1**\tC2\tC2*\tCO");
1190  strcat( save.chHeader[save.nsave], chTemp );
1191  }
1192  else if( nelem==ipOXYGEN )
1193  {
1194  sprintf( chTemp, "\tO1\tO1*\tO1**");
1195  strcat( save.chHeader[save.nsave], chTemp );
1196  }
1197 
1198  /* finally the new line */
1199  strcat( save.chHeader[save.nsave], "\n");
1200  }
1201 
1202  else if( p.nMatch("FITS") )
1203  {
1204 
1205 #ifdef FLT_IS_DBL
1206  fprintf( ioQQQ, "Saving FITS files is not currently supported in double precision.\n" );
1207  fprintf( ioQQQ, "Please recompile without the FLT_IS_DBL option.\n" );
1208  fprintf( ioQQQ, "Sorry.\n" );
1210 #else
1211  /* say that this is a FITS file output */
1212  save.lgFITS[save.nsave] = true;
1213  /* concatenating files in a grid run would be illegal FITS */
1215  save.lgPunLstIter[save.nsave] = true;
1217 
1218  strcpy( save.chSave[save.nsave], "FITS" );
1219 #endif
1220 
1221  }
1222 
1223  else if( p.nMatch("FRED") )
1224  {
1225  /* save out some stuff for Fred's dynamics project */
1226  sprintf( save.chHeader[save.nsave],
1227  "#Radius\tDepth\tVelocity(km/s)\tdvdr(cm/s)\thden\teden\tTemperature\tRadAccel line\tRadAccel con\t"
1228  "Force multiplier\ta(e thin)\t"
1229  "HI\tHII\tHeI\tHeII\tHeIII\tC2\tC3\tC4\tO1\t"
1230  "O2\tO3\tO4\tO5\tO6\tO7\tO8\t"
1231  "HI\tHII\tHeI\tHeII\tHeIII\tC2\tC3\tC4\tO1\t"
1232  "O2\tO3\tO4\tO5\tO6\tO7\tO8\tMg2\tMg2\tOVI(1034) TauIn\tTauCon\n");
1233 
1234  strcpy( save.chSave[save.nsave], "FRED" );
1235  }
1236 
1237  else if( p.nMatch("GAMM") )
1238  {
1239  /* save all photoionization rates for all subshells */
1240  sprintf( save.chHeader[save.nsave],
1241  "#Photoionization rates \n" );
1242  if( p.nMatch("ELEMENT") )
1243  {
1244  /* element keyword, find element name and stage of ionization,
1245  * will print photoionization rates for valence of that element */
1246  strcpy( save.chSave[save.nsave], "GAMe" );
1247 
1248  /* this returns element number on c scale */
1249  nelem = p.GetElem();
1250  /* this is the atomic number on the C scale */
1251  save.punarg[save.nsave][0] = (realnum)nelem;
1252 
1253  /* this will become the ionization stage on C scale */
1254  save.punarg[save.nsave][1] = (realnum)p.FFmtRead() - 1;
1255  if( p.lgEOL() )
1256  p.NoNumb("element ionization stage" );
1257  if( save.punarg[save.nsave][1]<0 || save.punarg[save.nsave][1]> nelem+1 )
1258  {
1259  fprintf(ioQQQ,"Bad ionization stage - please check Hazy.\nSorry.\n");
1261  }
1262  }
1263  else
1264  {
1265  /* no element - so make table of all rates */
1266  strcpy( save.chSave[save.nsave], "GAMt" );
1267  }
1268 
1269  }
1270  else if( p.nMatch("GRAI") )
1271  {
1272  /* save grain ... options */
1273  if( p.nMatch("OPAC") )
1274  {
1275  /* check for keyword UNITS on line, then scan wavelength or energy units,
1276  * sets save.chConPunEnr*/
1277  ChkUnits(p);
1278 
1279  strcpy( save.chSave[save.nsave], "DUSO" );
1280  /* save grain opacity command in twice, here and above in opacity */
1281  sprintf( save.chHeader[save.nsave],
1282  "#grain\tnu/%s\tabs+scat*(1-g)\tabs\tscat*(1-g)\tscat\tscat*(1-g)/[abs+scat*(1-g)]\n",
1284  }
1285  else if( p.nMatch("ABUN") )
1286  {
1287  /* save grain abundance */
1288  strcpy( save.chSave[save.nsave], "DUSA" );
1289  sprintf( save.chHeader[save.nsave],
1290  "#grain\tdepth\tabundance (g/cm^3)\n" );
1291  }
1292  else if( p.nMatch("D/G ") )
1293  {
1294  /* save grain dust/gas mass ratio */
1295  strcpy( save.chSave[save.nsave], "DUSD" );
1296  sprintf( save.chHeader[save.nsave],
1297  "#grain\tdepth\tdust/gas mass ratio\n" );
1298  }
1299  else if( p.nMatch("PHYS") )
1300  {
1301  /* save grain physical conditions */
1302  strcpy( save.chSave[save.nsave], "DUSP" );
1303  sprintf( save.chHeader[save.nsave],
1304  "#grain\tdepth\tpotential\n" );
1305  }
1306  else if( p.nMatch(" QS ") )
1307  {
1308  strcpy( save.chSave[save.nsave], "DUSQ" );
1309  sprintf( save.chHeader[save.nsave],
1310  "#grain\tnu\tQ_abs\tQ_scat\n" );
1311  }
1312  else if( p.nMatch("TEMP") )
1313  {
1314  /* save temperatures of each grain species */
1315  strcpy( save.chSave[save.nsave], "DUST" );
1316  /* cannot save grain labels since they are not known yet */
1317  sprintf( save.chHeader[save.nsave],
1318  "#grain temperature\n" );
1319  }
1320  else if( p.nMatch("DRIF") )
1321  {
1322  /* save drift velocity of each grain species */
1323  strcpy( save.chSave[save.nsave], "DUSV" );
1324  /* cannot save grain labels since they are not known yet */
1325  sprintf( save.chHeader[save.nsave],
1326  "#grain drift velocity\n" );
1327  }
1328  else if( p.nMatch("EXTI") )
1329  {
1330  /* save grain extinction */
1331  strcpy( save.chSave[save.nsave], "DUSE" );
1332  /* cannot save grain labels since they are not known yet */
1333  sprintf( save.chHeader[save.nsave],
1334  "#depth\tA_V(extended)\tA_V(point)\n" );
1335  }
1336  else if( p.nMatch("CHAR") )
1337  {
1338  /* save charge per grain (# elec/grain) for each grain species */
1339  strcpy( save.chSave[save.nsave], "DUSC" );
1340  /* cannot save grain labels since they are not known yet */
1341  sprintf( save.chHeader[save.nsave],
1342  "#grain charge\n" );
1343  }
1344  else if( p.nMatch("HEAT") )
1345  {
1346  /* save heating due to each grain species */
1347  strcpy( save.chSave[save.nsave], "DUSH" );
1348  /* cannot save grain labels since they are not known yet */
1349  sprintf( save.chHeader[save.nsave],
1350  "#grain heating\n" );
1351  }
1352  else if( p.nMatch("POTE") )
1353  {
1354  /* save floating potential of each grain species */
1355  strcpy( save.chSave[save.nsave], "DUSP" );
1356  /* cannot save grain labels since they are not known yet */
1357  sprintf( save.chHeader[save.nsave],
1358  "#grain\tdepth\tpotential\n" );
1359  }
1360  else if( p.nMatch("H2RA") )
1361  {
1362  /* save grain H2rate - H2 formation rate for each type of grains */
1363  strcpy( save.chSave[save.nsave], "DUSR" );
1364  /* cannot save grain labels since they are not known yet */
1365  sprintf( save.chHeader[save.nsave],
1366  "#grain H2 formation rates\n" );
1367  }
1368  else
1369  {
1370  fprintf( ioQQQ, "There must be a second key on this GRAIN command; The options I know about follow (required key in CAPS):\n");
1371  fprintf( ioQQQ, "OPACity, ABUNdance, D/G mass ratio, PHYSical conditions, QS , TEMPerature, DRIFt velocity, EXTInction, CHARge, HEATing, POTEntial, H2RAtes\nSorry.\n" );
1373  }
1374  }
1375 
1376  else if( p.nMatch("GAUN") )
1377  {
1378  strcpy( save.chSave[save.nsave], "GAUN" );
1379  sprintf( save.chHeader[save.nsave],
1380  "#Gaunt factors.\n" );
1381  }
1382  else if( p.nMatch("GRID") )
1383  {
1384  strcpy( save.chSave[save.nsave], "GRID" );
1385  /* automatically generate no hash option */
1386  save.lgHashEndIter[save.nsave] = false;
1387  }
1388  else if( p.nMatch( "HIST" ) )
1389  {
1390  /* save pressure history of current zone */
1391  if( p.nMatch( "PRES") )
1392  {
1393  /* save pressure history - density - pressure for this zone */
1394  strcpy( save.chSave[save.nsave], "HISp" );
1395  sprintf( save.chHeader[save.nsave],
1396  "#iter zon\tdensity\tpres cur\tpres error\n" );
1397  }
1398  /* save temperature history of current zone */
1399  else if( p.nMatch( "TEMP" ) )
1400  {
1401  /* save pressure history - density - pressure for this zone */
1402  strcpy( save.chSave[save.nsave], "HISt" );
1403  sprintf( save.chHeader[save.nsave],
1404  "#iter zon\ttemperature\theating\tcooling\n" );
1405  }
1406  }
1407 
1408  else if( p.nMatch("HTWO") )
1409  {
1410  fprintf(ioQQQ," Sorry, this command has been replaced with the "
1411  "SAVE H2 CREATION and SAVE H2 DESTRUCTION commands.\n");
1413  }
1414 
1415  /* QHEAT has to come before HEAT... */
1416  else if( p.nMatch("QHEA") )
1417  {
1418  /* this is just a dummy clause, do the work below after parsing is over.
1419  * this is a no-nothing, picked up to stop optimizer */
1420  ((void)0);
1421  }
1422 
1423  else if( p.nMatch("HEAT") )
1424  {
1425  /* save heating */
1426  strcpy( save.chSave[save.nsave], "HEAT" );
1427  /*>>chng 06 jun 06, revise to be same as save cooling */
1428  sprintf( save.chHeader[save.nsave],
1429  "#depth cm\tTemp K\tHtot erg/cm3/s\tCtot erg/cm3/s\theat fracs\n" );
1430  }
1431 
1432  else if( p.nMatch("HELI") &&!( p.nMatch("IONI")))
1433  {
1434  /* save helium & helium-like iso sequence, but not save helium ionization rate
1435  * save helium line wavelengths */
1436  if( p.nMatch("LINE") && p.nMatch("WAVE") )
1437  {
1438  strcpy( save.chSave[save.nsave], "HELW" );
1439  sprintf( save.chHeader[save.nsave],
1440  "#wavelengths of lines from He-like ions\n" );
1441  }
1442  else
1443  {
1444  fprintf( ioQQQ, "save helium has options: LINE WAVElength.\nSorry.\n" );
1446  /* no key */
1447  }
1448  }
1449 
1450  else if( p.nMatch("HUMM") )
1451  {
1452  strcpy( save.chSave[save.nsave], "HUMM" );
1453  sprintf( save.chHeader[save.nsave],
1454  "#input to DHs routine.\n" );
1455  }
1456 
1457  else if( p.nMatch("HYDR") )
1458  {
1459  /* save hydrogen physical conditions */
1460  if( p.nMatch("COND") )
1461  {
1462  strcpy( save.chSave[save.nsave], "HYDc" );
1463  sprintf( save.chHeader[save.nsave],
1464  "#depth\tTe\tHDEN\tEDEN\tHI/H\tHII/H\tH2/H\tH2+/H\tH3+/H\tH-/H\n" );
1465  /* save hydrogen ionization */
1466  }
1467 
1468  /* save information on 21 cm excitation processes - accept either keyword 21cm or 21 cm */
1469  else if( p.nMatch("21 CM") ||p.nMatch("21CM"))
1470  {
1471  /* save information about 21 cm line */
1472  strcpy( save.chSave[save.nsave], "21CM" );
1473  sprintf( save.chHeader[save.nsave],
1474  "#depth\tT(spin)\tT(kin)\tT(Lya/21cm)\tnLo\tnHi\tOccLya\ttau(21cm)"
1475  "\ttau(Lya)\topac(21 cm)\tn/Ts\ttau(21)\tTex(Lya)\tN(H0)/Tspin"
1476  "\tSum_F0\tSum_F1\tSum_T21\n" );
1477  }
1478 
1479  else if( p.nMatch("IONI") )
1480  {
1481  /* save hydrogen ionization */
1482  strcpy( save.chSave[save.nsave], "HYDi" );
1483  sprintf( save.chHeader[save.nsave],
1484  "#hion\tzn\tgam1\tcoll ion1\tRecTot\tHRecCaB\thii/hi\tSim hii/hi"
1485  "\time_Hrecom_long(esc)\tdec2grd\texc pht\texc col\trec eff\tsec ion\n" );
1486  }
1487  else if( p.nMatch("POPU") )
1488  {
1489  /* save hydrogen populations */
1490  strcpy( save.chSave[save.nsave], "HYDp" );
1491  sprintf( save.chHeader[save.nsave],
1492  "#depth\tn(H0)\tn(H+)\tn(1s)\tn(2s)\tn(2p)\tetc\n" );
1493  }
1494  else if( p.nMatch("LINE") )
1495  {
1496  /* save hydrogen lines
1497  * hydrogen line intensities and optical depths */
1498  strcpy( save.chSave[save.nsave], "HYDl" );
1499  sprintf( save.chHeader[save.nsave],
1500  "#nHi\tlHi\tnLo\tlLo\tE(ryd)\ttau\n" );
1501  }
1502  else if( p.nMatch(" LYA") )
1503  {
1504  /* save hydrogen Lya some details about Lyman alpha */
1505  strcpy( save.chSave[save.nsave], "HYDL" );
1506  sprintf( save.chHeader[save.nsave],
1507  "#depth\tTauIn\tTauTot\tn(2p)/n(1s)\tTexc\tTe\tTex/T\tPesc\tPdes\tpump\topacity\talbedo\n" );
1508  }
1509  else
1510  {
1511  fprintf( ioQQQ, "Save hydrogen has options: CONDitions, 21 CM, LINE, POPUlations, and IONIzation.\nSorry.\n" );
1513  }
1514  }
1515 
1516  else if( p.nMatch("IONI") )
1517  {
1518  if( p.nMatch("RATE") )
1519  {
1520  /* save ionization rates, search for the name of an element */
1521  if( (nelem = p.GetElem() ) < 0 )
1522  {
1523  fprintf( ioQQQ, "There must be an element name on the ionization rates command. Sorry.\n" );
1525  }
1526  save.punarg[save.nsave][0] = (realnum)nelem;
1527  strcpy( save.chSave[save.nsave], "IONR" );
1528  sprintf( save.chHeader[save.nsave],
1529  "#%s depth\teden\tdynamics.Rate\tabund\tTotIonize\tTotRecom\tSource\t ... \n",
1530  elementnames.chElementSym[nelem]);
1531  }
1532  else
1533  {
1534  /* save table giving ionization means */
1535  strcpy( save.chSave[save.nsave], "IONI" );
1536  sprintf( save.chHeader[save.nsave],
1537  "#Mean ionization distribution\n" );
1538  }
1539  }
1540 
1541  else if( p.nMatch(" IP ") )
1542  {
1543  strcpy( save.chSave[save.nsave], " IP " );
1544  sprintf( save.chHeader[save.nsave],
1545  "#ionization potentials, valence shell\n" );
1546  }
1547 
1548  else if( p.nMatch("LEID") )
1549  {
1550  if( p.nMatch( "LINE" ) )
1551  {
1552  /* save Leiden lines
1553  * final intensities of the Leiden PDR models */
1554  strcpy( save.chSave[save.nsave], "LEIL" );
1555  sprintf( save.chHeader[save.nsave], "#ion\twl\tInt\trel int\n");
1556  }
1557  else
1558  {
1559  /* save Leiden structure
1560  * structure of the Leiden PDR models */
1561  strcpy( save.chSave[save.nsave], "LEIS" );
1562  sprintf( save.chHeader[save.nsave],
1563  /* 1-17 */
1564  "#Leid depth\tA_V(extentd)\tA_V(point)\tTe\tH0\tH2\tCo\tC+\tOo\tCO\tO2\tCH\tOH\te\tHe+\tH+\tH3+\t"
1565  /* 18 - 30 */
1566  "N(H0)\tN(H2)\tN(Co)\tN(C+)\tN(Oo)\tN(CO)\tN(O2)\tN(CH)\tN(OH)\tN(e)\tN(He+)\tN(H+)\tN(H3+)\t"
1567  /* 31 - 32 */
1568  "H2(Sol)\tH2(FrmGrn)\tH2(photodiss)\t"
1569  /* 33 - 46*/
1570  "G0(DB96)\trate(CO)\trate(C)\theat\tcool\tGrnP\tGr-Gas-Cool\tGr-Gas-Heat\tCOds\tH2dH\tH2vH\tChaT\tCR H\tMgI\tSI\t"
1571  "Si\tFe\tNa\tAl\tC\tC610\tC370\tC157\tC63\tC146\n" );
1572  }
1573  }
1574 
1575 
1576  // save results for NLTE series of plasma comparison meetings,
1577  // specifically NLTE7 2011 Dec
1578  else if( p.nMatch("NLTE") )
1579  {
1580 # ifdef USE_NLTE7
1581  strcpy( save.chSave[save.nsave], "NLTE" );
1582 # else
1583  fprintf(ioQQQ," PROBLEM You must enable the USE_NLTE7 macro at compile-time to use this command.\n");
1584  fprintf(ioQQQ," To do so, add EXTRA=\"-DUSE_NLTE7\" to the end of the make command.\n");
1585  fprintf(ioQQQ," An example for a quad core machine:\n make -j 4 EXTRA=\"-DUSE_NLTE7\" \n");
1586  fprintf(ioQQQ," in the sys_XXX folder that you want to use.\n\n\n");
1588 # endif
1589  }
1590 
1591  /* FE2NRG, FE2TP, and FE2COLL write the internal Fe II data into Stout format */
1592  else if (p.nMatch("FE2NRG"))
1593  {
1594  strcpy( save.chSave[save.nsave], "LY1" );
1595  }
1596 
1597  else if (p.nMatch("FE2TP"))
1598  {
1599  strcpy( save.chSave[save.nsave], "LY2" );
1600  }
1601 
1602  else if (p.nMatch("FE2COLL"))
1603  {
1604  strcpy( save.chSave[save.nsave], "LY3" );
1605  }
1606 
1607  else if( (p.nMatch("LINE") && p.nMatch("LIST")) || p.nMatch("LINELIST") )
1608  {
1609  /* save line list "output file" "Line List file" */
1610  strcpy( save.chSave[save.nsave], "LLST" );
1611 
1612  /*
1613  * we parsed off the second file name at start of this routine
1614  * check if file was found, use it if it was, else abort
1615  */
1616  if( !lgSecondFilename )
1617  {
1618  fprintf(ioQQQ , "There must be a second file name between "
1619  "double quotes on the SAVE LINE LIST command. This second"
1620  " file contains the input line list. I did not find it.\nSorry.\n");
1622  }
1623 
1624  /* actually get the lines, and malloc the space in the arrays
1625  * cdGetLineList will look on path */
1626  if( save.ipPnunit[save.nsave] == NULL )
1627  {
1628  /* make sure we free any allocated space from a previous call */
1630 
1631  save.nLineList[save.nsave] = cdGetLineList(chSecondFilename,
1634 
1635  if( save.nLineList[save.nsave] < 0 )
1636  {
1637  fprintf(ioQQQ,"DISASTER could not open SAVE LINE LIST file %s \n",
1638  chSecondFilename );
1640  }
1641  }
1642 
1643  // check whether intrinsic or emergent line emissivity
1644  save.lgEmergent[save.nsave] = false;
1645  if( p.nMatch("EMER") )
1646  save.lgEmergent[save.nsave] = true;
1647 
1648  // check whether cumulative or specific line emission
1649  save.lgCumulative[save.nsave] = false;
1650  if( p.nMatch("CUMU") )
1651  save.lgCumulative[save.nsave] = true;
1652 
1653  /* ratio option, in which pairs of lines form ratios, first over
1654  * second */
1655  if( p.nMatch("RATI") )
1656  {
1657  save.lgLineListRatio[save.nsave] = true;
1658  if( save.nLineList[save.nsave]%2 )
1659  {
1660  /* odd number of lines - cannot take ratio */
1661  fprintf(ioQQQ , "There must be an even number of lines to"
1662  " take ratios of lines. There were %li, an odd number."
1663  "\nSorry.\n", save.nLineList[save.nsave]);
1665  }
1666  }
1667  else
1668  {
1669  /* no ratio */
1670  save.lgLineListRatio[save.nsave] = false;
1671  }
1672 
1673  /* keyword absolute says to do absolute rather than relative intensities
1674  * relative intensities are the default */
1675  if( p.nMatch("ABSO") )
1676  {
1677  save.punarg[save.nsave][0] = 1;
1678  }
1679  else
1680  {
1681  save.punarg[save.nsave][0] = 0;
1682  }
1683 
1684  // check whether column or row (default)
1685  if( p.nMatch("COLUMN") )
1686  {
1687  save.punarg[save.nsave][1] = 1;
1688  }
1689  else
1690  {
1691  save.punarg[save.nsave][1] = 0;
1692  }
1693 
1694  /* give header line */
1695  sprintf( save.chHeader[save.nsave], "#lineslist" );
1696  // do header now if reporting rows of lines
1697  if( !save.punarg[save.nsave][1] )
1698  {
1699  for( long int j=0; j<save.nLineList[save.nsave]; ++j )
1700  {
1701  /* if taking ratio then put div sign between pairs */
1702  if( save.lgLineListRatio[save.nsave] && is_odd(j) )
1703  strcat( save.chHeader[save.nsave] , "/" );
1704  else
1705  strcat( save.chHeader[save.nsave] , "\t" );
1706  sprintf( chTemp, "%s ", save.chLineListLabel[save.nsave][j] );
1707  strcat( save.chHeader[save.nsave], chTemp );
1708  sprt_wl( chTemp, save.wlLineList[save.nsave][j] );
1709  strcat( save.chHeader[save.nsave], chTemp );
1710  }
1711  }
1712  strcat( save.chHeader[save.nsave], "\n" );
1713  }
1714 
1715  else if( p.nMatch("LINE") && !p.nMatch("XSPE") && !p.nMatch("NEAR"))
1716  {
1717  /* save line options -
1718  * this is not save xspec lines and not linear option
1719  * check for keyword UNITS on line, then scan wavelength or energy units,
1720  * sets save.chConPunEnr*/
1721  ChkUnits(p);
1722 
1723  /* save line emissivity, line intensity, line array,
1724  * and line data */
1725  if( p.nMatch("STRU") )
1726  {
1727  fprintf(ioQQQ," The SAVE LINES STRUCTURE command is now SAVE LINES "
1728  "EMISSIVITY.\n Sorry.\n\n");
1730  }
1731 
1732  else if( p.nMatch("PRES") )
1733  {
1734  /* save contributors to line pressure */
1735  strcpy( save.chSave[save.nsave], "PREL" );
1736  sprintf( save.chHeader[save.nsave],
1737  "#P depth\tPtot\tPline/Ptot\tcontributors to line pressure\n" );
1738  }
1739 
1740  else if( p.nMatch("EMIS") )
1741  {
1742  /* this used to be the save lines structure command, is now
1743  * the save lines emissivity command
1744  * give line emissivity vs depth */
1745  // check whether intrinsic or emergent line emissivity
1746  save.lgEmergent[save.nsave] = false;
1747  if( p.nMatch("EMER") )
1748  save.lgEmergent[save.nsave] = true;
1749  strcpy( save.chSave[save.nsave], "LINS" );
1750  sprintf( save.chHeader[save.nsave],
1751  "#");
1752  /* read in the list of lines to examine */
1753  parse_save_line(p,false, chTemp );
1754  strcat( save.chHeader[save.nsave], chTemp );
1755  }
1756 
1757  else if( p.nMatch(" RT " ) )
1758  {
1759  /* save line RT */
1760  strcpy( save.chSave[save.nsave], "LINR" );
1761  /* save some details needed for line radiative transfer
1762  * routine in save_line.cpp */
1763  Parse_Save_Line_RT(p);
1764  }
1765 
1766  else if( p.nMatch("CUMU") )
1767  {
1768  bool lgEOL;
1769  /* save lines cumulative
1770  * this will be integrated line intensity, function of depth */
1771  strcpy( save.chSave[save.nsave], "LINC" );
1772  // option for intrinsic (default) or emergent
1773  save.lgEmergent[save.nsave] = false;
1774  if( p.nMatch("EMER") )
1775  save.lgEmergent[save.nsave] = true;
1776  /* option for either relative intensity or abs luminosity */
1777  if( p.nMatch("RELA") )
1778  {
1779  lgEOL = true;
1780  sprintf( save.chHeader[save.nsave], "#" );
1781  }
1782  else
1783  {
1784  sprintf( save.chHeader[save.nsave], "#" );
1785  lgEOL = false;
1786  }
1787  /* read in the list of lines to examine */
1788  parse_save_line(p, lgEOL, chTemp );
1789  strcat( save.chHeader[save.nsave], chTemp );
1790  }
1791 
1792  else if( p.nMatch("DATA") )
1793  {
1794  /* save line data, done in SaveLineData */
1795 
1796  /* the default will be to make wavelengths like in the printout, called labels,
1797  * if units appears then other units will be used instead */
1798  save.chConPunEnr[save.nsave] = "labl";
1799 
1800  /* check for keyword UNITS on line, then scan wavelength or energy units if present,
1801  * units are copied into save.chConPunEnr */
1802  if( p.nMatch("UNIT") )
1803  ChkUnits(p);
1804  strcpy( save.chSave[save.nsave], "LIND" );
1805  sprintf( save.chHeader[save.nsave],
1806  "#Emission line data.\n" );
1807  }
1808 
1809  else if( p.nMatch("ARRA") )
1810  {
1811  /* save line array -
1812  * output energies and luminosities of predicted lines */
1813  strcpy( save.chSave[save.nsave], "LINA" );
1814  sprintf( save.chHeader[save.nsave],
1815  "#enr\tID\tI(intrinsic)\tI(emergent)\ttype\n" );
1816  }
1817 
1818  else if( p.nMatch("LABE") )
1819  {
1820  /* save line labels */
1821  strcpy( save.chSave[save.nsave], "LINL" );
1822  sprintf( save.chHeader[save.nsave],
1823  "#index\tlabel\twavelength\tcomment\n" );
1824  /* this controls whether we will print lots of redundant
1825  * info labels for transferred lines - if keyword LONG appears
1826  * then do so, if does not appear then do not - this is default */
1827  if( p.nMatch("LONG") )
1828  save.punarg[save.nsave][0] = 1;
1829  else
1830  save.punarg[save.nsave][0] = 0;
1831  }
1832 
1833  else if( p.nMatch("OPTI") )
1834  {
1835  /* save line optical depths, done in SaveLineStuff */
1836  strcpy( save.chSave[save.nsave], "LINO" );
1837 
1838  /* the default will be to make wavelengths line in the printout, called labels,
1839  * if units appears then other units will be used instead */
1840  save.chConPunEnr[save.nsave] = "labl";
1841 
1842  /* check for keyword UNITS on line, then scan wavelength or energy units if present,
1843  * units are copied into save.chConPunEnr */
1844  if( p.nMatch("UNIT") )
1845  ChkUnits(p);
1846 
1847  sprintf( save.chHeader[save.nsave],
1848  "#species\tenergy/%s\topt depth\tdamp\n",
1850 
1851  /* this is optional limit to smallest optical depths */
1852  save.punarg[save.nsave][0] = (realnum)pow(10.,p.FFmtRead());
1853  /* this is default of 0.1 napier */
1854  if( p.lgEOL() )
1855  {
1856  save.punarg[save.nsave][0] = 0.1f;
1857  }
1858  }
1859 
1860  else if( p.nMatch("POPU") )
1861  {
1862  /* save line populations command - first give index and inforamtion
1863  * for all lines, then populations for lines as a function of
1864  * depth, using this index */
1865  strcpy( save.chSave[save.nsave], "LINP" );
1866  sprintf( save.chHeader[save.nsave],
1867  "#population information\n" );
1868  /* this is optional limit to smallest population to save - always
1869  * interpreted as a log */
1870  save.punarg[save.nsave][0] = (realnum)pow(10.,p.FFmtRead());
1871 
1872  /* this is default - all positive populations */
1873  if( p.lgEOL() )
1874  save.punarg[save.nsave][0] = 0.f;
1875 
1876  if( p.nMatch(" OFF") )
1877  {
1878  /* no lower limit - print all lines */
1879  save.punarg[save.nsave][0] = -1.f;
1880  }
1881  }
1882 
1883  else if( p.nMatch("INTE") )
1884  {
1885  /* this will be full set of line intensities */
1886  strcpy( save.chSave[save.nsave], "LINI" );
1887  sprintf( save.chHeader[save.nsave],
1888  "#Emission line intrinsic intensities per unit inner area\n" );
1889  if( p.nMatch("COLU") )
1890  /* column is key to save single column */
1891  strcpy( save.chPunRltType, "column" );
1892  else
1893  /* array is key to save large array */
1894  strcpy( save.chPunRltType, "array " );
1895 
1896  save.punarg[save.nsave][0] = 0.;
1897  // ALL option - all lines, even zero intensities
1898  if( p.nMatch( " ALL" ) )
1899  save.punarg[save.nsave][0] = -1.;
1900 
1901  // check whether intrinsic or emergent line emissivity
1902  save.lgEmergent[save.nsave] = false;
1903  if( p.nMatch("EMER") )
1904  save.lgEmergent[save.nsave] = true;
1905 
1906  if( p.nMatch("EVER") )
1907  {
1908  save.LinEvery = (long int)p.FFmtRead();
1909  save.lgLinEvery = true;
1910  if( p.lgEOL() )
1911  {
1912  fprintf( ioQQQ,
1913  "There must be a second number, the number of zones to print.\nSorry.\n" );
1915  }
1916  }
1917  else
1918  {
1919  save.LinEvery = geometry.nend[0];
1920  save.lgLinEvery = false;
1921  }
1922  }
1923  else
1924  {
1925  fprintf( ioQQQ,
1926  "This option for SAVE LINE is something that I do not understand. Sorry.\n" );
1928  }
1929  }
1930 
1931  else if( p.nMatch(" MAP") )
1932  {
1933  strcpy( save.chSave[save.nsave], "MAP " );
1934  sprintf( save.chHeader[save.nsave],
1935  "#te, heating, cooling.\n" );
1936  /* do cooling space map for specified zones
1937  * if no number, or <0, do map and save out without doing first zone
1938  * does map by calling punt(" map")
1939  */
1940  hcmap.MapZone = (long)p.FFmtRead();
1941  if( p.lgEOL() )
1942  {
1943  hcmap.MapZone = 1;
1944  }
1945 
1946  if( p.nMatch("RANG") )
1947  {
1948  bool lgLogOn;
1949  hcmap.RangeMap[0] = (realnum)p.FFmtRead();
1950  if( hcmap.RangeMap[0] <= 10. && !p.nMatch("LINE") )
1951  {
1952  hcmap.RangeMap[0] = (realnum)pow((realnum)10.f,hcmap.RangeMap[0]);
1953  lgLogOn = true;
1954  }
1955  else
1956  {
1957  lgLogOn = false;
1958  }
1959 
1960  hcmap.RangeMap[1] = (realnum)p.FFmtRead();
1961  if( lgLogOn )
1962  hcmap.RangeMap[1] = (realnum)pow((realnum)10.f,hcmap.RangeMap[1]);
1963 
1964  if( p.lgEOL() )
1965  {
1966  fprintf( ioQQQ, "There must be a zone number, followed by two temperatures, on this line. Sorry.\n" );
1968  }
1969  }
1970  }
1971 
1972  else if( p.nMatch("MOLE") )
1973  {
1974  /* molecules, especially for PDR calculations */
1975  strcpy( save.chSave[save.nsave], "MOLE" );
1976  }
1977 
1978  else if( p.nMatch("MONI") )
1979  {
1980  /* save monitors */
1981  strcpy( save.chSave[save.nsave], "MONI" );
1982  }
1983 
1984  else if( p.nMatch("OPTICAL") && p.nMatch("DEPTH") )
1985  {
1986  /* check for keyword UNITS on line, then scan wavelength or energy units if present,
1987  * units are copied into save.chConPunEnr */
1988  ChkUnits(p);
1989 
1990  /* "every" option to save this on every zone -
1991  * not present then only last zone is saved */
1992  if( p.nMatch("EVER" ) )
1993  {
1994  /* save every zone */
1995  save.lgSaveEveryZone[save.nsave] = true;
1997  }
1998  else
1999  {
2000  /* only save last zone */
2001  save.lgSaveEveryZone[save.nsave] = false;
2003  }
2004 
2005  if( p.nMatch("FINE") )
2006  {
2007  /* save fine continuum optical depths */
2008  rfield.lgSaveOpacityFine = true;
2009  strcpy( save.chSave[save.nsave], "OPTf" );
2010  sprintf( save.chHeader[save.nsave], "#energy/%s\tTau tot\topacity\n",
2012  /* range option - important since so much data */
2013  if( p.nMatch("RANGE") )
2014  {
2015  /* get lower and upper range, eventually must be in Ryd */
2016  double Energy1 = p.FFmtRead();
2017  double Energy2 = p.FFmtRead();
2018  if( p.lgEOL() )
2019  {
2020  fprintf(ioQQQ,"There must be two numbers, the lower and upper energy range in Ryd.\nSorry.\n");
2022  }
2023  if( p.nMatch("UNIT" ) )
2024  {
2025  // apply units to range option
2026  const char *energyUnits = p.StandardEnergyUnit();
2027  Energy unitChange;
2028  unitChange.set(Energy1, energyUnits );
2029  Energy1 = unitChange.Ryd();
2030  unitChange.set(Energy2, energyUnits );
2031  Energy2 = unitChange.Ryd();
2032  }
2033  /* get lower and upper rang in Ryd */
2034  save.punarg[save.nsave][0] = (realnum)MIN2( Energy1 , Energy2 );
2035  save.punarg[save.nsave][1] = (realnum)MAX2( Energy1 , Energy2 );
2036  //fprintf(ioQQQ , "DEBUG units change fine %.3e %.3e\n" , save.punarg[save.nsave][0] ,
2037  // save.punarg[save.nsave][1] );
2038  //cdEXIT(EXIT_FAILURE);
2039  }
2040  else
2041  {
2042  /* these mean full energy range */
2043  save.punarg[save.nsave][0] = 0.;
2044  save.punarg[save.nsave][1] = 0.;
2045  }
2046  /* optional last parameter - how many points to bring together */
2047  save.punarg[save.nsave][2] = (realnum)p.FFmtRead();
2048  /* default is to bring together ten */
2049  if( p.lgEOL() )
2050  save.punarg[save.nsave][2] = 10;
2051  if( save.punarg[save.nsave][2] < 1 )
2052  {
2053  fprintf(ioQQQ,"The number of fine opacities to skip must be > 0 \nSorry.\n");
2055  }
2056  }
2057  else
2058  {
2059  /* save coarse continuum optical depths */
2060  strcpy( save.chSave[save.nsave], "OPTc" );
2061  sprintf( save.chHeader[save.nsave],
2062  "#energy/%s\ttotal\tabsorp\tscat\n",
2064  }
2065 
2066  }
2067  else if( p.nMatch(" OTS") )
2068  {
2069  strcpy( save.chSave[save.nsave], " OTS" );
2070  sprintf( save.chHeader[save.nsave],
2071  "#otscon, lin, conOpac LinOpc\n" );
2072  }
2073 
2074  else if( p.nMatch("OVER") && p.nMatch(" OVE") )
2075  {
2076  /* save overview of model results */
2077  strcpy( save.chSave[save.nsave], "OVER" );
2078  sprintf( save.chHeader[save.nsave],
2079  "#depth\tTe\tHtot\thden\teden\t2H_2/H\tHI\tHII\tHeI\tHeII\tHeIII\tCO/C\tC1\tC2\tC3\tC4\tO1\tO2\tO3\tO4\tO5\tO6\tH2O/O\tAV(point)\tAV(extend)\n" );
2080  }
2081 
2082  else if( p.nMatch(" PDR") )
2083  {
2084  strcpy( save.chSave[save.nsave], " PDR" );
2085  sprintf( save.chHeader[save.nsave],
2086  "#depth\tH colden\tTe\tHI/HDEN\tH2/HDEN\tH2*/HDEN\tCI/C\tCO/C\tH2O/O\tG0\tAV(point)\tAV(extend)\tTauV(point)\n" );
2087  }
2088 
2089  else if( p.nMatch("PERF") )
2090  {
2091  /* output performance characteristics per zone */
2092  strcpy( save.chSave[save.nsave], "PERF" );
2093  sprintf( save.chHeader[save.nsave],
2094  "#zone\tdTime\tElapsed t\tnPres2Ioniz\n" );
2095  }
2096 
2097  else if( p.nMatch("PHYS") )
2098  {
2099  /* save physical conditions */
2100  strcpy( save.chSave[save.nsave], "PHYS" );
2101  sprintf( save.chHeader[save.nsave],
2102  "#PhyC depth\tTe\tn(H)\tn(e)\tHtot\taccel\tfillfac\n" );
2103  }
2104 
2105  else if( p.nMatch("POIN") )
2106  {
2107  /* save out the pointers */
2108  save.lgPunPoint = true;
2109  /* this does not count as a save option (really) */
2110  strcpy( save.chSave[save.nsave], "" );
2111  save.lgRealSave[save.nsave] = false;
2112  }
2113 
2114  else if( p.nMatch("PRES") )
2115  {
2116  /* the save pressure command */
2117  strcpy( save.chSave[save.nsave], "PRES" );
2118  sprintf( save.chHeader[save.nsave],
2119  "#P depth\tPerror%%\tPcurrent\tPIn+Pinteg\tPgas(r0)\tPgas\tPram"
2120  "\tPrad(line)\tPinteg\tV(wind km/s)\tcad(wind km/s)\tP(mag)\tV(turb km/s)"
2121  "\tP(turb)\tPgr_Int\tint thin elec\tconv?\n" );
2122  }
2123 
2124  else if( p.nMatch("RADI") )
2125  {
2126  /* the save radius command */
2127  sprintf( save.chHeader[save.nsave], "#NZONE\tradius\tdepth\tdr\n" );
2128  /* option to only save the outer radius */
2129  if( p.nMatch( "OUTE" ) )
2130  {
2131  /* only outer radius */
2132  strcpy( save.chSave[save.nsave], "RADO" );
2133  }
2134  else
2135  {
2136  /* all radii */
2137  strcpy( save.chSave[save.nsave], "RADI" );
2138  }
2139  }
2140 
2141  else if( p.nMatch("RECO") )
2142  {
2143  if( p.nMatch("COEF") )
2144  {
2145  /* recombination coefficients for everything */
2146 
2147  /* this is logical flag used in routine ion_recom to create the save output */
2148  save.lgioRecom = true;
2149  /* this does not count as a save option (really) */
2150  strcpy( save.chSave[save.nsave], "" );
2151  save.lgRealSave[save.nsave] = false;
2152  }
2153 
2154  else if( p.nMatch("EFFI") )
2155  {
2156  /* save recombination efficiency */
2157  strcpy( save.chSave[save.nsave], "RECE" );
2158  sprintf( save.chHeader[save.nsave],
2159  "#Recom effic H, Heo, He+\n" );
2160  }
2161 
2162  else
2163  {
2164  fprintf( ioQQQ, "No option recognized on this save recombination command\n" );
2165  fprintf( ioQQQ, "Valid options are COEFFICIENTS, AGN, and EFFICIENCY\nSorry.\n" );
2167  }
2168  }
2169 
2170  /* save results command, either as single column or wide array */
2171  else if( p.nMatch("RESU") )
2172  {
2173  strcpy( save.chSave[save.nsave], "RESU" );
2174  if( p.nMatch("COLU") )
2175  {
2176  /* column is key to save single column */
2177  strcpy( save.chPunRltType, "column" );
2178  }
2179  else
2180  {
2181  /* array is key to save large array */
2182  strcpy( save.chPunRltType, "array " );
2183  }
2184 
2185  /* do not change following, is used as flag in getlines */
2186  sprintf( save.chHeader[save.nsave],
2187  "#results of calculation\n" );
2188  }
2189 
2190  else if( p.nMatch("SECO") )
2191  {
2192  /* save secondary ionization rate */
2193  strcpy( save.chSave[save.nsave], "SECO" );
2194  sprintf( save.chHeader[save.nsave],
2195  "#depth\tIon(H^0)\tDiss(H_2)\tExcit(Lya)\n" );
2196  }
2197 
2198  else if( p.nMatch("SOUR") )
2199  {
2200 
2201  /* check for keyword UNITS on line, then scan wavelength or energy units if present,
2202  * units are copied into save.chConPunEnr */
2203  ChkUnits(p);
2204 
2205  if( p.nMatch("DEPT") )
2206  {
2207  /* print continuum source function as function of depth */
2208  strcpy( save.chSave[save.nsave], "SOUD" );
2209  sprintf( save.chHeader[save.nsave],
2210  "#continuum source function vs depth\n" );
2211  }
2212  else if( p.nMatch("SPEC") )
2213  {
2214  /* print spectrum continuum source function at 1 depth */
2215  strcpy( save.chSave[save.nsave], "SOUS" );
2216  sprintf( save.chHeader[save.nsave],
2217  "#continuum source function nu/%s\tConEmitLocal/widflx"
2218  "\tabs opac\tConSourceFcnLocal\tConSourceFcnLocal/plankf\tConSourceFcnLocal/flux\n",
2220  }
2221  else
2222  {
2223  fprintf( ioQQQ, "A second keyword must appear on this line.\n" );
2224  fprintf( ioQQQ, "They are DEPTH and SPECTRUM.\n" );
2225  fprintf( ioQQQ, "Sorry.\n" );
2227  }
2228  }
2229 
2230 
2231  /* save spectrum the new form of the save continuum, will eventually replace the standard
2232  * save continuum command */
2233  else if( p.nMatch("SPECTRUM") && !p.nMatch("XSPE") )
2234  {
2235  /* this flag is checked in PrtComment to generate a caution
2236  * if continuum is saved but iterations not performed */
2237  save.lgPunContinuum = true;
2238 
2239  /* set flag for spectrum */
2240  strcpy( save.chSave[save.nsave], "CONN" );
2241 
2242  /* check for keyword UNITS on line, then scan wavelength or energy units if present,
2243  * units are copied into save.chConPunEnr */
2244  ChkUnits(p);
2245 
2246  sprintf( save.chHeader[save.nsave],
2247  "#Cont Enr/%s\tincid nFn\ttrans\tdiff\tlines \n",
2249  }
2250 
2251  else if( p.nMatch("SPECIAL") )
2252  {
2253  /* save special, will call routine SaveSpecial */
2254  strcpy( save.chSave[save.nsave], "SPEC" );
2255  sprintf( save.chHeader[save.nsave], "#Special.\n" );
2256  }
2257 
2258  else if( p.nMatch("SPECIES") )
2259  {
2260  strcpy( save.chSave[save.nsave], "SPCS" );
2261 
2262  // option to save information about a particular species,
2263  // the "second filename" may really be the species label. Rename here for clarity
2264  strcpy( chLabel, chSecondFilename );
2265  if( !lgSecondFilename )
2266  {
2267  strcpy( save.chSaveSpecies[save.nsave], "" );
2268  }
2269  else if( strlen(chLabel) >= CHARS_SPECIES )
2270  {
2271  fprintf( ioQQQ,"Species string is limited to %li characters.\nSorry.\n", (long)CHARS_SPECIES );
2273  }
2274  else
2275  strncpy( save.chSaveSpecies[save.nsave], chLabel, CHARS_SPECIES );
2276 
2277  if (p.nMatch( "COLUMN" ) )
2278  {
2279  /* column densities*/
2280  strcpy( save.chSaveArgs[save.nsave], "COLU" );
2281  }
2282  else if (p.nMatch( "ENERG" ) )
2283  {
2284  /* energy levels, default Rydbergs but option to change units */
2285  ChkUnits(p);
2286  strcpy( save.chSaveArgs[save.nsave], "ENER" );
2287  }
2288  else if( p.nMatch("LABELS") )
2289  {
2290  strcpy( save.chSaveArgs[save.nsave], "LABE" );
2291  }
2292  else if (p.nMatch( "LEVELS" ) )
2293  {
2294  /* the number of levels in this zone */
2295  strcpy( save.chSaveArgs[save.nsave], "LEVL" );
2296  }
2297  else if (p.nMatch( "POPUL" ) )
2298  {
2299  /* save species population fraction for 1 level*/
2300  strcpy( save.chSaveArgs[save.nsave], "POPU" );
2301  }
2302  else
2303  {
2304  fprintf( ioQQQ, "ParseSave cannot find a recognized keyword on this SAVE SPECIES command line.\n" );
2305  fprintf( ioQQQ, "I know about the keywords COLUMN DENSITIES, LABELS, LEVELS, and POPULATIONS.\nSorry.\n" );
2307  }
2308  }
2309 
2310  else if( p.nMatch("TEMP") )
2311  {
2312  /* save temperature command */
2313  strcpy( save.chSave[save.nsave], "TEMP" );
2314  sprintf( save.chHeader[save.nsave],
2315  "#depth\tTe\tcC/dT\tdt/dr\td^2T/dr^2\n" );
2316  }
2317 
2318  else if( p.nMatch("TIME") && p.nMatch("DEPE") )
2319  {
2320  /* information about time dependent solutions */
2321  strcpy( save.chSave[save.nsave], "TIMD" );
2322  /* do not want to separate iterations with special character */
2324  /* write header */
2325  sprintf( save.chHeader[save.nsave] ,
2326  "#elapsed time\ttime step \tscale cont\tn(H)\t<T>\t<H+/H rad>\t<H0/H rad>\t<H2/H rad>\t<He+/He rad>\t<CO/H>\t<redshift>\t<ne/nH>\n" );
2327  }
2328 
2329  else if( p.nMatch("TPRE") )
2330  {
2331  /* debug output from the temperature predictor in zonestart,
2332  * set with save tpred command */
2333  strcpy( save.chSave[save.nsave], "TPRE" );
2334  sprintf( save.chHeader[save.nsave],
2335  "#zone old temp, guess Tnew, new temp delta \n" );
2336  }
2337 
2338  else if( p.nMatch("WIND") )
2339  {
2340  strcpy( save.chSave[save.nsave], "WIND" );
2341  sprintf( save.chHeader[save.nsave],
2342  "#radius\tdepth\tvel [cm/s]\tTot accel [cm s-2]\tLin accel [cm s-2]"
2343  "\tCon accel [cm s-2]\tforce multiplier\ta_gravity\n" );
2344  if( p.nMatch( "TERM" ) )
2345  {
2346  /* only save for last zone, the terminal velocity, for grids */
2347  save.punarg[save.nsave][0] = 0.;
2348  }
2349  else
2350  {
2351  /* one means save every zone */
2352  save.punarg[save.nsave][0] = 1.;
2353  }
2354  }
2355 
2356  else if( p.nMatch("XSPE") )
2357  {
2358  /* say that this is a FITS file output */
2359  save.lgFITS[save.nsave] = true;
2360 
2361  /* the save xspec commands */
2362  save.lgPunLstIter[save.nsave] = true;
2363 
2364  /* remember that a save xspec command was entered */
2365  grid.lgSaveXspec = true;
2366 
2367  /* range option - important since so much data */
2368  if( p.nMatch("RANGE") )
2369  {
2370  /* get lower and upper range, must be in keV */
2371  save.punarg[save.nsave][0] = (realnum)p.FFmtRead();
2372  save.punarg[save.nsave][1] = (realnum)p.FFmtRead();
2373  if( p.lgEOL() )
2374  {
2375  fprintf(ioQQQ,"There must be two numbers, the lower and upper energy range in keV.\nSorry.\n");
2377  }
2378  if( save.punarg[save.nsave][0] >=save.punarg[save.nsave][1] )
2379  {
2380  fprintf(ioQQQ,"The two energies for the range must be in increasing order.\nSorry.\n");
2382  }
2383 
2386  }
2387  else
2388  {
2389  /* these mean full energy range */
2390  save.punarg[save.nsave][0] = 0;
2391  save.punarg[save.nsave][1] = 0;
2392  }
2393 
2394  if( p.nMatch("ATAB") )
2395  {
2396  /* save xspec atable command */
2397 
2398  if( p.nMatch("TOTA") )
2399  {
2400  /* total spectrum */
2401  strcpy( save.chSave[save.nsave], "XTOT" );
2402  grid.lgOutputTypeOn[0] = true;
2403  save.FITStype[save.nsave] = 0;
2404  }
2405  else if( p.nMatch("INCI") )
2406  {
2407  if( p.nMatch("ATTE") )
2408  {
2409  /* attenuated incident continuum */
2410  strcpy( save.chSave[save.nsave], "XATT" );
2411  grid.lgOutputTypeOn[2] = true;
2412  save.FITStype[save.nsave] = 2;
2413  }
2414  else if( p.nMatch("REFL") )
2415  {
2416  /* reflected incident continuum */
2417  strcpy( save.chSave[save.nsave], "XRFI" );
2418  grid.lgOutputTypeOn[3] = true;
2419  save.FITStype[save.nsave] = 3;
2420  }
2421  else
2422  {
2423  /* incident continuum */
2424  strcpy( save.chSave[save.nsave], "XINC" );
2425  grid.lgOutputTypeOn[1] = true;
2426  save.FITStype[save.nsave] = 1;
2427  }
2428  }
2429  else if( p.nMatch("DIFF") )
2430  {
2431  if( p.nMatch("REFL") )
2432  {
2433  /* reflected diffuse continuous emission */
2434  strcpy( save.chSave[save.nsave], "XDFR" );
2435  grid.lgOutputTypeOn[5] = true;
2436  save.FITStype[save.nsave] = 5;
2437  }
2438  else
2439  {
2440  /* diffuse continuous emission outward */
2441  strcpy( save.chSave[save.nsave], "XDFO" );
2442  grid.lgOutputTypeOn[4] = true;
2443  save.FITStype[save.nsave] = 4;
2444  }
2445  }
2446  else if( p.nMatch("LINE") )
2447  {
2448  if( p.nMatch("REFL") )
2449  {
2450  /* reflected lines */
2451  strcpy( save.chSave[save.nsave], "XLNR" );
2452  grid.lgOutputTypeOn[7] = true;
2453  save.FITStype[save.nsave] = 7;
2454  }
2455  else
2456  {
2457  /* outward lines */
2458  strcpy( save.chSave[save.nsave], "XLNO" );
2459  grid.lgOutputTypeOn[6] = true;
2460  save.FITStype[save.nsave] = 6;
2461  }
2462  }
2463  else if( p.nMatch("SPEC") )
2464  {
2465  if( p.nMatch("REFL") )
2466  {
2467  /* reflected spectrum */
2468  strcpy( save.chSave[save.nsave], "XREF" );
2469  grid.lgOutputTypeOn[9] = true;
2470  save.FITStype[save.nsave] = 9;
2471  }
2472  else
2473  {
2474  /* transmitted spectrum */
2475  strcpy( save.chSave[save.nsave], "XTRN" );
2476  grid.lgOutputTypeOn[8] = true;
2477  save.FITStype[save.nsave] = 8;
2478  }
2479  }
2480  else
2481  {
2482  /* transmitted spectrum */
2483  strcpy( save.chSave[save.nsave], "XTRN" );
2484  grid.lgOutputTypeOn[8] = true;
2485  save.FITStype[save.nsave] = 8;
2486  }
2487  }
2488  else if( p.nMatch("MTAB") )
2489  {
2490  /* save xspec mtable */
2491  strcpy( save.chSave[save.nsave], "XSPM" );
2492  grid.lgOutputTypeOn[10] = true;
2493  save.FITStype[save.nsave] = 10;
2494  }
2495  else
2496  {
2497  fprintf( ioQQQ, "Support only for xspec atable and xspec mtable.\n" );
2498  cdEXIT( EXIT_FAILURE );
2499  }
2500  }
2501 
2502  /* save column density has to come last so do not trigger specific column
2503  * densities, H2, FeII, etc.
2504  * Need both keywords since column is also the keyword for one line per line */
2505  else if( p.nMatch("COLU") && p.nMatch("DENS") )
2506  {
2507  if( p.nMatch("SOME" ))
2508  {
2509  /* flag saying save some column densities */
2510  strcpy( save.chSave[save.nsave], "COLS" );
2512  }
2513  else
2514  {
2515  /* save column densities table */
2516  strcpy( save.chSave[save.nsave], "COLU" );
2517  }
2518  }
2519  else
2520  {
2521  fprintf( ioQQQ,
2522  "ParseSave cannot find a recognized keyword on this SAVE command line.\nSorry.\n" );
2524  }
2525 
2526  /* only open if file has not already been opened during a previous call */
2527  if( save.ipPnunit[save.nsave] == NULL )
2528  {
2529  string file_name;
2530  file_name += chFilename;
2531  string mode = "w";
2532  if( save.lgFITS[save.nsave] )
2533  mode += "b";
2534 
2535  /* open the file with the name and mode generated above */
2536  save.ipPnunit[save.nsave] = open_data( file_name.c_str(), mode.c_str(), AS_LOCAL_ONLY );
2537 
2538  /* option to set no buffering for this file. The setbuf command may
2539  * ONLY be issued right after the open of the file. Giving it after
2540  * i/o has been done may result in loss of the contents of the buffer, PvH */
2541  if( p.nMatch("NO BUFFER") )
2542  setbuf( save.ipPnunit[save.nsave] , NULL );
2543  }
2544 
2545  /***************************************************************/
2546  /* */
2547  /* The following are special save options and must be done */
2548  /* after the parsing and file opening above. */
2549  /* */
2550  /* NB: these are ALSO parsed above. Here we DO something. */
2551  /* */
2552  /***************************************************************/
2553 
2554  if( p.nMatch("CONV") && p.nMatch("REAS") )
2555  {
2556  /* save reason model declared not converged
2557  * not a true save command, since done elsewhere */
2560  save.lgPunConv = true;
2561  fprintf( save.ipPunConv,
2562  "# reason for continued iterations\n" );
2563  strcpy( save.chSave[save.nsave], "" );
2564  save.lgRealSave[save.nsave] = false;
2565  }
2566 
2567  else if( p.nMatch("CONV") && p.nMatch("BASE") )
2568  {
2569  /* save some quantities we are converging */
2570  save.lgTraceConvergeBase = true;
2571  /* the second save occurrence - file has been opened,
2572  * copy handle, also pass on special no hash option */
2573  if( p.nMatch("NO HA") )
2574  save.lgTraceConvergeBaseHash = false;
2576  /* set save last flag to whatever it was above */
2578  static bool lgPrtHeader = true;
2579  if( lgPrtHeader )
2580  fprintf( save.ipTraceConvergeBase,
2581  "#zone\theat\tcool\teden\n" );
2582  lgPrtHeader = false;
2583  }
2584 
2585  else if( p.nMatch(" DR ") )
2586  {
2587  static bool lgPrtHeader = true;
2588  /* the second save dr occurrence - file has been opened,
2589  * copy handle to ipDRout, also pass on special no hash option */
2590  if( p.nMatch("NO HA") )
2591  save.lgDRHash = false;
2593  /* set save last flag to whatever it was above */
2596  if( lgPrtHeader )
2597  fprintf( save.ipDRout,
2598  "#zone\tdepth\tdr\tdr 2 go\treason \n" );
2599  lgPrtHeader = false;
2600  strcpy( save.chSave[save.nsave], "" );
2601  save.lgRealSave[save.nsave] = false;
2602  }
2603 
2604  else if( p.nMatch("QHEA") )
2605  {
2609  fprintf( gv.QHSaveFile,
2610  "#Probability distributions from quantum heating routine.\n" );
2611  save.lgRealSave[save.nsave] = false;
2612  }
2613 
2614  else if( p.nMatch("POIN") )
2615  {
2616  /* save out the pointers */
2619  save.lgPunPoint = true;
2620  fprintf( save.ipPoint,
2621  "#pointers. \n" );
2622  strcpy( save.chSave[save.nsave], "" );
2623  save.lgRealSave[save.nsave] = false;
2624  }
2625 
2626  else if( p.nMatch("RECO") && p.nMatch("COEF") )
2627  {
2628  /* recombination coefficients for everything
2629  * save.lgioRecom set to false in routine zero, non-zero value
2630  * is flag to save recombination coefficients. the output is actually
2631  * produced by a series of routines, as they generate the recombination
2632  * coefficients. these include
2633  * diel supres, helium, hydrorecom, iibod, and makerecomb*/
2636  /* this is logical flag used in routine ion_recom to create the save output */
2637  save.lgioRecom = true;
2638  fprintf( save.ioRecom,
2639  "#recombination coefficients cm3 s-1 for current density and temperature\n" );
2640  strcpy( save.chSave[save.nsave], "" );
2641  save.lgRealSave[save.nsave] = false;
2642  }
2643 
2644  else if( p.nMatch("GRID") )
2645  {
2646  /* this enables saving GRID output outside the main SaveDo() loop */
2649  }
2650 
2651  else if( p.nMatch(" MAP") )
2652  {
2653  /* say output goes to special save */
2655  }
2656 
2657  /* check that string written into save.chHeader[save.nsave] can actually fit there
2658  * we may have overrun this buffer, an internal error */
2659  /* check that there are less than nChar characters in the line */
2660  char *chEOL = strchr_s(save.chHeader[save.nsave] , '\0' );
2661 
2662  /* return null if input string longer than nChar, the longest we can read.
2663  * Print and return null but chLine still has as much of the line as
2664  * could be placed in cdLine */
2665  if( (chEOL==NULL) || (chEOL - save.chHeader[save.nsave])>=MAX_HEADER_SIZE-1 )
2666  {
2667  fprintf( ioQQQ, "DISASTER save.chHeader[%li] has been overwritten "
2668  "with a line too long to be read.\n", save.nsave );
2670  }
2671 
2672  /* if lgPunHeader true and cdHeader has been set to a string then print header
2673  * logic to prevent more than one header in grid calculation */
2675  {
2676  fprintf( save.ipPnunit[save.nsave], "%s", save.chHeader[save.nsave] );
2677  save.lgPunHeader[save.nsave] = false;
2678  }
2679 
2680  /* increment total number of save commands, */
2681  ++save.nsave;
2682  return;
2683 }
2684 
2685 /*SaveFilesInit initialize save file pointers, called from InitCoreload
2686  * called one time per core load
2687  * NB KEEP THIS ROUTINE SYNCHED UP WITH THE NEXT ONE, CloseSaveFiles */
2689 {
2690  long int i;
2691  static bool lgFIRST = true;
2692 
2693  DEBUG_ENTRY( "SaveFilesInit()" );
2694 
2695  ASSERT( lgFIRST );
2696  lgFIRST = false;
2697 
2698  /* set lgNoClobber to not overwrite files, reset with clobber on save line
2699  * if we are running a grid (grid command entered in cdRead) grid.lgGrid
2700  * true, is false if single sim. For grid we want to not clobber files
2701  * by default, do clobber for optimizer since this was behavior before */
2702  bool lgNoClobberDefault = false;
2703  if( grid.lgGrid )
2704  {
2705  /* cdRead encountered grid command - do not want to clobber files */
2706  lgNoClobberDefault = true;
2707  }
2708 
2709  for( i=0; i < LIMPUN; i++ )
2710  {
2711  save.lgNoClobber[i] = lgNoClobberDefault;
2712  }
2713  save.lgPunConv_noclobber = lgNoClobberDefault;
2714  save.lgDROn_noclobber = lgNoClobberDefault;
2715  save.lgTraceConvergeBase_noclobber = lgNoClobberDefault;
2716  save.lgPunPoint_noclobber = lgNoClobberDefault;
2717  save.lgioRecom_noclobber = lgNoClobberDefault;
2718  save.lgQHSaveFile_noclobber = lgNoClobberDefault;
2719  save.lgSaveGrid_noclobber = lgNoClobberDefault;
2720 
2721  /* initialize chHeader strings with nonsense, compare later to see if we have any actual headers. */
2722  save.chNONSENSE = "ArNdY38dZ9us4N4e12SEcuQ";
2723 
2724  for( i=0; i < LIMPUN; i++ )
2725  {
2726  save.ipPnunit[i] = NULL;
2727 
2728  // is this a real save command? set false with the dummy
2729  // save commands like save dr
2730  save.lgRealSave[i] = true;
2731 
2732  // do we need to save header?
2733  save.lgPunHeader[i] = true;
2734  strcpy( save.chHeader[i], save.chNONSENSE );
2735  }
2736 
2737  save.lgTraceConvergeBase = false;
2738 
2739  save.ipDRout = NULL;
2740  save.lgDROn = false;
2741 
2742  save.ipTraceConvergeBase = NULL;
2743  save.lgTraceConvergeBase = false;
2744 
2745  save.ipPunConv = NULL;
2746  save.lgPunConv = false;
2747 
2748  save.ipPoint = NULL;
2749  save.lgPunPoint = false;
2750 
2751  gv.QHSaveFile = NULL;
2752 
2753  save.ioRecom = NULL;
2754  save.lgioRecom = false;
2755 
2756  grid.pnunit = NULL;
2757 
2758  ioMAP = NULL;
2759 
2760  return;
2761 }
2762 
2763 /*CloseSaveFiles close save files called from cdEXIT upon termination,
2764  * from cloudy before returning
2765  * NB - KEEP THIS ROUTINE SYNCHED UP WITH THE PREVIOUS ONE, SaveFilesInit */
2766 void CloseSaveFiles( bool lgFinal )
2767 {
2768  long int i;
2769 
2770  DEBUG_ENTRY( "CloseSaveFiles()" );
2771 
2772  /* close all save units cloudy opened with save command,
2773  * lgNoClobber is set false with CLOBBER option on save, says to
2774  * overwrite the files */
2775  for( i=0; i < save.nsave; i++ )
2776  {
2777  /* if lgFinal is true, we close everything, no matter what.
2778  * this means ignoring "no clobber" options */
2779  if( save.ipPnunit[i] != NULL && ( !save.lgNoClobber[i] || lgFinal ) )
2780  {
2781  /* Test that any FITS files are the right size! */
2782  if( save.lgFITS[i] )
2783  {
2784  /* \todo 2 This overflows for file sizes larger (in bytes) than
2785  * a long int can represent (about 2GB on most 2007 systems) */
2786  fseek(save.ipPnunit[i], 0, SEEK_END);
2787  long file_size = ftell(save.ipPnunit[i]);
2788  if( file_size%2880 )
2789  {
2790  fprintf( ioQQQ, " PROBLEM FITS file is wrong size!\n" );
2791  }
2792  }
2793 
2794  fclose( save.ipPnunit[i] );
2795  save.ipPnunit[i] = NULL;
2796  }
2797  }
2798 
2799  /* following file handles are aliased to ipPnunit which was already closed above */
2800  if( save.ipDRout != NULL && ( !save.lgDROn_noclobber || lgFinal ) )
2801  {
2802  save.ipDRout = NULL;
2803  save.lgDROn = false;
2804  }
2805 
2806  if( save.ipTraceConvergeBase != NULL && ( !save.lgTraceConvergeBase_noclobber || lgFinal ) )
2807  {
2808  save.ipTraceConvergeBase = NULL;
2809  save.lgTraceConvergeBase = false;
2810  }
2811 
2812  if( save.ipPunConv != NULL && ( !save.lgPunConv_noclobber || lgFinal ) )
2813  {
2814  save.ipPunConv = NULL;
2815  save.lgPunConv = false;
2816  }
2817  if( save.ipPoint != NULL && ( !save.lgPunPoint_noclobber || lgFinal ) )
2818  {
2819  save.ipPoint = NULL;
2820  save.lgPunPoint = false;
2821  }
2822  if( gv.QHSaveFile != NULL && ( !save.lgQHSaveFile_noclobber || lgFinal ) )
2823  {
2824  gv.QHSaveFile = NULL;
2825  }
2826  if( save.ioRecom != NULL && ( !save.lgioRecom_noclobber || lgFinal ) )
2827  {
2828  save.ioRecom = NULL;
2829  save.lgioRecom = false;
2830  }
2831  if( grid.pnunit != NULL && ( !save.lgSaveGrid_noclobber || lgFinal ) )
2832  {
2833  grid.pnunit = NULL;
2834  }
2835  ioMAP = NULL;
2836 
2837  return;
2838 }
2839 
2840 /*ChkUnits check for keyword UNITS on line, then scan wavelength or energy units if present,
2841  * units are copied into save.chConPunEnr - when doing output, the routine call
2842  * AnuUnit( energy ) will automatically return the energy in the right units,
2843  * when called to do save output */
2845 {
2846 
2847  DEBUG_ENTRY( "ChkUnits()" );
2848 
2849  /* option to set units for continuum energy in save output */
2850  if( p.nMatch("UNITS") )
2851  {
2852  // p.StandardEnergyUnit() will terminate if no unit was recognized
2854  }
2855  else
2856  {
2858  }
2859  return;
2860 }
#define MIN2
Definition: cddefines.h:765
char chSpeciesDominantRates[LIMPUN][CHARS_SPECIES]
Definition: save.h:368
bool lgPunLstIter[LIMPUN]
Definition: save.h:271
void Parse_Save_Line_RT(Parser &p)
Definition: save_line.cpp:291
bool nMatch(const char *chKey) const
Definition: parser.h:135
FILE * ioMAP
Definition: cdinit.cpp:9
realnum fe2ener[2]
Definition: atomfeii.h:233
realnum punarg[LIMPUN][3]
Definition: save.h:254
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition: cpu.cpp:616
bool lgPunContinuum
Definition: save.h:251
FILE * pnunit
Definition: grid.h:58
bool lgSaveOpacityFine
Definition: rfield.h:423
bool is_odd(int j)
Definition: cddefines.h:718
realnum HiEnergy_keV
Definition: grid.h:56
double FFmtRead(void)
Definition: parser.cpp:353
FILE * QHSaveFile
Definition: grainvar.h:571
bool lgNoClobber[LIMPUN]
Definition: save.h:200
NORETURN void TotalInsanity(void)
Definition: service.cpp:886
t_input input
Definition: input.cpp:12
string chGridPrefix
Definition: save.h:302
bool lgGrid
Definition: grid.h:40
void parse_save_colden(Parser &p, char chHeader[])
Definition: save_colden.cpp:17
bool lgDRHash
Definition: save.h:335
bool lgShortFe2
Definition: atomfeii.h:225
string chFilenamePrefix
Definition: save.h:306
FILE * ipPnunit[LIMPUN]
Definition: save.h:197
char chHeader[LIMPUN][MAX_HEADER_SIZE]
Definition: save.h:231
long int MapZone
Definition: hcmap.h:20
long int * nend
Definition: geometry.h:80
void set(double energy)
Definition: energy.h:19
const int ipOXYGEN
Definition: cddefines.h:316
#define MAX2
Definition: cddefines.h:786
bool lgSaveGrid_noclobber
Definition: save.h:209
long nMatch(const char *chKey, const char *chCard)
Definition: service.cpp:451
realnum LoEnergy_keV
Definition: grid.h:56
char chOpcTyp[LIMPUN][5]
Definition: save.h:229
realnum emm
Definition: rfield.h:49
#define NFE2LEVN
Definition: atomfeii.h:180
void CloseSaveFiles(bool lgFinal)
bool nMatchErase(const char *chKey)
Definition: parser.h:158
vector< char * > chLineListLabel[LIMPUN]
Definition: save.h:178
char chVersion[INPUT_LINE_LENGTH]
Definition: version.h:19
bool lgioRecom_noclobber
Definition: save.h:206
FILE * ioRecom
Definition: save.h:345
FILE * ioQQQ
Definition: cddefines.cpp:7
char chPunRltType[7]
Definition: save.h:319
char chTitle[INPUT_LINE_LENGTH]
Definition: input.h:32
bool lgQHPunLast
Definition: grainvar.h:569
void ParseSave(Parser &p)
Definition: parse_save.cpp:51
bool lgRealSave[LIMPUN]
Definition: save.h:213
Definition: parser.h:31
bool lgTraceConvergeBase_noclobber
Definition: save.h:208
bool lgEmergent[LIMPUN]
Definition: save.h:216
const char * chNONSENSE
Definition: save.h:233
long int nsave
Definition: save.h:222
static t_version & Inst()
Definition: cddefines.h:179
t_elementnames elementnames
Definition: elementnames.cpp:5
FILE * ipPunConv
Definition: save.h:329
realnum egamry
Definition: rfield.h:52
bool lgCumulative[LIMPUN]
Definition: save.h:219
string optname[LIMPUN]
Definition: save.h:257
int GetQuote(char *chLabel, bool lgABORT)
Definition: parser.h:209
long int LinEvery
Definition: save.h:350
bool lg_separate_iterations[LIMPUN]
Definition: save.h:242
t_geometry geometry
Definition: geometry.cpp:5
bool lgPunConv
Definition: save.h:328
long int cdGetLineList(const char chFile[], vector< char * > &chLabels, vector< realnum > &wl)
FILE * ipTraceConvergeBase
Definition: save.h:342
static const long LIMPUN
Definition: save.h:11
bool lgSaveXspec
Definition: grid.h:37
static const long MAX_HEADER_SIZE
Definition: save.h:12
#define STATIC
Definition: cddefines.h:101
const double WAVNRYD
Definition: physconst.h:173
char chSaveArgs[LIMPUN][5]
Definition: save.h:265
t_rfield rfield
Definition: rfield.cpp:8
bool lgPunPoint
Definition: save.h:325
float realnum
Definition: cddefines.h:107
const char * StandardEnergyUnit(void) const
Definition: parser.cpp:174
bool lgPunPoint_noclobber
Definition: save.h:205
const int NUM_OUTPUT_TYPES
Definition: grid.h:21
#define EXIT_FAILURE
Definition: cddefines.h:144
const int INPUT_LINE_LENGTH
Definition: cddefines.h:258
long int nSaveEveryZone[LIMPUN]
Definition: save.h:262
bool lgQHSaveFile_noclobber
Definition: save.h:207
#define cdEXIT(FAIL)
Definition: cddefines.h:438
bool lgLineListRatio[LIMPUN]
Definition: save.h:182
bool lgDRPLst
Definition: save.h:335
FILE * ipPoint
Definition: save.h:324
NORETURN void NoNumb(const char *chDesc) const
Definition: parser.cpp:233
const char * strchr_s(const char *s, int c)
Definition: cddefines.h:1443
long int GetElem(void) const
Definition: parser.cpp:209
diatomics h2("h2", 4100.,&hmi.H2_total, Yan_H2_CS)
char chElementNameShort[LIMELM][CHARS_ELEMENT_NAME_SHORT]
Definition: elementnames.h:21
t_grid grid
Definition: grid.cpp:5
char chElementSym[LIMELM][CHARS_ELEMENT_SYM]
Definition: elementnames.h:25
STATIC void ChkUnits(Parser &p)
void SaveLineListFree(long i)
Definition: save.h:157
bool lgSaveEveryZone[LIMPUN]
Definition: save.h:261
const char * chConPunEnr[LIMPUN]
Definition: save.h:284
bool lgLinEvery
Definition: save.h:351
realnum RangeMap[2]
Definition: hcmap.h:23
char chSaveSpecies[LIMPUN][CHARS_SPECIES]
Definition: save.h:267
#define ASSERT(exp)
Definition: cddefines.h:582
void sprt_wl(char *chString, realnum wl)
Definition: prt.cpp:25
int FITStype[LIMPUN]
Definition: save.h:277
bool lgFITS[LIMPUN]
Definition: save.h:274
const char * StandardEnergyUnit(const char *chCard)
Definition: energy.cpp:47
vector< realnum > wlLineList[LIMPUN]
Definition: save.h:180
bool lgHashEndIter[LIMPUN]
Definition: save.h:291
const int LIMELM
Definition: cddefines.h:262
double Ryd() const
Definition: energy.h:26
Definition: energy.h:7
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:688
const int ipHELIUM
Definition: cddefines.h:310
bool lgTraceConvergeBase
Definition: save.h:340
bool lgDROn
Definition: save.h:335
t_FeII FeII
Definition: atomfeii.cpp:5
diatomics hd("hd", 4100.,&hmi.HD_total, Yan_H2_CS)
bool lgOutputTypeOn[NUM_OUTPUT_TYPES]
Definition: grid.h:53
bool lgEOL(void) const
Definition: parser.h:98
bool lgioRecom
Definition: save.h:346
const int ipCARBON
Definition: cddefines.h:314
t_hcmap hcmap
Definition: hcmap.cpp:21
void parse_save_line(Parser &p, bool lgLog3, char *chHeader)
Definition: save_line.cpp:33
void parse_save_average(Parser &p, long int ipPun, char *chHeader)
GrainVar gv
Definition: grainvar.cpp:5
bool lgTraceConvergeBaseHash
Definition: save.h:340
bool lgPunConv_noclobber
Definition: save.h:203
t_save save
Definition: save.cpp:5
FILE * ipDRout
Definition: save.h:334
const int ipHYDROGEN
Definition: cddefines.h:309
Energy emisfreq[LIMPUN]
Definition: save.h:371
bool lgSaveToSeparateFiles[LIMPUN]
Definition: save.h:238
long nLineList[LIMPUN]
Definition: save.h:176
char chSave[LIMPUN][5]
Definition: save.h:225
void H2_ParseSave(Parser &p, char *chHeader)
Definition: mole_h2_io.cpp:111
bool lgDROn_noclobber
Definition: save.h:204
void SaveFilesInit(void)
realnum fe2thresh
Definition: atomfeii.h:236
bool lgPunHeader[LIMPUN]
Definition: save.h:246