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