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