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