iter_startend.cpp

Go to the documentation of this file.

/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
 * others.  For conditions of distribution and use see copyright notice in license.txt */
/*IterStart set and save values of many variables at start of iteration after initial temperature set*/
/*IterRestart restart iteration */
#include "cddefines.h"
#include "cddrive.h"
#include "physconst.h"
#include "iso.h"
#include "taulines.h"
#include "hydrogenic.h"
#include "struc.h"
#include "dynamics.h"
#include "prt.h"
#include "hyperfine.h"
#include "dense.h"
#include "magnetic.h"
#include "continuum.h"
#include "geometry.h"
#include "h2.h"
#include "co.h"
#include "he.h"
#include "grains.h"
#include "atomfeii.h"
#include "pressure.h"
#include "stopcalc.h"
#include "conv.h"
#include "mean.h"
#include "ca.h"
#include "thermal.h"
#include "atoms.h"
#include "wind.h"
#include "opacity.h"
#include "timesc.h"
#include "trace.h"
#include "colden.h"
#include "secondaries.h"
#include "hmi.h"
#include "radius.h"
#include "phycon.h"
#include "called.h"
#include "mole.h"
#include "rfield.h"
#include "ionbal.h"
#include "atmdat.h"
#include "lines.h"
#include "molcol.h"
#include "input.h"
#include "rt.h"
#include "iterations.h"
#include "cosmology.h"
#include "deuterium.h"
/* these are the static saved variables, set in IterStart and reset in IterRestart */

static double h2plus_heat_save, HeatH2Dish_used_save, HeatH2Dexc_used_save,
        hmihet_save, hmitot_save, H2_Solomon_dissoc_rate_used_H2g_save,deriv_HeatH2Dexc_used_save,
        H2_Solomon_dissoc_rate_used_H2s_save, H2_H2g_to_H2s_rate_used_save, H2_photodissoc_used_H2s_save,
        H2_photodissoc_used_H2g_save, 
        UV_Cont_rel2_Draine_DB96_face,
        UV_Cont_rel2_Draine_DB96_depth , UV_Cont_rel2_Habing_TH85_face ,
        **saveMoleSource , **saveMoleSink;
static realnum ***SaveMoleChTrRate;

/* arguments are atomic number, ionization stage*/
static realnum xIonFsave[LIMELM][LIMELM+1];
static realnum deutDenseSave0, deutDenseSave1;
static double HeatSave[LIMELM][LIMELM];
static realnum supsav[LIMELM][LIMELM],p2nit,d5200r;
/* save values of dr and drNext */
static double drSave , drNextSave;

/* also places to save them between iterations */
static long int IonLowSave[LIMELM], 
        IonHighSave[LIMELM];

/* these are used to reset the level populations of the h and he model atoms */
/*static realnum hnsav[LIMELM][LMHLVL+1], */
static bool lgHNSAV = false;

static double ***HOpacRatSav;

static double ortho_save , para_save;
static double ***hnsav,
  edsav;

static realnum gas_phase_save[LIMELM];
static double *den_save;

/*IterStart set and save values of many variables at start of iteration after initial temperature set*/
void IterStart(void)
{
        long int i, 
          ion, 
                ion2,
          ipISO,
          n ,
          nelem;
        double fhe1nx, 
          ratio;

        DEBUG_ENTRY( "IterStart()" );

        /* allocate two matrices if first time in this core load 
         * these variables are malloced here because they are file static - 
         * used to save values at start of first zone, and reset to this value at
         * start of new iteration */
        if( !lgHNSAV )
        {
                /* set flag so we never do this again */
                lgHNSAV = true;

                HOpacRatSav = (double ***)MALLOC(sizeof(double **)*NISO );

                hnsav = (double ***)MALLOC(sizeof(double **)*NISO );

                for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
                {

                        HOpacRatSav[ipISO] = (double **)MALLOC(sizeof(double *)*LIMELM );

                        hnsav[ipISO] = (double **)MALLOC(sizeof(double *)*LIMELM );

                        /* now do the second dimension */
                        for( nelem=ipISO; nelem<LIMELM; ++nelem )
                        {
                                HOpacRatSav[ipISO][nelem] = (double *)MALLOC(sizeof(double)*(unsigned)(iso_sp[ipISO][nelem].numLevels_max+1) );

                                hnsav[ipISO][nelem] = (double *)MALLOC(sizeof(double)*(unsigned)(iso_sp[ipISO][nelem].numLevels_max+1) );
                        }
                }
                saveMoleSource = 
                        (double**)MALLOC(sizeof(double*)*(unsigned)LIMELM );
                saveMoleSink = 
                        (double**)MALLOC(sizeof(double*)*(unsigned)LIMELM );
                SaveMoleChTrRate = 
                        (realnum***)MALLOC(sizeof(realnum**)*(unsigned)LIMELM );
                for(nelem=0; nelem<LIMELM; ++nelem )
                {
                        /* chemistry source and sink terms for ionization ladders */
                        saveMoleSource[nelem] = 
                                (double*)MALLOC(sizeof(double)*(unsigned)(nelem+2) );
                        saveMoleSink[nelem] = 
                                (double*)MALLOC(sizeof(double)*(unsigned)(nelem+2) );
                        SaveMoleChTrRate[nelem] = 
                                (realnum**)MALLOC(sizeof(realnum*)*(unsigned)(nelem+2) );
                        for( ion=0; ion<nelem+2; ++ion )
                        {
                                SaveMoleChTrRate[nelem][ion] = 
                                        (realnum*)MALLOC(sizeof(realnum)*(unsigned)(nelem+2) );
                        }
                }
                den_save = (double*)MALLOC(sizeof(double)*(mole_global.num_calc) );
        }

        /* IterStart is called at the start of EVERY iteration */
        if( trace.lgTrace )
        {
                fprintf( ioQQQ, " IterStart called.\n" );
        }

        /* check if this is the last iteration */
        if( iteration > iterations.itermx )
                iterations.lgLastIt = true;

        /* flag set true in RT_line_one_tauinc if maser ever capped */
        rt.lgMaserCapHit = false;
        /* flag for remembering if maser ever set dr in any zone */
        rt.lgMaserSetDR = false;

        /* zero out charge transfer heating and cooling fractions */
        atmdat.HCharHeatMax = 0.;
        atmdat.HCharCoolMax = 0.;

        /* these are set false if limits of Case B tables is ever exceeded */
        for(nelem=0; nelem<HS_NZ; ++nelem )
        {
                atmdat.lgHCaseBOK[0][nelem] = true;
                atmdat.lgHCaseBOK[1][nelem] = true;
        }

        /* reset the largest number of levels in the database species */
        for( long ipSpecies=0; ipSpecies<nSpecies; ++ipSpecies )
        {
                dBaseSpecies[ipSpecies].numLevels_local = dBaseSpecies[ipSpecies].numLevels_max;
                dBaseSpecies[ipSpecies].CoolTotal = 0.;
        }

        /* the reason this calculation stops */
        strncpy( StopCalc.chReasonStop, "reason not specified.", 
                sizeof(StopCalc.chReasonStop) );

        /* zero fractions of He0 destruction due to 23S */
        he.nzone = 0;
        he.frac_he0dest_23S = 0.;
        he.frac_he0dest_23S_photo = 0.;

        dense.EdenMax = 0.;
        dense.EdenMin = BIGFLOAT;

        timesc.time_H2_Dest_longest = 0.;
        timesc.time_H2_Form_longest = 0.;
        /* remains neg if not evaluated */
        timesc.time_H2_Dest_here = -1.;
        timesc.time_H2_Form_here = 0.;

        timesc.BigCOMoleForm = 0.;
        timesc.TimeH21cm = 0.;
        thermal.CoolHeatMax = 0.;
        hydro.HCollIonMax = 0.;

        atmdat.HIonFracMax = 0.;

        /* will save highest n=2p/1s ratio for hydrogen atom*/
        hydro.pop2mx = 0.;
        hydro.lgHiPop2 = false;

        hydro.nLyaHot = 0;
        hydro.TLyaMax = 0.;

        /* evaluate the gas and radiation pressure */
        /* this sets values of pressure.PresTotlCurr */
        pressure.PresInteg = 0.;
        pressure.pinzon = 0.;
        pressure.PresIntegElecThin = 0.;
        PresTotCurrent();

        dynamics.HeatMax = 0.;
        dynamics.CoolMax = 0.;
        if( dynamics.lgAdvection )
                DynaIterStart();

        /* reset these since we now have a valid solution */
        pressure.pbeta = 0.;
        pressure.RadBetaMax = 0.;
        pressure.lgPradCap = false;
        pressure.lgPradDen = false;
        /* this flag will say we hit the sonic point */
        pressure.lgSonicPoint = false;
        pressure.lgStrongDLimbo = false;

        /* define two timescales for equilibrium, Compton and thermal */
        timesc.tcmptn = 1.e0/(rfield.qtot*6.65e-25*dense.eden);
        timesc.time_therm_long = 1.5*dense.pden*BOLTZMANN*phycon.te/thermal.ctot;

        /* this will be total mass in computed structure */
        dense.xMassTotal = 0.;

        for( nelem=0; nelem < LIMELM; nelem++ )
        {

                if( dense.lgElmtOn[nelem] )
                {

                        /* now zero out the ionic fractions */
                        for( ion=0; ion < (nelem + 2); ion++ )
                        {
                                xIonFsave[nelem][ion] = dense.xIonDense[nelem][ion];
                        }

                        for( ion=0; ion < LIMELM; ion++ )
                        {
                                HeatSave[nelem][ion] = thermal.heating[nelem][ion];
                        }
                }
        }

        deutDenseSave0 = deut.xIonDense[0];
        deutDenseSave1 = deut.xIonDense[1];

        /* >>chng 02 jan 28, add ipISO loop */
        for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
        {
                /* >>chng 03 apr 11, from 0 to ipISO */
                for( nelem=ipISO; nelem < LIMELM; nelem++ )
                {
                        if( dense.lgElmtOn[nelem] )
                        {
                                for( n=0; n < iso_sp[ipISO][nelem].numLevels_max; n++ )
                                {
                                        HOpacRatSav[ipISO][nelem][n] = iso_sp[ipISO][nelem].fb[n].ConOpacRatio;
                                        hnsav[ipISO][nelem][n] = iso_sp[ipISO][nelem].st[n].Pop();
                                }
                        }
                        for( vector<two_photon>::iterator tnu = iso_sp[ipISO][nelem].TwoNu.begin(); tnu != iso_sp[ipISO][nelem].TwoNu.end(); ++tnu )
                                tnu->induc_dn_max = 0.;
                }
        }

        rfield.ipEnergyBremsThin = 0;
        rfield.EnergyBremsThin = 0.;

        p2nit = atoms.p2nit;
        d5200r = atoms.d5200r;
        atoms.nNegOI = 0;
        for( i=0; i< N_OI_LEVELS; ++i )
                atoms.popoi[i] = 0.;

        /* save molecular fractions and ionization range */
        for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
        {
                /* save molecular densities */
                gas_phase_save[nelem] = dense.gas_phase[nelem];
                IonLowSave[nelem] = dense.IonLow[nelem];
                IonHighSave[nelem] = dense.IonHigh[nelem];
        }
        /*fprintf(ioQQQ,"DEBUG IterStart save set to gas_phase hden %.4e \n",
                dense.gas_phase[0]);*/

        edsav = dense.eden;

        /* save molecular densities */
        for( i=0; i<mole_global.num_calc; i++) 
        {
                den_save[i] = mole.species[i].den;
        }
        ortho_save = h2.ortho_density;
        para_save = h2.para_density;

        for( nelem=ipHYDROGEN; nelem < LIMELM; ++nelem )
        {
                /* these have one more ion than above */
                for( ion=0; ion<nelem+2; ++ion )
                {
                        /* zero out the source and sink arrays */
                        saveMoleSource[nelem][ion] = mole.source[nelem][ion];
                        saveMoleSink[nelem][ion] = mole.sink[nelem][ion];
                        /*>>chng 06 jun 27, move from here, where leak occurs, to
                        * above where xMoleChTrRate first created */
                        /*mole.xMoleChTrRate[nelem][ion] = (realnum*)MALLOC(sizeof(realnum)*(unsigned)(nelem+2));*/
                        for( ion2=0; ion2<nelem+2; ++ion2 )
                        {
                                SaveMoleChTrRate[nelem][ion][ion2] = mole.xMoleChTrRate[nelem][ion][ion2];
                        }
                }
        }

        hmihet_save = hmi.hmihet;
        hmitot_save = hmi.hmitot;

        h2plus_heat_save = hmi.h2plus_heat;

        HeatH2Dish_used_save = hmi.HeatH2Dish_used;
        HeatH2Dexc_used_save = hmi.HeatH2Dexc_used;

        deriv_HeatH2Dexc_used_save = hmi.deriv_HeatH2Dexc_used;
        H2_Solomon_dissoc_rate_used_H2g_save = hmi.H2_Solomon_dissoc_rate_used_H2g;
        H2_Solomon_dissoc_rate_used_H2s_save = hmi.H2_Solomon_dissoc_rate_used_H2s;

        H2_H2g_to_H2s_rate_used_save = hmi.H2_H2g_to_H2s_rate_used;

        H2_photodissoc_used_H2s_save = hmi.H2_photodissoc_used_H2s;
        H2_photodissoc_used_H2g_save = hmi.H2_photodissoc_used_H2g;

        UV_Cont_rel2_Draine_DB96_face = hmi.UV_Cont_rel2_Draine_DB96_face;
        UV_Cont_rel2_Draine_DB96_depth = hmi.UV_Cont_rel2_Draine_DB96_depth;
        UV_Cont_rel2_Habing_TH85_face = hmi.UV_Cont_rel2_Habing_TH85_face;

        /* save zone thickness, and next zone thickness */
        drSave = radius.drad;
        drNextSave = radius.drNext;
        /* remember largest and smallest dr over iteration */
        radius.dr_min_last_iter = BIGFLOAT;
        radius.dr_max_last_iter = 0.;

        geometry.nprint = iterations.IterPrnt[iteration-1];

        colden.TotMassColl = 0.;
        colden.tmas = 0.;
        colden.wmas = 0.;
        colden.rjnmin = FLT_MAX;
        colden.ajmmin = FLT_MAX;

        thermal.nUnstable = 0;
        thermal.lgUnstable = false;

        rfield.extin_mag_B_point = 0.;
        rfield.extin_mag_V_point = 0.;
        rfield.extin_mag_B_extended = 0.;
        rfield.extin_mag_V_extended = 0.;

        /* plus 1 is to zero out point where unit integration occurs */
        for( i=0; i < rfield.nupper+1; i++ )
        {
                /* diffuse fields continuum */
                /* >>chng 03 nov 08, recover SummedDif */
                rfield.SummedDifSave[i] = rfield.SummedDif[i];
                rfield.ConEmitReflec[0][i] = 0.;
                rfield.ConEmitOut[0][i] = 0.;
                rfield.ConInterOut[i] = 0.;
                /* save OTS continuum for next iteration */
                rfield.otssav[i][0] = rfield.otscon[i];
                rfield.otssav[i][1] = rfield.otslin[i];
                rfield.outlin[0][i] = 0.;
                rfield.outlin_noplot[i] = 0.;
                rfield.ConOTS_local_OTS_rate[i] = 0.;
                rfield.ConOTS_local_photons[i] = 0.;
                rfield.DiffuseEscape[i] = 0.;

                /* save initial absorption, scattering opacities for next iteration */
                opac.opacity_abs_savzon1[i] = opac.opacity_abs[i];
                opac.opacity_sct_savzon1[i] = opac.opacity_sct[i];

                /* will accumulate optical depths through cloud */
                opac.TauAbsFace[i] = opac.taumin;
                opac.TauScatFace[i] = opac.taumin;
                /* >>chng 99 dec 04, having exactly 1 zone thickness for first zone caused discontinuity
                 * for heating in very high T models in func_map.in.  zone 1 and 2 were 20% different,
                 * since tau in is 1e-20, e2 is 0.9999, and so some H ots
                 * these were not here at all - changed to dr/2 */
                /* attenuation of flux by optical depths IN THIS ZONE 
                 * DirectionalCosin is 1/COS(theta), is usually 1, reset with illuminate command,
                 * option for illumination of slab at an angle */
                /* >>chng 04 oct 09, from drad to radius.drad_x_fillfac - include fill fac, PvH */
                opac.ExpZone[i] = sexp(opac.opacity_abs[i]*radius.drad_x_fillfac/2.*geometry.DirectionalCosin);

                /* e(-tau) in inward direction, up to illuminated face */
                opac.ExpmTau[i] = (realnum)opac.ExpZone[i];

                /* e2(tau) in inward direction, up to illuminated face, this is used to get the
                 * recombination escape probability */
                opac.E2TauAbsFace[i] = (realnum)e2(opac.TauAbsFace[i] );
        }

        t_fe2ovr_la::Inst().zero_opacity();

        /* this zeros out arrays to hold mean ionization fractions
         * later entered by mean
         * read out by setlim */
        mean.MeanZero();

        /* zero out the column densities */
        for( i=0; i < NCOLD; i++ )
        {
                colden.colden[i] = 0.;
        }
        colden.He123S = 0.;
        colden.H0_ov_Tspin = 0.;
        colden.OH_ov_Tspin = 0.;
        colden.coldenH2_ov_vel = 0.;

        /* upper and lower levels of H0 1s */
        colden.H0_21cm_upper =0;
        colden.H0_21cm_lower =0;

        for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
        {
                (*diatom)->ortho_colden = 0.;
                (*diatom)->para_colden = 0.;
        }

        for( i=0; i < mole_global.num_calc; i++ )
        {
                mole.species[i].column = 0.;
        }
        for( i=0; i < 5; i++ )
        {
                colden.C2Pops[i] = 0.;
                colden.C2Colden[i] = 0.;
                /* pops and column density for SiII atom */
                colden.Si2Pops[i] = 0.;
                colden.Si2Colden[i] = 0.;
        }
        for( i=0; i < 3; i++ )
        {
                /* pops and column density for CI atom */
                colden.C1Pops[i] = 0.;
                colden.C1Colden[i] = 0.;
                /* pops and column density for OI atom */
                colden.O1Pops[i] = 0.;
                colden.O1Colden[i] = 0.;
                /* pops and column density for CIII atom */
                colden.C3Pops[i] = 0.;
        }
        for( i=0; i < 4; i++ )
        {
                /* pops and column density for CIII atom */
                colden.C3Colden[i] = 0.;
        }

        /* these are some line of sight emission measures */
        colden.dlnenp = 0.;
        colden.dlnenHep = 0.;
        colden.dlnenHepp = 0.;
        colden.dlnenCp = 0.;
        colden.H0_ov_Tspin = 0.;
        colden.OH_ov_Tspin = 0.;

        // zero column densities of all states
        for( unsigned i = 0; i < mole.species.size(); ++i )
        {
                if( mole.species[i].levels != NULL )
                {
                        for( qList::iterator st = mole.species[i].levels->begin(); st != mole.species[i].levels->end(); ++st )
                        {
                                (*st).ColDen() = 0.;    
                        }
                }
        }

        /* now zero heavy element molecules */
        molcol("ZERO",ioQQQ);

        /* this will be sum of all free free heating over model */
        thermal.FreeFreeTotHeat = 0.;

        thermal.thist = 0.;
        thermal.tlowst = 1e20f;

        wind.AccelAver = 0.;
        wind.acldr = 0.;
        ionbal.ilt = 0;
        ionbal.iltln = 0;
        ionbal.ilthn = 0;
        ionbal.ihthn = 0;
        ionbal.ifail = 0;

        secondaries.SecHIonMax = 0.;
        for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
        {
                for( ion=0; ion<nelem+1; ++ion )
                {
                        supsav[nelem][ion] = secondaries.csupra[nelem][ion];
                }
        }
        secondaries.x12sav = secondaries.x12tot;
        secondaries.hetsav = secondaries.HeatEfficPrimary;
        secondaries.savefi = secondaries.SecIon2PrimaryErg;
        for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
        {
                for( ion=0; ion<nelem+1; ++ion )
                {
                        ionbal.CompRecoilHeatRateSave[nelem][ion] = ionbal.CompRecoilHeatRate[nelem][ion];
                        ionbal.CompRecoilIonRateSave[nelem][ion] = ionbal.CompRecoilIonRate[nelem][ion];
                }
        }

        /* these will keep track of the number of convergence failures that occur */
        conv.nTeFail = 0;
        conv.nTotalFailures = 0;
        conv.nPreFail = 0;
        conv.failmx = 0.;
        conv.nIonFail = 0;
        conv.nPopFail = 0;
        conv.nNeFail = 0;
        conv.nGrainFail = 0;
        conv.nChemFail = 0;
        conv.dCmHdT = 0.;

        /* abort flag */
        lgAbort = false;

        GrainStartIter();

        rfield.comtot = 0.;

        co.codfrc = 0.;
        co.codtot = 0.;

        hmi.HeatH2DexcMax = 0.;
        hmi.CoolH2DexcMax = 0.;
        hmi.h2dfrc = 0.;
        hmi.h2line_cool_frac = 0.;
        hmi.h2dtot = 0.;
        thermal.HeatLineMax = 0.;
        thermal.GBarMax = 0.;
        hyperfine.cooling_max = 0.;
        hydro.cintot = 0.;
        geometry.lgZoneTrp = false;

        hmi.h2pmax = 0.;

        /************************************************************************
         *
         * allocate space for lines arrays 
         *
         ************************************************************************/



        /* this was set in call to lines above */
        ASSERT( LineSave.nsum > 0);
        ASSERT( LineSave.nsumAllocated == LineSave.nsum );

        /* zero emission line arrays - this has to be done on every iteration */
        for( i=0; i < LineSave.nsum; i++ )
        {
                LineSv[i].SumLine[0] = 0.;
                LineSv[i].SumLine[1] = 0.;
                LineSv[i].emslin[0] = 0.;
                LineSv[i].emslin[1] = 0.;
        }

        /* now zero out some variables set by last call to LINES */
        hydro.cintot = 0.;
        thermal.totcol = 0.;
        rfield.comtot = 0.;
        thermal.FreeFreeTotHeat = 0.;
        thermal.power = 0.;
        thermal.HeatLineMax = 0.;
        thermal.GBarMax = 0.;
        hyperfine.cooling_max = 0.;

        hmi.h2pmax = 0.;

        co.codfrc = 0.;
        co.codtot = 0.;

        hmi.h2dfrc = 0.;
        hmi.h2line_cool_frac = 0.;
        hmi.h2dtot = 0.;
        timesc.sound = 0.;

        if( LineSave.lgNormSet )
        {
                /* set normalization line index to zero from negative initial value so that
                 * we pass the assert in cdLine, in order to get the norm line index */
                LineSave.ipNormWavL = 0;
                LineSave.ipNormWavL = cdLine( LineSave.chNormLab , LineSave.WavLNorm , &fhe1nx , &ratio );
                if( LineSave.ipNormWavL < 0 )
                {
                        /* did not find the line if return is negative */
                        fprintf( ioQQQ, "PROBLEM could not find the normalisation line.\n");
                        fprintf( ioQQQ, 
                        "IterStart could not find a line with a wavelength of %f and a label of %s\n", 
                          LineSave.WavLNorm, LineSave.chNormLab );
                        fprintf( ioQQQ, "Please check the emission line output to find the correct line identification.\n");
                        fprintf( ioQQQ, "Sorry.\n");
                        cdEXIT(EXIT_FAILURE);
                }
        }
        else
        {
                /* set normalization line index to zero from negative initial value so that
                 * we pass the assert in cdLine, in order to get the norm line index */
                LineSave.ipNormWavL = 0;
                LineSave.ipNormWavL = cdLine( "TOTL" , 4861.36f , &fhe1nx , &ratio );
                if( LineSave.ipNormWavL < 0 )
                {
                        /* did not find the line if return is negative */
                        fprintf( ioQQQ, "PROBLEM could not find the default normalisation line.\n");
                        fprintf( ioQQQ, 
                        "IterStart could not find a line with a wavelength of 4861 and a label of TOTL\n" );
                        fprintf( ioQQQ, "Please check the emission line output to find the correct line identification.\n");
                        fprintf( ioQQQ, "Sorry.\n");
                        cdEXIT(EXIT_FAILURE);
                }
        }

        /* set up stop line command on first iteration 
         * find index for lines and save for future iterations 
         * StopCalc.nstpl is zero (false) if no stop line commands entered */
        if( iteration == 1 && StopCalc.nstpl )
        {
                /* nstpl is number of stop line commands, 0 if none entered */
                for( long int nStopLine=0; nStopLine < StopCalc.nstpl; nStopLine++ )
                {
                        double relint, absint ;
                        /* returns array index for line in array stack if we found the line, 
                         * return negative of total number of lines as debugging aid if line not found */
                        StopCalc.ipStopLin1[nStopLine] = cdLine( StopCalc.chStopLabel1[nStopLine], 
                                /* wavelength of line in angstroms, not format printed by code */
                                StopCalc.StopLineWl1[nStopLine], &relint, &absint );

                        if( StopCalc.ipStopLin1[nStopLine]<0 )
                        {
                                fprintf( ioQQQ, 
                                  " IterStart could not find first line in STOP LINE command, number %ld with label *%s* and wl %.1f\n", 
                                  StopCalc.ipStopLin1[nStopLine] ,
                                  StopCalc.chStopLabel1[nStopLine],
                                  StopCalc.StopLineWl1[nStopLine]);
                                cdEXIT(EXIT_FAILURE);
                        }

                        StopCalc.ipStopLin2[nStopLine] = cdLine( StopCalc.chStopLabel2[nStopLine], 
                                /* wavelength of line in angstroms, not format printed by code */
                                StopCalc.StopLineWl2[nStopLine], &relint, &absint );

                        if( StopCalc.ipStopLin2[nStopLine] < 0 )
                        {
                                fprintf( ioQQQ, 
                                " IterStart could not find second line in STOP LINE command, line number %ld with label *%s* and wl %.1f\n", 
                                StopCalc.ipStopLin1[nStopLine] ,
                                StopCalc.chStopLabel1[nStopLine],
                                StopCalc.StopLineWl1[nStopLine]);
                                cdEXIT(EXIT_FAILURE);
                        }

                        if( trace.lgTrace )
                        {
                                fprintf( ioQQQ, 
                                        " stop line 1 is number %5ld wavelength is %f label is %4.4s\n", 
                                        StopCalc.ipStopLin1[nStopLine], 
                                        LineSv[StopCalc.ipStopLin1[nStopLine]].wavelength, 
                                        LineSv[StopCalc.ipStopLin1[nStopLine]].chALab );
                                fprintf( ioQQQ, 
                                        " stop line 2 is number %5ld wavelength is %f label is %4.4s\n", 
                                  StopCalc.ipStopLin2[nStopLine], 
                                  LineSv[StopCalc.ipStopLin2[nStopLine]].wavelength, 
                                  LineSv[StopCalc.ipStopLin2[nStopLine]].chALab  );
                        }
                }
        }

        /* option to only print last iteration */
        if( prt.lgPrtLastIt )
        {
                if( iteration == 1 )
                {
                        called.lgTalk = false;
                }

                /* initial condition of TALK may be off if optimization used or not master rank
                 * sec part is for print last command 
                 * lgTalkForcedOff is set true when cdTalk is called with false
                 * to turn off printing */
                if( iterations.lgLastIt && !prt.lgPrtStart && !called.lgTalkForcedOff )
                {
                        called.lgTalk = called.lgTalkSave;
                }
        }

        if( opac.lgCaseB )
        {
                if( trace.lgTrace )
                {
                        fprintf( ioQQQ, " IterStart does not change mid-zone optical depth since CASE B\n" );
                }
        }

        /* check if induced recombination can be ignored */
        hydro.FracInd = 0.;
        hydro.fbul = 0.;

        /* remember some things about effects of induced rec on H only
         * don't do ground state since SPHERE turns it off */
        for( i=ipH2p; i < (iso_sp[ipH_LIKE][ipHYDROGEN].numLevels_max - 1); i++ )
        {
                if( iso_sp[ipH_LIKE][ipHYDROGEN].trans(i+1,i).Emis().Aul() <= iso_ctrl.SmallA )
                        continue;

                ratio = iso_sp[ipH_LIKE][ipHYDROGEN].fb[i].RecomInducRate*iso_sp[ipH_LIKE][ipHYDROGEN].fb[i].PopLTE/
                        SDIV(iso_sp[ipH_LIKE][ipHYDROGEN].fb[i].RecomInducRate*iso_sp[ipH_LIKE][ipHYDROGEN].fb[i].PopLTE + 
                        iso_sp[ipH_LIKE][ipHYDROGEN].fb[i].RadRecomb[ipRecRad]*
                        iso_sp[ipH_LIKE][ipHYDROGEN].fb[i].RadRecomb[ipRecNetEsc]);
                if( ratio > hydro.FracInd )
                {
                        hydro.FracInd = (realnum)ratio;
                        hydro.ndclev = i;
                }

                ratio = iso_sp[ipH_LIKE][ipHYDROGEN].trans(i+1,i).Emis().pump()/
                        (iso_sp[ipH_LIKE][ipHYDROGEN].trans(i+1,i).Emis().pump() + 
                        iso_sp[ipH_LIKE][ipHYDROGEN].trans(i+1,i).Emis().Aul());

                if( ratio > hydro.fbul )
                {
                        hydro.fbul = (realnum)ratio;
                        hydro.nbul = i;
                }
        }

        if( trace.lgTrace )
                fprintf( ioQQQ, " IterStart returns.\n" );
        return;
}

/*IterRestart reset many variables at the start of a new iteration 
 * called by cloudy after calculation is completed, when more iterations 
 * are needed - the iteration has been incremented when this routine is
 * called so iteration == 2 after first iteration, we are starting
 * the second iteration */
void IterRestart(void)
{
        long int i, 
          ion, 
                ion2,
          ipISO ,
          n, 
          nelem;
        double SumOTS;

        DEBUG_ENTRY( "IterRestart()" );

        /* this is case where temperature floor has been set, if it was hit
         * then we did a constant temperature calculation, and must go back to 
         * a thermal solution 
         * test on thermal.lgTemperatureConstantCommandParsed distinguishes
         * from temperature floor option, so not reset if constant temperature
         * was actually set */
        if( StopCalc.TeFloor > 0. && !thermal.lgTemperatureConstantCommandParsed )
        {
                thermal.lgTemperatureConstant = false;
                thermal.ConstTemp = 0.;
        }

        lgAbort = false;

        /* reset some parameters needed for magnetic field */
        Magnetic_reinit();

        opac.stimax[0] = 0.;
        opac.stimax[1] = 0.;

        for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
                (*diatom)->H2_Reset();

        for( nelem=ipHYDROGEN; nelem < LIMELM; ++nelem )
        {
                /* these have one more ion than above */
                for( ion=0; ion<nelem+2; ++ion )
                {
                        /* zero out the source and sink arrays */
                        mole.source[nelem][ion] = saveMoleSource[nelem][ion];
                        mole.sink[nelem][ion] = saveMoleSink[nelem][ion];
                        /*>>chng 06 jun 27, move from here, where leak occurs, to
                        * above where xMoleChTrRate first created */
                        /*mole.xMoleChTrRate[nelem][ion] = (realnum*)MALLOC(sizeof(realnum)*(unsigned)(nelem+2));*/
                        for( ion2=0; ion2<nelem+2; ++ion2 )
                        {
                                mole.xMoleChTrRate[nelem][ion][ion2] = SaveMoleChTrRate[nelem][ion][ion2];
                        }
                }
        }

        /* reset molecular abundances */
        for( i=0; i<mole_global.num_calc; i++) 
        {
                mole.species[i].den = den_save[i];
        }
        dense.updateXMolecules();
        hmi.H2_total = findspecieslocal("H2")->den + findspecieslocal("H2*")->den;
        hmi.HD_total = findspecieslocal("HD")->den;
        /*fprintf(ioQQQ," IterRestar sets H2 total to %.2e\n",hmi.H2_total );*/
        h2.ortho_density = ortho_save;
        h2.para_density = para_save;
        {
                hmi.H2_total_f = (realnum)hmi.H2_total;
                h2.ortho_density_f = (realnum)h2.ortho_density;
                h2.para_density_f = (realnum)h2.para_density;
        }

        /* zero out the column densities */
        for( i=0; i < NCOLD; i++ )
        {
                colden.colden[i] = 0.;
        }
        colden.He123S = 0.;
        colden.coldenH2_ov_vel = 0.;

        for( i=0; i < mole_global.num_calc; i++ )
        {
                /* column densities */
                mole.species[i].column = 0.;
                /* largest fraction of atoms in molecules */
                mole.species[i].xFracLim = 0.;
                mole.species[i].nAtomLim = -1;
        }


        /* close out this iteration if dynamics or time dependence is enabled */
        if( dynamics.lgAdvection )
                DynaIterEnd();

        rfield.extin_mag_B_point = 0.;
        rfield.extin_mag_V_point = 0.;
        rfield.extin_mag_B_extended = 0.;
        rfield.extin_mag_V_extended = 0.;

        hmi.hmihet = hmihet_save;
        hmi.hmitot = hmitot_save;

        hmi.h2plus_heat = h2plus_heat_save;
        hmi.HeatH2Dish_used = HeatH2Dish_used_save;
        hmi.HeatH2Dexc_used = HeatH2Dexc_used_save;

        hmi.deriv_HeatH2Dexc_used = (realnum)deriv_HeatH2Dexc_used_save;
        hmi.H2_Solomon_dissoc_rate_used_H2g = H2_Solomon_dissoc_rate_used_H2g_save;
        hmi.H2_Solomon_dissoc_rate_used_H2s = H2_Solomon_dissoc_rate_used_H2s_save;
        hmi.H2_H2g_to_H2s_rate_used = H2_H2g_to_H2s_rate_used_save;
        hmi.H2_photodissoc_used_H2s = H2_photodissoc_used_H2s_save;
        hmi.H2_photodissoc_used_H2g = H2_photodissoc_used_H2g_save;

        hmi.UV_Cont_rel2_Draine_DB96_face = (realnum)UV_Cont_rel2_Draine_DB96_face;
        hmi.UV_Cont_rel2_Draine_DB96_depth = (realnum)UV_Cont_rel2_Draine_DB96_depth;
        hmi.UV_Cont_rel2_Habing_TH85_face = (realnum)UV_Cont_rel2_Habing_TH85_face;

        rfield.ipEnergyBremsThin = 0;
        rfield.EnergyBremsThin = 0.;
        rfield.lgUSphON = false;

        radius.glbdst = 0.;
        /* zone thickness, and next zone thickness */
        radius.drad = drSave;
        radius.drad_mid_zone = drSave/2.;
        radius.drNext = drNextSave;

        /* was min dr ever used? */
        radius.lgDrMinUsed = false;

        continuum.cn4861 = continuum.sv4861;
        continuum.cn1216 = continuum.sv1216;

        /* total luminosity */
        continuum.totlsv = continuum.TotalLumin;

        /* set debugger on now if NZONE desired is 0 */
        if( (trace.nznbug == 0 && iteration >= trace.npsbug) && trace.lgTrOvrd )
        {
                if( trace.nTrConvg==0 )
                {
                        trace.lgTrace = true;
                }
                else
                        /* trace convergence entered = but with negative flag = make positive,
                         * abs and not mult by -1 since may trigger more than one time */
                        trace.nTrConvg = abs( trace.nTrConvg );

                fprintf( ioQQQ, " IterRestart called.\n" );
        }
        else
        {
                trace.lgTrace = false;
        }

        /* zero secondary suprathermals variables for first ionization attem */
        for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
        {
                for( ion=0; ion<nelem+1; ++ion )
                {
                        secondaries.csupra[nelem][ion] = supsav[nelem][ion];
                }
        }
        secondaries.x12tot = secondaries.x12sav;
        secondaries.HeatEfficPrimary = secondaries.hetsav;
        secondaries.SecIon2PrimaryErg = secondaries.savefi;
        for( nelem=0; nelem<LIMELM; ++nelem)
        {
                for( ion=0; ion<nelem+1; ++ion )
                {
                        ionbal.CompRecoilHeatRate[nelem][ion] = ionbal.CompRecoilHeatRateSave[nelem][ion];
                        ionbal.CompRecoilIonRate[nelem][ion] = ionbal.CompRecoilIonRateSave[nelem][ion];
                }
        }

        wind.lgVelPos = true;
        wind.AccelMax = 0.;
        wind.AccelAver = 0.;
        wind.acldr = 0.;
        wind.windv = wind.windv0;

        thermal.nUnstable = 0;
        thermal.lgUnstable = false;

        pressure.pbeta = 0.;
        pressure.RadBetaMax = 0.;
        pressure.lgPradCap = false;
        pressure.lgPradDen = false;
        /* this flag will say we hit the sonic point */
        pressure.lgSonicPoint = false;
        pressure.lgStrongDLimbo = false;

        pressure.PresInteg = 0.;
        pressure.pinzon = 0.;
        pressure.PresIntegElecThin = 0.;

        pressure.RhoGravity_dark = 0.;
        pressure.RhoGravity_self = 0.;
        pressure.RhoGravity_external = 0.;
        pressure.RhoGravity = 0.;
        pressure.IntegRhoGravity = 0.;

        EdenChange( edsav );
        dense.EdenHCorr = edsav;
        dense.EdenHCorr_f = (realnum)dense.EdenHCorr;
        dense.EdenTrue = edsav;

        for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
        {
                /* reset molecular densities */
                dense.SetGasPhaseDensity( nelem, gas_phase_save[nelem] );
                dense.IonLow[nelem] = IonLowSave[nelem];
                dense.IonHigh[nelem] = IonHighSave[nelem];
        }
        /*fprintf(ioQQQ,"DEBUG IterRestart gas_phase set to save  hden %.4e\n",
                dense.gas_phase[0]);*/

        for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
        {
                for( long nelem=ipISO; nelem < LIMELM; ++nelem )
                {
                        iso_sp[ipISO][nelem].Reset();
                }
        }

        for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
        {
                for( nelem=ipISO; nelem < LIMELM; nelem++ )
                {
                        if( dense.lgElmtOn[nelem] )
                        {
                                for( n=ipH1s; n < iso_sp[ipISO][nelem].numLevels_max; n++ )
                                {
                                        iso_sp[ipISO][nelem].fb[n].ConOpacRatio = HOpacRatSav[ipISO][nelem][n];
                                        iso_sp[ipISO][nelem].st[n].Pop() = hnsav[ipISO][nelem][n];
                                }
                        }
                }
        }

        if( trace.lgTrace && trace.lgNeBug )
        {
                fprintf( ioQQQ, "     EDEN set to%12.4e by IterRestart.\n", 
                  dense.eden );
        }

        for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
        {
                for( ion=0; ion < (nelem + 2); ion++ )
                {
                        dense.xIonDense[nelem][ion] = xIonFsave[nelem][ion];
                }
                for( i=0; i < LIMELM; i++ )
                {
                        thermal.heating[nelem][i] = HeatSave[nelem][i];
                }
        }

        deut.xIonDense[0] = deutDenseSave0;
        deut.xIonDense[1] = deutDenseSave1;

        GrainRestartIter();

        rfield.resetCoarseTransCoef();

        /* continuum was saved in flux_total_incident */
        for( i=0; i < rfield.nupper; i++ )
        {
                /* time-constant part of beamed continuum */
                rfield.flux_beam_const[i] = rfield.flux_beam_const_save[i];
                /* continuum flux_time_beam_save has the initial value of the
                 * time-dependent beamed continuum */
                rfield.flux_beam_time[i] = rfield.flux_time_beam_save[i]*
                        rfield.time_continuum_scale;

                if( cosmology.lgDo )
                {
                        double slope_ratio;
                        double fac_old = TE1RYD*rfield.anu[i]/(CMB_TEMP*(1. + cosmology.redshift_start  ));
                        double fac_new = TE1RYD*rfield.anu[i]/(CMB_TEMP*(1. + cosmology.redshift_current));

                        if( fac_old > log10(DBL_MAX) )
                        {
                                slope_ratio = 0.;
                        }
                        else if( fac_new > log10(DBL_MAX) )
                        {
                                slope_ratio = 0.;
                        }
                        else if( fac_old > 1.e-5 )
                        {
                                slope_ratio = (exp(fac_old) - 1.)/(exp(fac_new) - 1.);
                        }
                        else
                        {
                                slope_ratio = (fac_old*(1. + fac_old/2.))/(fac_new*(1. + fac_new/2.));
                        }

                        rfield.flux_isotropic[i] = rfield.flux_isotropic_save[i] * 
                                pow( ((realnum)1.f + cosmology.redshift_current)/((realnum)1.f + cosmology.redshift_start), (realnum)4.f ) * (realnum)slope_ratio;
                }
                else
                {
                        /* the isotropic continuum source is time steady */
                        rfield.flux_isotropic[i] = rfield.flux_isotropic_save[i];
                }

                rfield.flux[0][i] = rfield.flux_isotropic[i] + rfield.flux_beam_time[i] +
                        rfield.flux_beam_const[i];

                rfield.SummedDif[i] = rfield.SummedDifSave[i];
                rfield.SummedCon[i] = rfield.flux[0][i] + rfield.SummedDif[i];

                rfield.OccNumbIncidCont[i] = rfield.flux[0][i]*rfield.convoc[i];
                rfield.otscon[i] = rfield.otssav[i][0];
                rfield.otslin[i] = rfield.otssav[i][1];
                rfield.ConInterOut[i] = 0.;
                rfield.OccNumbBremsCont[i] = 0.;
                rfield.OccNumbDiffCont[i] = 0.;
                rfield.OccNumbContEmitOut[i] = 0.;
                rfield.outlin[0][i] = 0.;
                rfield.outlin_noplot[i] = 0.;
                rfield.ConOTS_local_OTS_rate[i] = 0.;
                rfield.ConOTS_local_photons[i] = 0.;
                rfield.DiffuseEscape[i] = 0.;

                opac.opacity_abs[i] = opac.opacity_abs_savzon1[i];
                opac.OldOpacSave[i] = opac.opacity_abs_savzon1[i];
                opac.opacity_sct[i] = opac.opacity_sct_savzon1[i];
                opac.albedo[i] = 
                        opac.opacity_sct[i]/MAX2(1e-30,opac.opacity_sct[i] + opac.opacity_abs[i]);
                opac.tmn[i] = 1.;
                /* >>chng 99 dec 04, having exactly 1 for first zone caused discontinuity
                 * for heating in very high T models in func_map.in.  zone 1 and 2 were 20% different,
                 * since tau in is 1e-20, e2 is 0.9999, and so some H ots
                opac.ExpmTau[i] = 1.;
                opac.E2TauAbsFace[i] = 1.;*/
                /* >>chng 99 dec 04, having exactly 1 for first zone caused discontinuity
                 * for heating in very high T models in func_map.in.  zone 1 and 2 were 20% different,
                 * since tau in is 1e-20, e2 is 0.9999, and so some H ots
                 * these were not here at all*/
                /* attenuation of flux by optical depths IN THIS ZONE 
                 * DirectionalCosin is 1/COS(theta), is usually 1, reset with illuminate command,
                 * option for illumination of slab at an angle */
                opac.ExpZone[i] = sexp(opac.opacity_abs[i]*radius.drad/2.*geometry.DirectionalCosin);

                /* e(-tau) in inward direction, up to illuminated face */
                opac.ExpmTau[i] = (realnum)opac.ExpZone[i];

                /* e2(tau) in inward direction, up to illuminated face */
                opac.E2TauAbsFace[i] = (realnum)e2(opac.TauAbsFace[i]);
                rfield.reflin[0][i] = 0.;
                rfield.ConEmitReflec[0][i] = 0.;
                rfield.ConEmitOut[0][i] = 0.;
                rfield.ConRefIncid[0][i] = 0.;

                /* escape in the outward direction
                 * on second and later iterations define outward E2 */
                if( iteration > 1 )
                {
                        /* e2 from current position to outer edge of shell */
                        realnum tau = MAX2(SMALLFLOAT , opac.TauAbsTotal[i] - opac.TauAbsFace[i] );
                        opac.E2TauAbsOut[i] = (realnum)e2( tau );
                        ASSERT( opac.E2TauAbsOut[i]>=0. && opac.E2TauAbsOut[i]<=1. );
                }
                else
                        opac.E2TauAbsOut[i] = 1.;

        }

        /* update continuum */
        RT_OTS_Update(&SumOTS);

        thermal.FreeFreeTotHeat = 0.;
        atoms.p2nit = p2nit;
        atoms.d5200r = d5200r;

        if( called.lgTalk )
        {
                fprintf( ioQQQ, "\f\n          Start Iteration Number %ld   %75.75s\n\n\n", 
                  iteration, input.chTitle );
        }

        /* reset variable to do with FeII atom, in FeIILevelPops */
        FeIIReset();
        ASSERT(lgElemsConserved());
        return;
}

/* do some work with ending the iteration */
void IterEnd(void)
{

        DEBUG_ENTRY( "IterEnd()" );

        if( lgAbort )
        {
                return;
        }

        /* give indication of geometry */
        double fac = radius.depth/radius.rinner;
        if( fac < 0.1 )
        {
                geometry.lgGeoPP = true;
        }
        else
        {
                geometry.lgGeoPP = false;
        }

        if( iteration > dynamics.n_initial_relax && dynamics.lgTimeDependentStatic )
        {
                // report cumulative lines per unit mass rather than flux (per unit
                // area), so total is meaningful when density isn't constant

                double cumulative_factor = (dynamics.timestep/
                                                                                         colden.TotMassColl);
                // save cumulative lines
                for( long n=0; n<LineSave.nsum; ++n )
                {
                        for( long nEmType=0; nEmType<2; ++nEmType )
                        {
                                LineSv[n].SumLine[nEmType+2] += (realnum) LineSv[n].SumLine[nEmType]*cumulative_factor;
                        }
                }
                // save cumulative continua
                for( long n=0; n<rfield.nflux; ++n)
                {
                        
                        rfield.flux[1][n] += (realnum) rfield.flux[0][n]*cumulative_factor;
                        rfield.ConEmitReflec[1][n] += (realnum) rfield.ConEmitReflec[0][n]*cumulative_factor;
                        rfield.ConEmitOut[1][n] += (realnum) rfield.ConEmitOut[0][n]*cumulative_factor;
                        rfield.ConRefIncid[1][n] += (realnum) rfield.ConRefIncid[0][n]*cumulative_factor;
                        rfield.flux_total_incident[1][n] += (realnum) rfield.flux_total_incident[0][n]*cumulative_factor;
                        rfield.reflin[1][n] += (realnum) rfield.reflin[0][n]*cumulative_factor;
                        rfield.outlin[1][n] += (realnum) rfield.outlin[0][n]*cumulative_factor;
                }
        }
        
        
        /* save previous iteration's results */
        struc.nzonePreviousIteration = nzone;
        for( long i=0; i<struc.nzonePreviousIteration; ++i )
        {
                struc.depth_last[i] = struc.depth[i];
                struc.drad_last[i] = struc.drad[i];
        }

        /* all continue were attenuated in last zone in radius_increment to represent the intensity
         * in the middle of the next zone - this was too much since we now want
         * intensity emergent from outer edge of last zone */
        for( long i=0; i < rfield.nflux; i++ )
        {
                {
                        enum{DEBUG_LOC=false};
                        if( DEBUG_LOC)
                        {
                                fprintf(ioQQQ,"i=%li opac %.2e \n", i, 
                                        (double)opac.opacity_abs[i]*radius.drad_x_fillfac/2.*(double)geometry.DirectionalCosin );
                        }
                }
                double tau = opac.opacity_abs[i]*radius.drad_x_fillfac/2.*geometry.DirectionalCosin;
                ASSERT( tau > 0. );
                fac = sexp( tau );

                /* >>chng 02 dec 14, add first test to see whether product of ratio is within
                 * range of a float - ConInterOut can be finite and fac just above a small float,
                 * so ratio exceeds largest size of a float */
                /*if( fac > SMALLFLOAT )*/
                if( (realnum)(fac/SDIV(rfield.ConInterOut[i]))>SMALLFLOAT && fac > SMALLFLOAT )
                {
                        rfield.ConInterOut[i] /= (realnum)fac;
                        rfield.outlin[0][i] /= (realnum)fac;
                        rfield.outlin_noplot[i] /= (realnum)fac;
                }
        }

        /* remember thickness of previous iteration */
        radius.StopThickness[iteration-1] = radius.depth;
        return;
}