cloudy  trunk
radius_next.cpp
Go to the documentation of this file.
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
00002  * others.  For conditions of distribution and use see copyright notice in license.txt */
00003 /*radius_next use adaptive logic to find next zone thickness */
00004 /*ContRate called by radius_next to find energy of maximum continuum-gas interaction */
00005 /*GrainRateDr called by radius_next to find grain heating rate dr */
00014 #include "cddefines.h"
00015 #include "lines_service.h"
00016 #include "iso.h"
00017 #include "geometry.h"
00018 #include "h2.h"
00019 #include "mole.h"
00020 #include "hyperfine.h"
00021 #include "opacity.h"
00022 #include "dense.h"
00023 #include "heavy.h"
00024 #include "grainvar.h"
00025 #include "elementnames.h"
00026 #include "conv.h"
00027 #include "rfield.h"
00028 #include "dynamics.h"
00029 #include "thermal.h"
00030 #include "hmi.h"
00031 #include "coolheavy.h"
00032 #include "timesc.h"
00033 #include "doppvel.h"
00034 #include "stopcalc.h"
00035 #include "colden.h"
00036 #include "phycon.h"
00037 #include "rt.h"
00038 #include "trace.h"
00039 #include "wind.h"
00040 #include "punch.h"
00041 #include "taulines.h"
00042 #include "pressure.h"
00043 #include "iterations.h"
00044 #include "struc.h"
00045 #include "radius.h"
00046 
00047 #if 0
00048 /*ChkRate called by radius_next to check how rates of destruction of various species changes */
00049 STATIC void ChkRate(
00050           /* element number on physical scale */
00051           long int nelem, 
00052           /* change in destruction rate */
00053           double *dDestRate, 
00054           /* old and new destruction rates */
00055           double *DestRateOld,
00056           double *DestRateNew,
00057           /* stage of ionization on the physical scale */
00058           long int *istage);
00059 #endif
00060 
00061 /*ContRate called by radius_next to find energy of maximum continuum-gas interaction */
00062 STATIC void ContRate(double *freqm, 
00063   double *opacm);
00064 
00065 /*GrainRateDr called by radius_next to find grain heating rate dr */
00066 STATIC void GrainRateDr(double *grfreqm, 
00067   double *gropacm);
00068 
00069 /*radius_next use adaptive logic to find next zone thickness 
00070  * return 0 if ok, 1 for abort */
00071 int radius_next(void)
00072 {
00073         char chLbl[11];
00074         bool lgDoPun;
00075 
00076         long int level , ipStrong;
00077 
00078         double thickness_total , drThickness , DepthToGo , AV_to_go;
00079         int mole_dr_change;
00080 
00081         double DrGrainHeat, 
00082           GlobDr, 
00083           SaveOHIIoHI, 
00084           SpecDr, 
00085           Strong, 
00086           TauDTau, 
00087           TauInwd, 
00088           drSolomon_BigH2 ,
00089           dEfrac, 
00090           dHdStep, 
00091           dRTaue, 
00092           dTdStep, 
00093           dnew, 
00094           drConPres, 
00095           drH2_heat_cool = 0. ,
00096           dH2_heat_cool = 0.,
00097           drH2_abund = 0. ,
00098           dr_mole_abund = 0.,
00099           dH2_abund=0.,
00100           dCO_abund=0.,
00101           drDepth, 
00102           drDest, 
00103           drEfrac, 
00104           drFail, 
00105           drFluc, 
00106           drHMase, 
00107           drHe1ovHe2,
00108           drHion, 
00109           drInter, 
00110           drLeiden_hack ,
00111           drLineHeat, 
00112           drSphere, 
00113           drTab, 
00114           drdHdStep, 
00115           drdTdStep, 
00116           drThermalFront ,
00117           drmax, 
00118           dVelRelative,
00119           fac2, 
00120           factor, 
00121           freqm, 
00122           grfreqm=0., 
00123           gropacm=0., 
00124           hdnew, 
00125           opacm, 
00126           recom_dr_last_iter ,
00127           OldDR ,
00128           winddr, 
00129           WindAccelDR,
00130           x;
00131 
00132         double change_heavy_frac=-1. , change_heavy_frac_big , dr_change_heavy ,
00133                 frac_change_heavy_frac_big, Efrac_old , Efrac_new;
00134         long int ichange_heavy_nelem=-1 , nelem , ion , ichange_heavy_ion=-1;
00135 
00136         static double OHIIoHI, 
00137           OldHeat = 0., 
00138           OldTe = 0.,
00139           OlddTdStep = 0.,
00140           OldH2Abund=0.,
00141           OldWindVelocity=0.,
00142           Old_He_atom_ov_ion = 0,
00143           Old_H2_heat_cool;
00144         static long int iteration_last=-1;
00145 
00146         static double BigRadius = 1e30;
00147         bool lgFirstCall;
00148 
00149         DEBUG_ENTRY( "radius_next()" );
00150 
00151 
00152         /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
00153          *
00154          * changes in logic
00155          * 95 oct 19, drSphere now 3% of radius, was 2%, fewer zone
00156          * 95 oct 19, subtracted grain opacity from total opacity used
00157          * to get thickness in routine ContRate
00158          *
00159          *>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
00160          *
00161          * free statement labels >= 13
00162          *
00163          *-----------------------------------------------------------------------
00164          *
00165          * this sub determines the thickness of the next zone
00166          * if is called one time for each zone
00167          * flag lgNxtDROff is true if this is initialization of radius_next,
00168          * is false if we are to use logic to find dr
00169          *
00170          *----------------------------------------------------------------------- */
00171 
00172         /* >>chng 03 sep 21 - decide whether this is the first call */
00173         if( (iteration != iteration_last) && (nzone==0) )
00174         {
00175                 /* this is the first call in this iteration */
00176                 iteration_last = iteration;
00177                 lgFirstCall = true;
00178         }
00179         else
00180                 lgFirstCall = false;
00181 
00182         if( trace.lgTrace )
00183         {
00184                 fprintf( ioQQQ, "   radius_next called\n" );
00185         }
00186 
00187         /* save current dr */
00188         OldDR = radius.drad;
00189 
00190         /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
00191         /*  1    '' radius_next keys from change in H ionization'',e11.3)')
00192          * check on change in hydrogen ionizaiton */
00193         if( lgFirstCall )
00194         {
00195                 if( dense.xIonDense[ipHYDROGEN][1] > 0. && dense.xIonDense[ipHYDROGEN][0] > 0. )
00196                 {
00197                         OHIIoHI = dense.xIonDense[ipHYDROGEN][1]/dense.xIonDense[ipHYDROGEN][0];
00198                 }
00199                 else
00200                 {
00201                         OHIIoHI = 0.;
00202                 }
00203                 SaveOHIIoHI = OHIIoHI;
00204                 drHion = BigRadius;
00205                 /* else if(hii.gt.0. .and. hi.gt.0. .and. OHIIoHI.gt.0. ) then
00206                  * >>chng 97 jul 9, for deep in PDR vastly now ionz H slowed down works */
00207         }
00208 
00209         else if( (dense.xIonDense[ipHYDROGEN][1] > 0. && dense.xIonDense[ipHYDROGEN][0] > 0.) && OHIIoHI > 1e-15 )
00210         {
00211                 double atomic_frac = (dense.xIonDense[ipHYDROGEN][1]/dense.xIonDense[ipHYDROGEN][0]);
00212                 /* ratio of current HII/HI to old value - < 1 when becoming more neutral */
00213                 /* this is now change in atomic fraction */
00214                 x = 1. - atomic_frac /OHIIoHI;
00215                 if( atomic_frac > 0.05 && atomic_frac < 0.9 )
00216                 {
00217                         /* >>chng 96 feb 23 from 0.7 to 0.3 because punched through H i-front
00218                          * >>chng 97 may 5, from 0.3 to 0.2 since punched through H i-front */
00219                         /* >>chng 02 jun 05 from 0.2 to 0.05 poorly resolved i-front, also added two-branch logic*/
00220                         drHion = radius.drad*MAX2( 0.2 , 0.05/MAX2(1e-10,x) );
00221                 }
00222                 else if( x > 0. )
00223                 {
00224                         /* >>chng 96 feb 23 from 0.7 to 0.3 because punched through H i-front
00225                          * >>chng 97 may 5, from 0.3 to 0.2 since punched through H i-front */
00226                         drHion = radius.drad*MAX2( 0.2 , 0.2/MAX2(1e-10,x) );
00227                 }
00228                 else
00229                 {
00230                         drHion = BigRadius;
00231                 }
00232                 SaveOHIIoHI = OHIIoHI;
00233                 OHIIoHI = dense.xIonDense[ipHYDROGEN][1]/dense.xIonDense[ipHYDROGEN][0];
00234         }
00235 
00236         else
00237         {
00238                 SaveOHIIoHI = OHIIoHI;
00239                 if( dense.xIonDense[ipHYDROGEN][1] > 0. && dense.xIonDense[ipHYDROGEN][0] > 0. )
00240                 {
00241                         OHIIoHI = dense.xIonDense[ipHYDROGEN][1]/dense.xIonDense[ipHYDROGEN][0];
00242                 }
00243                 else
00244                 {
00245                         OHIIoHI = 0.;
00246                 }
00247                 drHion = BigRadius;
00248         }
00249 
00250         /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
00251         /* "radius_next keys from H maser, dt, ij=" possible hydrogen maser action */
00252         if( rt.dTauMase < -0.01 )
00253         {
00254                 /* maser so powerful, must limit inc in tay
00255                  * >>chng 96 aug 08, from 0.3 to 0.1 due to large maser crash */
00256                 drHMase = radius.drad*MAX2(0.1,-0.2/rt.dTauMase);
00257         }
00258         else
00259         {
00260                 drHMase = BigRadius;
00261         }
00262 
00263         /* >>chng 03 nov 09, try doing he in following, not above */
00264         Old_He_atom_ov_ion = 0.;
00265         drHe1ovHe2 = BigRadius;
00266 
00267         /* check on N0 - 1 ionization changes,
00268          * >>chng 03 jun 06, add this test due to smashing into H ifront in blr89.in
00269          */
00270         dr_change_heavy = BigRadius;
00271         change_heavy_frac_big = -1.;
00272         frac_change_heavy_frac_big = -1.;
00273         /* >chng 03 jun 09, from 0.05 to 0.1, initial tests with zoning */
00274 #       define CHANGE_ION_HEAV  0.2f
00275 #       define CHANGE_ION_HHE   0.15f
00276         if( nzone > 4 )
00277         {
00278                 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
00279                 {
00280                         if( dense.lgElmtOn[nelem] )
00281                         {
00282                                 realnum change;
00283                                 /* this is the limit to the ionization we will check -
00284                                  * also used in prt_comment to check on whether oscillations occurred */
00285                                 realnum frac_limit;
00286                                 if( nelem<=ipHELIUM || nelem==ipCARBON )
00287                                 {
00288                                         frac_limit = 1e-4f;
00289                                         change = CHANGE_ION_HHE;
00290                                 }
00291                                 else
00292                                 {
00293                                         /* this var is used to print warnings,
00294                                          * make sure we converge below it */
00295                                         frac_limit = struc.dr_ionfrac_limit/2.f;
00296                                         change = CHANGE_ION_HEAV;
00297                                 }
00298                                 /* >>chng 03 dec 09, go up to full range of ion, not just =2 */
00299                                 for( ion=0; ion<=nelem+1; ++ion )
00300                                 {
00301                                         realnum abundnew = dense.xIonDense[nelem][ion]/SDIV( dense.gas_phase[nelem]);
00302                                         realnum abundold = struc.xIonDense[nelem][ion][nzone-3]/SDIV( struc.gas_phase[nelem][nzone-3]);
00303                                         if( abundnew > frac_limit && abundold > frac_limit )
00304                                         {
00305                                                 /* NB must make sure this test is not done when nzone-x <0 */
00306                                                 /* -2 because previous zone, and nzone is off by one (on physical, not C, scale) */
00307                                                 /*realnum abundold = struc.xIonDense[nelem][ion][nzone-3]/SDIV( struc.gas_phase[nelem][nzone-3]);*/
00308                                                 realnum abundolder = struc.xIonDense[nelem][ion][nzone-4]/SDIV( struc.gas_phase[nelem][nzone-4]);
00309                                                 realnum abundoldest = struc.xIonDense[nelem][ion][nzone-5]/SDIV( struc.gas_phase[nelem][nzone-5]);
00310                                                 /* this is fractional change */
00311                                                 /* >>chng 04 feb 28, take min of old and new abund, to try to pick up
00312                                                  * rapid changing Ca+ sooner */
00313                                                 change_heavy_frac = fabs(abundnew-abundold)/MIN2(abundold,abundnew);
00314                                                 /* want fractional change to be less than this factor */
00315                                                 if( (change_heavy_frac > change) && (change_heavy_frac > change_heavy_frac_big) &&
00316                                                         /* >>chng 03 dec 07, add test that abund is not oscillating */
00317                                                         /* also test that abundance is increasing - we are headed into a front */
00318                                                         ((abundnew-abundold)>0.)   && 
00319                                                         ((abundold-abundolder)>0.) && 
00320                                                         ((abundolder-abundoldest)>0.) )
00321                                                 {
00322                                                         ichange_heavy_nelem = nelem;
00323                                                         ichange_heavy_ion = ion;
00324                                                         change_heavy_frac_big = change_heavy_frac;
00325                                                         frac_change_heavy_frac_big = abundnew;
00326                                                         /* >>chng 03 dec 06, from min of 0.5 to min of 0.25, crash into He i-front 
00327                                                          * in hizqso.in */
00328                                                         /* >>chng 04 mar 03, min had become 0.1, forced oscillations in nova.in
00329                                                          * in silicon, zoning changed greatly, causing change in diffuse lin
00330                                                          * pumping.  put back to 0.25 */
00331                                                         dr_change_heavy = radius.drad * MAX2(0.25, change / change_heavy_frac );
00332                                                 }
00333                                         }
00334                                 }
00335                         }
00336                 }
00337         }
00338 
00339         /* if Leiden hacks are on then use increase in dust extinction as control 
00340          * >>chng 05 aug 13, add this */
00341         if(!co.lgUMISTrates)
00342         {
00343                 /* Draine field is only defined over narrow range in FUV - must not let change
00344                  * in extinction become too large - 
00345                  * prefactor is change in optical depth */
00346                 drLeiden_hack = MAX2( 0.02 , 0.05*rfield.extin_mag_V_point) / SDIV( geometry.FillFac * rfield.opac_mag_V_point );
00347         }
00348         else
00349         {
00350                 drLeiden_hack = BigRadius;
00351         }
00352         /* >>chng 04 feb 15, kill this block - not used */
00353 
00354         /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
00355         /* check how heating is changing
00356          * '' radius_next keys from change in heating; current, delta='', */
00357         if( nzone <= 1 || thermal.lgTemperatureConstant )
00358         {
00359                 drdHdStep = BigRadius;
00360                 dHdStep = FLT_MAX;
00361         }
00362         else
00363         {
00364                 dHdStep = fabs(thermal.htot-OldHeat)/thermal.htot;
00365                 if( dHdStep > 0. )
00366                 {
00367                         if( dense.gas_phase[ipHYDROGEN] >= 1e13 )
00368                         {
00369                                 /* drdHdStep = drad * MAX( 0.8 , 0.05/dHdStep ) */
00370                                 drdHdStep = radius.drad*MAX2(0.8,0.075/dHdStep);
00371                         }
00372                         else if( dense.gas_phase[ipHYDROGEN] >= 1e11 )
00373                         {
00374                                 /* drdHdStep = drad * MAX( 0.8 , 0.075/dHdStep ) */
00375                                 drdHdStep = radius.drad*MAX2(0.8,0.100/dHdStep);
00376                         }
00377                         else
00378                         {
00379                                 /* changed from .15 to .12 for outer edge of coolhii too steep dT
00380                                  * changed to .10 for symmetry, big change in some rates, 95nov14
00381                                  * changed from .10 to .125 since parispn seemed to waste zones
00382                                  * >>chng 96 may 21, from .125 to .15 since pn's still waste zones */
00383                                 drdHdStep = radius.drad*MAX2(0.8,0.15/dHdStep);
00384                         }
00385                 }
00386                 else
00387                 {
00388                         drdHdStep = BigRadius;
00389                 }
00390         }
00391         OldHeat = thermal.htot;
00392 
00393         /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
00394         /* pressure due to incident continuum if in eos */
00395         if( strcmp(dense.chDenseLaw,"CPRE") == 0 && pressure.lgContRadPresOn )
00396         {
00397                 if( nzone > 2 && pressure.pinzon > 0. )
00398                 {
00399                         /* pinzon is pressrue from acceleration onto previos zone
00400                          * want this less than some fraction of total pressure */
00401                         /* >>chng 06 feb 01, change from init pres to current total pressure
00402                          * in const press high U ulirgs current pressure may be quite larger
00403                          * than init pressure due to continuum absorption */
00404                         drConPres = 0.05*pressure.PresTotlCurr/(wind.AccelTot*
00405                           dense.xMassDensity*geometry.FillFac);
00406                 }
00407                 else
00408                 {
00409                         drConPres = BigRadius;
00410                 }
00411         }
00412         else
00413         {
00414                 drConPres = BigRadius;
00415         }
00416 
00417         /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
00418         /* check how temperature is changing
00419          *  1    '' radius_next keys from change in temperature; current, delta='', */
00420         if( nzone <= 1 )
00421         {
00422                 drdTdStep = BigRadius;
00423                 dTdStep = FLT_MAX;
00424                 OlddTdStep = dTdStep;
00425         }
00426         else
00427         {
00428                 /* change in temperature; current=      */
00429                 dTdStep = (phycon.te-OldTe)/phycon.te;
00430                 /* >>chng 02 dec 08, desired change in temperature per zone must not
00431                  * be smaller than allower error in temperature.  For now use relative error
00432                  * in heating - - cooling balance.  Better would be to also use c-h deriv wrt t
00433                  * to get slope */
00434                 x = conv.HeatCoolRelErrorAllowed*2.;
00435                 x = MAX2( 0.01 , x ); 
00436                 x = MIN2( 0.05 , x );
00437                 /* >>chng 02 dec 11 rjrw change back to 0.03 -- improve dynamics.dRad criterion instead */
00438                 x = 0.03;
00439                 /* >>chng 02 dec 07, do not do this if there is mild te jitter, 
00440                  * so that dT is changing sign - this happens
00441                  * in ism.ini, very stable temperature with slight noise up and down */
00442                 if( dTdStep*OlddTdStep > 0. )
00443                 {
00444                         /* >>chng 96 may 30, variable depending on temp
00445                          * >>chng 96 may 31, allow 0.7 smaller, was 0.8
00446                          * >>chng 97 may 05, from 0.7 to 0.5 stop from punching through thermal front
00447                          * >>chng 04 mar 30, from 0.7 to 0.5 stop from punching through thermal front,
00448                          * for some reason factor was 0.7, not 0.5 as per previous change */
00449                         double absdt = fabs(dTdStep);
00450                         drdTdStep = radius.drad*MAX2(0.5,x/absdt);
00451                 }
00452                 else
00453                 {
00454                         drdTdStep = BigRadius;
00455                 }
00456         }
00457         OlddTdStep = dTdStep;
00458         OldTe = phycon.te;
00459 
00460         /* >>chng 02 oct 06, only check on opacity - interaction if not
00461          * constant temperature - there were constant temperature models that
00462          * extended to infinity but hung with last few photons and this test.
00463          * better to ignore this test which is really for thermal models */
00464         /* >>chng 06 mar 20, do not call if recombination logic in place */
00465         if( !thermal.lgTemperatureConstant && !dynamics.lgRecom )
00466         {
00467                 /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
00468                 /* find freq where opacity largest and interaction rate is fastest
00469                 * "cont inter nu=%10.2e opac=%10.2e\n" */
00470                 ContRate(&freqm,&opacm);
00471 
00472                 /* use optical depth at max interaction energy
00473                 * >>chng 96 jun 06 was drChange=0.15 changed to 0.3 for high Z models
00474                 * taking too many zones
00475                 * drInter = drChange / MAX(1e-30,opacm*FillFac) */
00476 
00477                 drInter = 0.3/MAX2(1e-30,opacm*geometry.FillFac*geometry.DirectionalCosin);
00478         }
00479         else
00480         {
00481                 drInter = BigRadius;
00482                 freqm = 0.;
00483                 opacm = 1.;
00484         }
00485 
00486 #       if 0
00487         /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
00488         /* check on changes in destruction rates for various atoms */
00489         ChkRate(ipCARBON,&dDRCarbon,&DestOldCarb, &DestNewCarb,&icarstag);
00490         ChkRate(ipNITROGEN,&dDRNitrogen,&DestOldNit, &DestNewNit,&initstag);
00491         ChkRate(ipOXYGEN,&dDROxygen,&DestOldOxy, &DestNewOxy,&ioxystag);
00492         ChkRate(ipIRON,&dDRIron,&DestOldIron, &DestNewIron,&iironstag);
00493 
00494         dDestRate = MAX4(dDROxygen,dDRIron,dDRCarbon,dDRNitrogen);
00495 
00496         if( fp_equal( dDRCarbon, dDestRate ) )
00497         {
00498                 dDestRate = dDRCarbon;
00499                 DestOld = DestOldCarb;
00500                 DestNew = DestNewCarb;
00501                 istage = icarstag;
00502                 strcpy( chDestAtom, "Carbon  " );
00503         }
00504 
00505         else if( fp_equal( dDRNitrogen, dDestRate ) )
00506         {
00507                 dDestRate = dDRNitrogen;
00508                 DestOld = DestOldNit;
00509                 DestNew = DestNewNit;
00510                 istage = initstag;
00511                 strcpy( chDestAtom, "Nitrogen" );
00512         }
00513 
00514         else if( fp_equal( dDROxygen, dDestRate ) )
00515         {
00516                 dDestRate = dDROxygen;
00517                 DestOld = DestOldOxy;
00518                 DestNew = DestNewOxy;
00519                 strcpy( chDestAtom, "Oxygen  " );
00520                 istage = ioxystag;
00521         }
00522 
00523         else if( fp_equal( dDRIron, dDestRate ) )
00524         {
00525                 dDestRate = dDRIron;
00526                 DestOld = DestOldIron;
00527                 DestNew = DestNewIron;
00528                 istage = iironstag;
00529                 strcpy( chDestAtom, "Iron    " );
00530         }
00531 
00532         else
00533         {
00534                 fprintf( ioQQQ, " insanity following ChkRate\n" );
00535                 ShowMe();
00536                 cdEXIT(EXIT_FAILURE);
00537         }
00538 
00539         /*  radius_next keys from change in dest rates, atom= */
00540         if( dDestRate > 0. )
00541         {
00542                 /* if( te.gt.40 000. ) then
00543                  * added different accuracy for hot gas since tend to jump over
00544                  * big te range for small change in heat (intrinsically unstable)
00545                  * drDest = drad * MAX( 0.5 , 0.10/dDestRate )
00546                  * else
00547                  * was drChange, changed to .15 for parishii go through HeII=HeI I front
00548                  * drDest = drad * MAX( 0.5 , 0.15/dDestRate )
00549                  * >>chng 95 dec 18 from above to below to stop oscillations
00550                  * >>chng 95 dec 27 from min of .5 to .75 to stop zone size from changing rapidly
00551                  * drDest = drad * MAX( 0.8 , 0.20 /dDestRate )
00552                  * >>chng 96 jan 14 from .2 to .25 to use less zones
00553                  * >>chng 96 may 30 from min of 0.8 to 0.5 to prevent crashing into He i-front */
00554                 drDest = radius.drad*MAX2(0.5,0.20/dDestRate);
00555                 /* endif */
00556         }
00557         else
00558         {
00559                 drDest = BigRadius;
00560         }
00561 #       endif
00562         drDest = BigRadius;/*03 dec 12 is this needed? */
00563 
00564         /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
00565         /* check whether change in wind velocity constrains DRAD 
00566          * WJH 22 May 2004: disable when we are near the sonic point since
00567          * the velocity may be jumping all over the place but we just want
00568          * to push through it as quickly as we can */
00569         if( wind.windv!=0. && !pressure.lgSonicPoint && !pressure.lgStrongDLimbo )
00570         {
00571                 double v = fabs(wind.windv);
00572                 /* this is relative change in velocity */
00573                 dVelRelative = fabs(wind.windv-OldWindVelocity)/
00574                         MAX2(v,0.1*timesc.sound_speed_isothermal);
00575 
00576 #               define WIND_CHNG_VELOCITY_RELATIVE      0.01
00577                 /*fprintf(ioQQQ,"DEBUG rad %.3f vel %.2e dVelRelative %.2e", 
00578                         log10(radius.Radius) , wind.windv ,  dVelRelative );*/
00579 
00580                 /* do not use this on first zone since do not have old velocity */
00581                 if( dVelRelative > WIND_CHNG_VELOCITY_RELATIVE  && nzone>1 )
00582                 {
00583                         /* factor less than one if change larger than WIND_CHNG_VELOCITY_RELATIVE */
00584                         double factor = WIND_CHNG_VELOCITY_RELATIVE / dVelRelative;
00585                         /*fprintf(ioQQQ," DEBUG factor %.2e", factor );*/
00586                         factor = MAX2(0.8 , factor );
00587                         winddr = radius.drad * factor;
00588                 }
00589                 else
00590                 {
00591                         winddr = BigRadius;
00592                 }
00593                 /*fprintf(ioQQQ," DEBUG \n" );*/
00594 
00595                 WindAccelDR = BigRadius;
00596 
00597                 /* set dr from advective term,
00598                  * the 1/20 came from looking at one set of structure plots */
00599                 if( dynamics.lgAdvection && iteration > dynamics.n_initial_relax+1)
00600                 {
00601                         /* >>chng 02 dec 11, from 5 to 20 */
00602                         winddr = MIN2( winddr , dynamics.dRad / 20. );
00603                         /*>>chng 04 oct 06, set dVelRelative to dynamics.dRad since dVelRelative is printed as part
00604                          * of reason for choosing this criteria, want to reflect valid reason. */
00605                         dVelRelative = dynamics.dRad;
00606                 }
00607         }
00608         else
00609         {
00610                 winddr = BigRadius;
00611                 WindAccelDR = BigRadius;
00612                 dVelRelative = 0.;
00613         }
00614         /* remember old velocity */
00615         OldWindVelocity = wind.windv;
00616 
00617         /* inside out globule */
00618         if( strcmp(dense.chDenseLaw,"GLOB") == 0 )
00619         {
00620 #               define  DNGLOB  0.10
00621                 if( radius.glbdst < 0. )
00622                 {
00623                         fprintf( ioQQQ, " Globule distance is negative, internal overflow has occured,  sorry.\n" );
00624                         fprintf( ioQQQ, " This is routine radius_next, GLBDST is%10.2e\n", 
00625                           radius.glbdst );
00626                         cdEXIT(EXIT_FAILURE);
00627                 }
00628                 factor = radius.glbden*pow(radius.glbrad/radius.glbdst,radius.glbpow);
00629                 fac2 = radius.glbden*pow(radius.glbrad/(radius.glbdst - (realnum)radius.drad),radius.glbpow);
00630                 if( fac2/factor > 1. + DNGLOB )
00631                 {
00632                         /* DNGLOB is relative change in density allowed this zone, 0.10 */
00633                         GlobDr = radius.drad*DNGLOB/(fac2/factor - 1.);
00634                 }
00635                 else
00636                 {
00637                         GlobDr = BigRadius;
00638                 }
00639                 /* GlobDr = GLBDST * DNGLOB * (GLBRAD/GLBDST)**(-GLBPOW) / GLBPOW */
00640                 GlobDr = MIN2(GlobDr,radius.glbdst/20.);
00641         }
00642         else
00643         {
00644                 GlobDr = BigRadius;
00645         }
00646 
00647         /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
00648         hdnew = 0.;
00649         if( strncmp( dense.chDenseLaw , "DLW" , 3) == 0 )
00650         {
00651                 /* one of the special density laws, first get density at possible next dr */
00652                 if( strcmp(dense.chDenseLaw,"DLW1") == 0 )
00653                 {
00654                         hdnew = dense_fabden(radius.Radius+radius.drad,radius.depth+
00655                           radius.drad);
00656                 }
00657                 else if( strcmp(dense.chDenseLaw,"DLW2") == 0 )
00658                 {
00659                         hdnew = dense_tabden(radius.Radius+radius.drad,radius.depth+
00660                           radius.drad);
00661                 }
00662                 else
00663                 {
00664                         fprintf( ioQQQ, " dlw insanity in radius_next\n" );
00665                         cdEXIT(EXIT_FAILURE);
00666                 }
00667                 drTab = fabs(hdnew-dense.gas_phase[ipHYDROGEN])/MAX2(hdnew,dense.gas_phase[ipHYDROGEN]);
00668                 drTab = radius.drad*MAX2(0.2,0.10/MAX2(0.01,drTab));
00669         }
00670         else
00671         {
00672                 drTab = BigRadius;
00673         }
00674 
00675         /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
00676         /* special density law */
00677         if( strcmp(dense.chDenseLaw,"DLW1") == 0 )
00678         {
00679                 dnew = fabs(dense_fabden(radius.Radius+radius.drad,radius.depth+
00680                   radius.drad)/dense.gas_phase[ipHYDROGEN]-1.);
00681                 /* DNGLOB is relative change in density allowed this zone, 0.10 */
00682                 if( dnew == 0. )
00683                 {
00684                         SpecDr = radius.drad*3.;
00685                 }
00686                 else if( dnew/DNGLOB > 1.0 )
00687                 {
00688                         SpecDr = radius.drad/(dnew/DNGLOB);
00689                 }
00690                 else
00691                 {
00692                         SpecDr = BigRadius;
00693                 }
00694         }
00695         else
00696         {
00697                 SpecDr = BigRadius;
00698         }
00699 
00700         /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
00701         /* check grain line heating dominates
00702          * this is important in PDR and HII region calculations
00703          * >>chng 97 jul 03, added following check */
00704         if( thermal.heating[0][13]/thermal.htot > 0.2 )
00705         {
00706                 /* >>chng 01 jan 03, following returns 0 when NO light at all,
00707                  * had failed with botched assert */
00708                 GrainRateDr(&grfreqm,&gropacm);
00709                 DrGrainHeat = 1.0/MAX2(1e-30,gropacm*geometry.FillFac*geometry.DirectionalCosin);
00710         }
00711         else
00712         {
00713                 gropacm = 0.;
00714                 grfreqm = 0.;
00715                 DrGrainHeat = BigRadius;
00716         }
00717 
00718         /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
00719         /* check if line heating dominates
00720          * this is important in high metallicity models */
00721         if( thermal.heating[0][22]/thermal.htot > 0.2 )
00722         {
00723                 FndLineHt(&level,&ipStrong,&Strong);
00724                 if( Strong/thermal.htot > 0.1 )
00725                 {
00726                         if( level == 1 )
00727                         {
00728                                 /* a level1 line was the heat source (usual case)
00729                                  * drLineHeat = MAX(1.0,TauLines(ipLnTauIn,ipStrong)*0.2) /
00730                                  *  1      TauLines(ipLnDTau,ipStrong) */
00731                                 TauInwd = TauLines[ipStrong].Emis->TauIn;
00732                                 TauDTau = TauLines[ipStrong].Emis->PopOpc * TauLines[ipStrong].Emis->opacity / 
00733                                         DoppVel.doppler[TauLines[ipStrong].Hi->nelem-1];
00734                         }
00735                         else if( level == 2 )
00736                         {
00737                                 /* a atom_level2 line was the heat source
00738                                  * (bad case since not as well treated)
00739                                  * drLineHeat = MAX(1.0,WindLine(ipLnTauIn,ipStrong)*0.2) /
00740                                  *  1      WindLine(ipLnDTau,ipStrong) */
00741                                 TauInwd = TauLine2[ipStrong].Emis->TauIn;
00742                                 TauDTau = TauLine2[ipStrong].Emis->PopOpc * TauLine2[ipStrong].Emis->opacity / 
00743                                         DoppVel.doppler[TauLine2[ipStrong].Hi->nelem-1];
00744                         }
00745                         else if( level == 3 )
00746                         {
00747                                 /* a 12CO line was the heat source
00748                                  * (bad case since not as well treated)
00749                                  * drLineHeat = MAX(1.0,WindLine(ipLnTauIn,ipStrong)*0.2) /
00750                                  *  1      WindLine(ipLnDTau,ipStrong) */
00751                                 TauInwd = C12O16Rotate[ipStrong].Emis->TauIn;
00752                                 TauDTau = C12O16Rotate[ipStrong].Emis->PopOpc * C12O16Rotate[ipStrong].Emis->opacity / 
00753                                         DoppVel.doppler[C12O16Rotate[ipStrong].Hi->nelem-1];
00754                         }
00755                         else if( level == 4 )
00756                         {
00757                                 /* a 13CO line was the heat source
00758                                  * (bad case since not as well treated)
00759                                  * drLineHeat = MAX(1.0,WindLine(ipLnTauIn,ipStrong)*0.2) /
00760                                  *  1      WindLine(ipLnDTau,ipStrong) */
00761                                 TauInwd = C13O16Rotate[ipStrong].Emis->TauIn;
00762                                 TauDTau = C13O16Rotate[ipStrong].Emis->PopOpc * C13O16Rotate[ipStrong].Emis->opacity / 
00763                                         DoppVel.doppler[C13O16Rotate[ipStrong].Hi->nelem-1];
00764                         }
00765                         else if( level == 5 )
00766                         {
00767                                 /* >>chng 03 dec 07, this branch had been left off, caught by Hiroaki Oyaizu */
00768                                 /* a hyperfine transition */
00769                                 TauInwd = HFLines[ipStrong].Emis->TauIn;
00770                                 TauDTau = HFLines[ipStrong].Emis->PopOpc * HFLines[ipStrong].Emis->opacity / 
00771                                         DoppVel.doppler[HFLines[ipStrong].Hi->nelem-1];
00772                         }
00773                         else
00774                         {
00775                                 /* this is insane, since Strong was set, but not level */
00776                                 fprintf( ioQQQ, " PROBLEM radius_next Strong line heating set, but not level.\n" );
00777                                 TotalInsanity();
00778                         }
00779                         /* in following logic cannot use usual inverse opacity,
00780                          * since line heating competes with escape probability,
00781                          * so is effective at surprising optical depths */
00782                         if( TauDTau > 0. )
00783                         {
00784                                 drLineHeat = MAX2(1.,TauInwd)*0.4/TauDTau;
00785                         }
00786                         else
00787                         {
00788                                 drLineHeat = BigRadius;
00789                         }
00790                 }
00791                 else
00792                 {
00793                         TauInwd = 0.;
00794                         drLineHeat = BigRadius;
00795                         ipStrong = 0;
00796                         Strong = 0.;
00797                 }
00798         }
00799         else
00800         {
00801                 TauInwd = 0.;
00802                 drLineHeat = BigRadius;
00803                 ipStrong = 0;
00804                 level = 0;
00805                 Strong = 0.;
00806         }
00807 
00808         /* >>chng 03 mar 03, add this logic */
00809         /* do not let change in cooling/heating due to H2 become too large */
00810         drH2_heat_cool = BigRadius;
00811         if( lgFirstCall )
00812         {
00813                 Old_H2_heat_cool = 0.;
00814         }
00815         else if( !thermal.lgTemperatureConstant )
00816         {
00817                 /* this is case where temperature has not been set */
00818                 /* compare total heating - cooling due to h2 with total due to everything */
00819                 double H2_heat_cool = (fabs(hmi.HeatH2Dexc_used)+fabs(hmi.HeatH2Dish_used)) / thermal.htot;
00820                 if( H2_heat_cool > 0.1 )
00821                 {
00822                         dH2_heat_cool = fabs( H2_heat_cool - Old_H2_heat_cool );
00823                         dH2_heat_cool = SDIV(dH2_heat_cool);
00824                         /* >>chng 05 oct 27, had been taking 20% of original radius - this caused zoning
00825                          * to become very fine and may not have been needed - change from 0.2 to 0.5 */
00826                         /*drH2_heat_cool = radius.drad*MAX2(0.2,0.05/dH2_heat_cool);*/
00827                         drH2_heat_cool = radius.drad*MAX2(0.5,0.05/dH2_heat_cool);
00828                 }
00829                 else
00830                 {
00831                         drH2_heat_cool = BigRadius;
00832                 }
00833         }
00834         Old_H2_heat_cool = (fabs(hmi.HeatH2Dexc_used)+fabs(hmi.HeatH2Dish_used)) / thermal.htot;
00835 
00836         /* >>chng 03 mar 04, add this logic */
00837         /* do not let change in H2 and heavy element molecular abundances become too large 
00838          * "change in heav ele mole abundance, d(mole)/elem" */
00839         drH2_abund = BigRadius;
00840         dr_mole_abund = BigRadius;
00841         mole_dr_change = -1;
00842         if( nzone>=4 )
00843         {
00844                 /* first do H2 abundance */
00845                 double Old_H2_abund;
00846                 /* >>chng 04 dec 15, do not use special logic when large h2 is turned on */
00847                 int i;
00848                 Old_H2_abund = struc.H2_molec[ipMH2g][nzone-3] + struc.H2_molec[ipMH2s][nzone-3];
00849                 /* >>chng 03 jun 16, limit from 0.01 to 0.001, some models fell over H2 front due to
00850                  * large zoning, when large H2 was just this caused oscillations in solomon process */
00851                 /* >>chng 03 nov 22, from > 0.001 to > 3e-4, models that start almost in H2 have
00852                  * rapid increase in H2 at shallow depths, start sensing this sooner */
00853                 /* >>chng 03 dec 10, from 3e-4 to 1e-4, to get smaller chagned in leiden1.in test */
00854                 /* radius_next keys from change in H2 abundance, d(H2)/H */
00855                 /* >>chng 04 mar 11, start sensing H2 earlier since otherwise step size
00856                  * needs to become way too small tooo quickly.  limit changed from 1e-4 to 1e-6 */
00857                 /* >>chng 04 jun 29, fromo > 1e-6 to >1e-8, some pdr's had too large a change in H2 */
00858 
00859                 if( 2.*hmi.H2_total/dense.gas_phase[ipHYDROGEN] > 1e-8 )
00860                 {
00861                         double fac = 0.2;
00862                         /* this is percentage change in H2 density - "change in H2 abundance" */
00863                         dH2_abund = 2.*fabs( hmi.H2_total - Old_H2_abund ) / hmi.H2_total;
00864                         /* in testing th85ism the dH2_abund did come out zero */
00865                         /* >>chng 03 jun 16, change d(H2) from 0.05 to 0.1, fine resolution of H2/H exposed
00866                          * small numerical oscillations in Solomon process */
00867                         /* >>chng 03 nov 22, smallest possible ratio of dr(next)/dr changed from
00868                          * 0.2 to 0.05, models that started almost in H2 front need to be able to sense it */
00869                         /*drH2_abund = radius.drad*MAX2(0.05,fac/SDIV(dH2_abund) );*/
00870                         /* >>chng 04 mar 15, with such small possible changes in dr, only 0.05, a thermal front
00871                          * can easily cause large changes in H2 and T that are not due to zoning, but to the
00872                          * discontinuity.  make smallest change larger to prevent hald due to dr */
00873                         /* >>chng 05 aug 23, thermal front allowed dr to become much too small
00874                          * chng from 0.02 to .6 */
00875                         dH2_abund = SDIV(dH2_abund);
00876                         drH2_abund = radius.drad*MAX2(0.2,fac/dH2_abund );
00877                 }
00878                 else
00879                 {
00880                         drH2_abund = BigRadius;
00881                 }
00882 
00883                 /* check how molecular fractions of all heavy elements are chaning relative 
00884                  * to total gas phase abundance */
00885                 dr_mole_abund = BigRadius;
00886                 /* >>chng 04 jun 02, upper limit had been all species but now limit to real
00887                  * molecules since do not want this logic to work with the ions */
00888                 for(i=0; i<mole.num_comole_calc; ++i)
00889                 {
00890                         realnum abund,
00891                                 abund_limit;
00892                         if( !dense.lgElmtOn[COmole[i]->nelem_hevmol] || COmole[i]->n_nuclei <= 1 )
00893                                 continue;
00894                         /* >>chng 03 sep 21, add CO logic */
00895                         /* >>chng 04 mar 30, generalize to any molecule at all */
00896                         /* >>chng 04 mar 31 lower limit to abund had been 0.01, change
00897                          * to 0.001 to pick up approach to molecular fronts */
00898                         abund = COmole[i]->hevmol/dense.gas_phase[COmole[i]->nelem_hevmol];
00899                         /* is this an ice?  need special logic for ice since density increases
00900                          * exponentially when grain temp falls below sublimation temperature 
00901                          * >>chng 05 dec 06 - detect changes for smaller abundances for ices
00902                          * due to large changes in ice abundances */
00903                         if( COmole[i]->lgGas_Phase )
00904                         {
00905                                 abund_limit = 1e-3f;
00906                         }
00907                         else
00908                         {
00909                                 /* this is an ice - track its abundance at lower values so that
00910                                  * we resolve the sublimation transition region */
00911                                 abund_limit = 1e-5f;
00912                         }
00913 
00914                         if( abund > abund_limit )
00915                         {
00916                                 double drmole_one, relative_change, relative_change_denom;
00917                                 /* >>chng 05 dec 08, use smaller abundance for the denominator since just taking
00918                                  * current abundance will overlook case where current density is vastly large
00919                                  * than old density */
00920                                 if( struc.CO_molec[i][nzone-3]>SMALLFLOAT )
00921                                 {
00922                                         relative_change_denom = MIN2( COmole[i]->hevmol , struc.CO_molec[i][nzone-3] );
00923                                 }
00924                                 else
00925                                 {
00926                                         relative_change_denom = COmole[i]->hevmol;
00927                                 }
00928                                 /* the relative change in the abundance */
00929                                 relative_change = 
00930                                         /*fabs( COmole[i]->hevmol - struc.CO_molec[i][nzone-3] ) / COmole[i]->hevmol;*/
00931                                         /* >>chng 05 dec 08, use smaller abundance */
00932                                         fabs( COmole[i]->hevmol - struc.CO_molec[i][nzone-3] ) / relative_change_denom;
00933                                 /* >>chng 03 jun 16, change from 0.05 to 0.1, fine resolution of H2/H exposed
00934                                  * small numerical oscillations in Solomon process */
00935                                 /* >>chng 04 jun 02, from 0.1 back to 0.05, more extensive CO etc network
00936                                  * caused oscillations in SiO abundance and Si Si+ density. */
00937                                 /* >>chng 04 aug 03, from 0.05 to 0.035, leiden pdr model v2 had major
00938                                  * jump in eden */
00939                                 /* >>chng 04 oct 18, from 0.035 back to 0.05, leiden pdr v2 actually due to having
00940                                  * PAHs in fully molecular limit (??), this caused cool flow pdr grid to trip on
00941                                  * too small dr */
00942                                 relative_change = SDIV(relative_change);
00943                                 /*>>chng 05 dec 08, relative_change must be less than one - with
00944                                  * revised logic above can be bigger than one */
00945                                 if( relative_change > 1. )
00946                                         relative_change = 1./relative_change;
00947                                 /*drmole_one = radius.drad*MAX2(0.2,0.035/relative_change );*/
00948                                 /* >>chng 05 aug 23, thermal front allowed dr to become much too small
00949                                  * chng from 0.02 to .6 */
00950                                 /*drmole_one = radius.drad*MAX2(0.2,0.05/relative_change );*/
00951                                 drmole_one = radius.drad*MAX2(0.6,0.05/relative_change );
00952                                 /* final dr will be the smallest we encounter */
00953                                 if( drmole_one < dr_mole_abund )
00954                                 {
00955                                         /* this is the dr used to set next dr - keep track of which moe was changing */
00956                                         dr_mole_abund = drmole_one;
00957                                         mole_dr_change = i;
00958                                         dCO_abund = relative_change;
00959                                 }
00960                         }
00961                 }
00962         }
00963 
00964         /* some consideration due to big H2 molecule */
00965         drSolomon_BigH2 = H2_DR();
00966 
00967         /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
00968         /* can't make drmax large deep in model since causes feedback
00969          * oscillations with changes in heating or destruction rates
00970          * >>chng 96 oct 15, change from 5 to 10 */
00971         if( nzone < 5 )
00972         {
00973                 /* >>chng 96 nov 11, had been 4 * drad up to 11, change to following
00974                  * to be similar to c90.01, max of 1.3 between 5 and 11 
00975                  * >>chng 04 oct 29  geometry.DirectionalCosin was ioncorrect applied to this factor  */
00976                 drmax = 4.*radius.drad;
00977         }
00978         else
00979         {
00980                 drmax = 1.3*radius.drad;
00981         }
00982 
00983         /* >>chng 05 apr 05, do not sense temp oscillation, so that we can move past this
00984          * point if it occurs */
00985 #       if 0
00986         /* look for oscillations in electron density or tempeature - freeze dr if these occur */
00987         dr_ne_oscil = BigRadius;
00988         dr_te_oscil = BigRadius;
00989         if( nzone >= 11 )
00990         {
00991                 /* >>chng 96 oct 15, do not let zones increase if oscillations present */
00992                 /* >>chng 96 oct 31, error to declare oscillation propto toler, the 
00993                  *heating cooling tolerance */
00994                 realnum errorHC = POW2(conv.HeatCoolRelErrorAllowed);
00995                 realnum errorNe = (realnum)POW2(conv.EdenErrorAllowed );
00996                 limit = nzone -2;
00997                 ASSERT( limit < struc.nzlim );
00998                 for( k=nzone - 10; k < limit; k++ )
00999                 {
01000                         /* check that square of change both chng sign and is
01001                          * greater than square of heat-tool error */
01002                         if( (struc.testr[k-1] - struc.testr[k])/struc.testr[k]*
01003                           (struc.testr[k] - struc.testr[k+1])/struc.testr[k] < 
01004                           -errorHC )
01005                         {
01006                                 dr_te_oscil = radius.drad;
01007                         }
01008                         /* small residiual is to allow 0.01 rel error */
01009                         if( (struc.ednstr[k-1] - struc.ednstr[k])/struc.ednstr[k]*
01010                           (struc.ednstr[k] - struc.ednstr[k+1])/struc.ednstr[k] < 
01011                           -errorNe )
01012                         {
01013                                 /* >>chng 96 oct 15, do not let zones increase if oscillations present */
01014                                 /* radius_next keys from electron density oscillation*/
01015                                 dr_ne_oscil = radius.drad;
01016                         }
01017                 }
01018         }
01019         /* >>chng 05 apr 05, do not sense temp oscillation, so that we can move past this
01020          * point if it occurs */
01021         dr_ne_oscil = BigRadius;
01022         dr_te_oscil = BigRadius;
01023 #       endif
01024 
01025         /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
01026         /* check on several convergence criteria */
01027 
01028         /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
01029         if( !conv.lgConvTemp )
01030         {
01031                 drFail = radius.drad/2.;
01032         }
01033         else
01034         {
01035                 drFail = BigRadius;
01036         }
01037 
01038         /* time dependent recombination - conditions become very homogeneous and 
01039          * crash into I fronts - use last iteration's  radius as guess of current case */
01040         recom_dr_last_iter = BigRadius;
01041         if( dynamics.lgStatic && dynamics.lgRecom )
01042         {
01043                 /* first time through nzone == 1 */
01044                 static long int nzone_recom = -1;
01045                 if( nzone<=1 )
01046                 {
01047                         /* initialization case */
01048                         nzone_recom = 0;
01049                 }
01050                 else if( radius.depth < struc.depth_last[struc.nzone_last-1] && 
01051                         nzone_recom < struc.nzone_last )
01052                 {
01053                         ASSERT( nzone_recom>=0 && nzone_recom<struc.nzone_last );
01054                         /* case where we are within previous computed structure 
01055                          * first possibly increase nzone_recom so that it points deeper
01056                          * than this zone */
01057                         while( struc.depth_last[nzone_recom] < radius.depth &&
01058                                 nzone_recom < struc.nzone_last-1 )
01059                                 ++nzone_recom;
01060                         recom_dr_last_iter = struc.drad_last[nzone_recom]*3.;
01061                 }
01062                 else
01063                 {
01064                         /* case where we have overrun the previous iteration's saved structure */
01065                         recom_dr_last_iter = BigRadius;
01066                 }
01067         }
01068 
01069         /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
01070         /* change in electron density - radius_next keys from change in elec den,
01071          * remember old electron density during first call */
01072         /* this is low ionization solution */
01073         if( nzone > 2 )
01074         {
01075                 /* next is-2 since nzone is on physics not c scale, and we want previous zone */
01076                 Efrac_old = struc.ednstr[nzone-3]/struc.gas_phase[ipHYDROGEN][nzone-3];
01077                 Efrac_new = dense.eden/dense.gas_phase[ipHYDROGEN];
01078                 dEfrac = fabs(Efrac_old-Efrac_new) / Efrac_new;
01079 
01080                 if( dEfrac > SMALLFLOAT )
01081                 {
01082                         double fac = 0.04;
01083                         /* >>chng 03 dec 25, use smaller rel change in elec frac when most elec in ipMole or grains */
01084                         /* >>chng 04 sep 14, change to from from metals but comment out */
01085                         /* >>chng 04 sep 17, change to from from metals - uncomment */
01086                         if( dense.eden_from_metals > 0.5 )
01087                         {
01088                                 /* >>chng 04 sep 18, change 0.02 from 0.01 */
01089                                 /* >>chng 04 sep 18, change 0.02 from 0.04 */
01090                                 fac = 0.04;
01091                         }
01092                         /* >>chng 04 feb 23, add test on hydrogen being predomintly
01093                          * recombined due to three-body recom, which is very sensitive
01094                          * to the electron density - but only do this in partially ionized medium */
01095                         else if( iso.RecomCollisFrac[ipH_LIKE][ipHYDROGEN] > 0.8 && 
01096                                 dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN]>0.1 &&
01097                                 dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN]<0.8 )
01098 
01099                         {
01100                                 fac = 0.02;
01101                         }
01102                         /* >>chng 04 sep 17, change to 0.1 from 0.2 */
01103                         drEfrac = radius.drad*MAX2(0.1,fac/dEfrac);
01104                 }
01105                 else
01106                 {
01107                         drEfrac = BigRadius;
01108                 }
01109         }
01110         else
01111         {
01112                 dEfrac = 0.;
01113                 drEfrac = BigRadius;
01114                 Efrac_old = 0.;
01115                 Efrac_new = 0.;
01116         }
01117 
01118         /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
01119         /* do not let thickness get too large
01120          *  1    '' radius_next keys from relative depth'',e11.3)') */
01121         if( nzone > 20 )
01122         {
01123                 /*drDepth = radius.depth/20.;*/
01124                 /* >>chng 02 nov 05, change from 1/20 to 1/10 wasted zones early on */
01125                 drDepth = radius.depth/10.;
01126         }
01127         else
01128         {
01129                 drDepth = BigRadius;
01130         }
01131 
01132         /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
01133         /* case where stopping thickness or edge specified, need to approach slowly */
01134         thickness_total = BigRadius;
01135         DepthToGo = BigRadius;
01136         if( StopCalc.HColStop < 5e29 )
01137         {
01138                 double coleff = SDIV( dense.gas_phase[ipHYDROGEN]*geometry.FillFac);
01139                 DepthToGo = MIN2(DepthToGo ,
01140                         (StopCalc.HColStop-colden.colden[ipCOL_HTOT]) / coleff );
01141                 /* >>chng 03 oct 28, forgot to div col den above by eff density */
01142                 thickness_total = MIN2(thickness_total , StopCalc.HColStop / coleff );
01143         }
01144 
01145         if( StopCalc.colpls < 5e29 )
01146         {
01147                 double coleff = (double)SDIV( (dense.xIonDense[ipHYDROGEN][1])*geometry.FillFac);
01148                 DepthToGo = MIN2(DepthToGo ,
01149                         (StopCalc.colpls-colden.colden[ipCOL_Hp]) / coleff );
01150                 thickness_total = MIN2(thickness_total , StopCalc.colpls / coleff );
01151         }
01152 
01153         if( StopCalc.col_h2 < 5e29 )
01154         {
01155                 /* >>chng 03 apr 15, add this molecular hydrogen */
01156                 double coleff = (double)SDIV( hmi.H2_total*geometry.FillFac);
01157                 DepthToGo = MIN2(DepthToGo ,
01158                         (StopCalc.col_h2-colden.colden[ipCOL_H2g]-colden.colden[ipCOL_H2s]) / coleff );
01159                 thickness_total = MIN2(thickness_total , StopCalc.col_h2 / coleff );
01160         }
01161 
01162         if( StopCalc.col_h2_nut < 5e29 )
01163         {
01164                 /* >>chng 03 apr 15, add this molecular hydrogen */
01165                 double coleff = (double)SDIV( (2*hmi.H2_total+dense.xIonDense[ipHYDROGEN][1])*geometry.FillFac);
01166                 DepthToGo = MIN2(DepthToGo ,
01167                         (StopCalc.col_h2_nut-(2*(colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s])+dense.xIonDense[ipHYDROGEN][1])) / coleff  );
01168                 thickness_total = MIN2(thickness_total , StopCalc.col_h2_nut / coleff );
01169         }
01170 
01171         if( StopCalc.col_H0_ov_Tspin < 5e29 )
01172         {
01173                 /* >>chng 05 jan 09, add n(H0)/Tspin */
01174                 double coleff = (double)SDIV( dense.xIonDense[ipHYDROGEN][0] / hyperfine.Tspin21cm*geometry.FillFac );
01175                 DepthToGo = MIN2(DepthToGo ,
01176                         (StopCalc.col_H0_ov_Tspin - colden.H0_ov_Tspin) / coleff  );
01177                 thickness_total = MIN2(thickness_total , StopCalc.col_H0_ov_Tspin / coleff );
01178         }
01179 
01180         if( StopCalc.col_monoxco < 5e29 )
01181         {
01182                 /* >>chng 03 apr 15, add this, CO */
01183                 double coleff = (double)SDIV( (findspecies("CO")->hevmol)*geometry.FillFac);
01184                 DepthToGo = MIN2(DepthToGo ,
01185                         (StopCalc.col_monoxco-findspecies("CO")->hevcol) / coleff );
01186                 thickness_total = MIN2(thickness_total , StopCalc.col_monoxco / coleff );
01187         }
01188 
01189         if( StopCalc.colnut < 5e29 )
01190         {
01191                 double coleff = (double)SDIV( (dense.xIonDense[ipHYDROGEN][0])*geometry.FillFac);
01192                 DepthToGo = MIN2(DepthToGo ,
01193                         (StopCalc.colnut - colden.colden[ipCOL_H0]) / coleff );
01194                 thickness_total = MIN2(thickness_total , StopCalc.colnut / coleff );
01195         }
01196 
01197         /* this is case where outer radius is set */
01198         if( radius.router[iteration-1] < 5e29 )
01199         {
01200                 thickness_total = MIN2(thickness_total , radius.router[iteration-1] );
01201                 DepthToGo = MIN2(DepthToGo ,
01202                         radius.router[iteration-1] - radius.depth );
01203         }
01204 
01205         /* this is case where stopping optical depth was specified */
01206         if( StopCalc.iptnu != rfield.nupper )
01207         {
01208                 /* end optical depth has been specified */
01209                 double dt = SDIV(opac.opacity_abs[StopCalc.iptnu-1]*geometry.FillFac);
01210                 DepthToGo = MIN2(DepthToGo ,
01211                         (StopCalc.tauend-opac.TauAbsGeo[0][StopCalc.iptnu-1])/dt);
01212         }
01213         /* stop AV - usually this is dust, but we consider all opacity sources,
01214          * so always include this */
01215         /* compute some average grain properties */
01216         AV_to_go = BigRadius;
01217         if( rfield.opac_mag_V_extended > SMALLFLOAT && rfield.opac_mag_V_point>SMALLFLOAT )
01218         {
01219                 /* by default stop av is very large, and opacity can be very small, so ratio
01220                  * goes to inf - work with logs to see how big the number is */
01221                 double ave = log10(StopCalc.AV_extended - rfield.extin_mag_V_extended) - 
01222                         log10(rfield.opac_mag_V_extended);
01223                 double avp = log10(StopCalc.AV_point - rfield.extin_mag_V_point) - 
01224                         log10(rfield.opac_mag_V_point);
01225                 AV_to_go = MIN2( ave , avp );
01226                 if( AV_to_go < 37. )
01227                 {
01228                         AV_to_go = pow(10., AV_to_go );
01229                         /* this is to make sure that we go slightly deeper than AV so that
01230                          * we ttrigger this stop */
01231                         AV_to_go *= 1.0001;
01232                 }
01233                 else
01234                         AV_to_go = BigRadius;
01235                 /*fprintf(ioQQQ,"DEBUG next dr %.3e %.3e %.3e\n", AV_to_go , ave, avp );*/
01236         }
01237 
01238         /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
01239         /* set dr if one of above tests have triggered */
01240         if( DepthToGo <= 0. )
01241         {
01242                 TotalInsanity();
01243         }
01244         else if( DepthToGo < BigRadius )
01245         {
01246                 /* want to approach the outer edge slowly,
01247                  * the need for this logic is most evident in brems.in - 
01248                  * HI fraction varies across coronal model */
01249                 drThickness = MIN2( thickness_total/10. , DepthToGo );
01250         }
01251         else
01252         {
01253                 drThickness = BigRadius;
01254         }
01255 
01256         /*fprintf(ioQQQ,
01257                 "DEBUG depth2go z%li drThickness2 = %e %e %e\n",
01258                 nzone , drThickness , thickness_total , DepthToGo );*/
01259 
01260         /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
01261         /* spherical models, do not want delta R/R big */
01262         drSphere = radius.Radius*0.04;
01263 
01264         /* optical depth to electron scattering */
01265         /* >>chng 04 jun 16, add filling factor, was missing */
01266         dRTaue = radius.drChange/(dense.eden*6.65e-25*geometry.FillFac);
01267         /* >>chng 02 oct 06, increase dr when constant temperature,
01268          * to prevent some ct models from taking too many cells */
01269         if( thermal.lgTemperatureConstant ) 
01270                 dRTaue *= 3.;
01271 
01272         /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
01273         if( dense.flong != 0. )
01274         {
01275                 drFluc = 0.628/2./dense.flong;
01276                 /* >>chng 04 sep 18, caused cautions that ionization jumps occurred.
01277                  * set to have the value */
01278                 drFluc /= 2.;
01279         }
01280         else
01281         {
01282                 drFluc = BigRadius;
01283         }
01284 
01285         /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
01286         /* if density fluctuations in place then override change in heat
01287          * for dr set */
01288         if( strcmp(dense.chDenseLaw,"SINE") == 0 && dense.lgDenFlucOn )
01289         {
01290                 drdHdStep = BigRadius;
01291         }
01292 
01293         /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
01294         /*active dr sets */
01295         /* we are deep into model, use logic since we have several zones
01296          * of old data */
01297         radius.drNext = MIN2( MIN4( drmax, drInter, drLineHeat, winddr ),
01298                               MIN4( drFluc, GlobDr, DrGrainHeat, dr_change_heavy ) );
01299         radius.drNext = MIN2( MIN4( radius.drNext, SpecDr, drFail, WindAccelDR ),
01300                               MIN3( drSphere, radius.sdrmax, dRTaue ) );
01301         radius.drNext = MIN3( MIN3( radius.drNext, drDest, drdTdStep ),
01302                               MIN3( drdHdStep, drConPres, drTab ),
01303                               MIN3( drSolomon_BigH2, drLeiden_hack, recom_dr_last_iter ) );
01304         radius.drNext = MIN3( MIN4( radius.drNext, drHion, drDepth, dr_mole_abund ),
01305                               MIN4( AV_to_go, drEfrac, drHMase, drThickness ),
01306                               MIN3( drHe1ovHe2, drH2_heat_cool, drH2_abund ) );
01307 
01308         /*fprintf(ioQQQ,
01309                 "DEBUG depth2go drNext = %e \n",                radius.drNext );*/
01310 
01311         /* keep dr constant in first two zones, if it wants to increase,
01312          * but always allow it to decrease.
01313          * to guard against large increases in efrac for balmer cont photo dominated models,
01314          */
01315         if( nzone <= 1 && radius.drNext > OldDR )
01316         {
01317                 radius.drNext = OldDR;
01318         }
01319 
01320         /* option to force min drad */
01321         if( radius.drNext < radius.sdrmin )
01322         {
01323                 radius.drNext = radius.sdrmin;
01324         }
01325 
01328         /* a pressure failure has occurred - keep zone the same time, hoping to pass through
01329          * troubled region */
01330         if( !conv.lgConvPres || !conv.lgConvTemp )
01331         {
01332                 radius.drNext = radius.drad;
01333         }
01334 
01335         /* >>chng 05 aug 05, in case of thermal front, where temp is falling quickly and
01336          * conditions change very fast, the zone thickness does not really affect the change
01337          * in conditions and can cause zoning to become very very thin, which causes
01338          * an abort.  this occurs between 200 and 1000K.  if we are doing temp soln,
01339          * temp is between these values, and temp is changing rapidly, do not make sone
01340          * thickness much smaller */
01341         drThermalFront = BigRadius;
01342         /* do not use thermal front logic in case of recombination time dep cloud */
01343         if( nzone >=5 && !dynamics.lgStatic )
01344         {
01345                 /* >>chng 05 aug 23, upper bound of thermal from from 1000K to 4000K */
01346                 /*if( phycon.te > 200. && phycon.te < 1000. && */
01347                 if( phycon.te > 200. && phycon.te < 3000. && 
01348                         /* >>chng 05 aug 23, from > 10% in zone to to two zones > 5%,
01349                          * to fix leiden v3 with large H2 */
01350                         (struc.testr[nzone-3] - phycon.te)/phycon.te > 0.02 &&
01351                         (struc.testr[nzone-4] - phycon.te)/phycon.te > 0.02 &&
01352                         (struc.testr[nzone-5] - phycon.te)/phycon.te > 0.02 )
01353                 {
01354                         /* the 0.91 is to make dr unique, so that print statement that
01355                          * follows will identify this reason */
01356                         drThermalFront = radius.drad * 0.91;
01357                         radius.drNext = drThermalFront;
01358                 }
01359         }
01360 
01361         /* dr = zero is a logical mistake */
01362         if( radius.drNext <= 0. )
01363         {
01364                 fprintf( ioQQQ, " radius_next finds insane drNext:%10.2e\n", 
01365                   radius.drNext );
01366                 fprintf( ioQQQ, " all drs follow:%10.2e%10.2e%10.2e%10.2e\n all drs follow:%10.2e%10.2e%10.2e%10.2e%10.2e\n all drs follow:%10.2e%10.2e%10.2e%10.2e\n", 
01367                   drmax, drInter, drLineHeat, winddr, drFluc, GlobDr, SpecDr, 
01368                   drFail, drSphere, radius.sdrmax, dRTaue, 
01369                   OldH2Abund, drDepth );
01370                 cdEXIT(EXIT_FAILURE);
01371         }
01372 
01373         /* set flag if dr set by maser */
01374         if( fp_equal( radius.drNext, drHMase ) )
01375         {
01376                 rt.lgMaserSetDR = true;
01377         }
01378 
01379         /* all this is to only punch on last iteration
01380          * the punch dr command is not really a punch command, making this necessary
01381          * lgDRon is set true if "punch dr" entered */
01382         if( punch.lgDROn )
01383         {
01384                 if( punch.lgDRPLst )
01385                 {
01386                         /* lgDRPLst was set true if "punch" had "last" on it */
01387                         if( iterations.lgLastIt )
01388                         {
01389                                 lgDoPun = true;
01390                         }
01391                         else
01392                         {
01393                                 lgDoPun = false;
01394                         }
01395                 }
01396                 else
01397                 {
01398                         lgDoPun = true;
01399                 }
01400         }
01401         else
01402         {
01403                 lgDoPun = false;
01404         }
01405         if( (trace.lgTrace && trace.lgDrBug) || lgDoPun )
01406         {
01407                 if( !conv.lgConvTemp && nzone > 0 )
01408                 {
01409                         fprintf( punch.ipDRout, "#>>>> A temperature failure occured.\n" );
01410                 }
01411                 if( !conv.lgConvPres && nzone > 0 )
01412                 {
01413                         fprintf( punch.ipDRout, "#>>>> A pressure failure occured.\n" );
01414                 }
01415 
01416                 /* this is common part of each line, the zone count, depth, chosen dr, and depth2go */
01417                 fprintf( punch.ipDRout , "%ld\t%.5e\t%.3e\t%.3e\t", nzone, radius.depth, radius.drNext, radius.Depth2Go );
01418 
01419                 /*=======begin active dr sets */
01420                 if( fp_equal( radius.drNext, drLineHeat ) )
01421                 {
01422                         if( level == 1 )
01423                         {
01424                                 strcpy( chLbl, chLineLbl(&TauLines[ipStrong]) );
01425                                 fprintf( punch.ipDRout, "level 1 line heating,%10.10s TauIn%10.2e%10.2e%10.2e\n", 
01426                                   chLbl, TauInwd, TauLines[ipStrong].Emis->pump, 
01427                                   TauLines[ipStrong].Emis->Pesc );
01428                         }
01429                         else if( level == 2 )
01430                         {
01431                                 strcpy( chLbl, chLineLbl(&TauLine2[ipStrong]));
01432                                 fprintf( punch.ipDRout, "level 2 line heating,%10.10s TauIn%10.2e%10.2e%10.2e\n", 
01433                                   chLbl, TauInwd, TauLine2[ipStrong].Emis->pump, 
01434                                   TauLine2[ipStrong].Emis->Pesc );
01435                         }
01436                         else
01437                         {
01438                                 fprintf( ioQQQ, " insanity pr line heat\n" );
01439                                 cdEXIT(EXIT_FAILURE);
01440                         }
01441                 }
01442 
01443                 else if( fp_equal( radius.drNext, drDepth ) )
01444                 {
01445                         fprintf( punch.ipDRout, "relative depth\n");
01446                 }
01447 
01448                 else if( fp_equal( radius.drNext, drThermalFront ) )
01449                 {
01450                         fprintf( punch.ipDRout, "thermal front logic\n");
01451                 }
01452 
01453                 else if( fp_equal( radius.drNext, dr_change_heavy ) )
01454                 {
01455                         ASSERT( ichange_heavy_nelem < LIMELM+3 );
01456                         fprintf( punch.ipDRout, 
01457                                 "change in ionization, element %s ion (C scale) %li rel change %.2e ion frac %.2e\n",
01458                                 elementnames.chElementName[ichange_heavy_nelem],
01459                                 ichange_heavy_ion , 
01460                                 change_heavy_frac_big ,
01461                                 frac_change_heavy_frac_big);
01462                 }
01463 
01464                 else if( fp_equal( radius.drNext, drThickness ) )
01465                 {
01466                         fprintf( punch.ipDRout, "depth to go\n");
01467                 }
01468 
01469                 else if( fp_equal( radius.drNext, AV_to_go ) )
01470                 {
01471                         fprintf( punch.ipDRout, "A_V to go\n");
01472                 }
01473 
01474                 else if( fp_equal( radius.drNext, drTab ) )
01475                 {
01476                         fprintf( punch.ipDRout, "spec den law, new old den%10.2e%10.2e\n", 
01477                           hdnew, dense.gas_phase[ipHYDROGEN] );
01478                 }
01479 
01480                 else if( fp_equal( radius.drNext, drHMase ) )
01481                 {
01482                         fprintf( punch.ipDRout, "H maser dTauMase=%10.2e %li %li %li %li\n", 
01483                                 rt.dTauMase,
01484                                 rt.mas_species,
01485                                 rt.mas_ion,
01486                                 rt.mas_hi,
01487                                 rt.mas_lo );
01488                 }
01489 
01490                 else if( fp_equal( radius.drNext, drHe1ovHe2 ) )
01491                 {
01492                         /* radius_next keys from change in He2/He1 ionization, old=%11.3e sv=%11.3e FrLya*/
01493                         fprintf( punch.ipDRout, "change in He0/He+ ionization, ratio %.2e\n", 
01494                           Old_He_atom_ov_ion );
01495                 }
01496 
01497                 else if( fp_equal( radius.drNext, drH2_heat_cool ) )
01498                 {
01499                         /* radius_next keys from change in H2 heating, old=%11.3e sv=%11.3e FrLya*/
01500                         fprintf( punch.ipDRout, "change in H2 heating/cooling, d(c,h)/H %.2e\n", 
01501                           dH2_heat_cool );
01502                 }
01503 
01504                 else if( fp_equal( radius.drNext, drH2_abund ) )
01505                 {
01506                         /* radius_next keys from change in H2 abundance, old=%11.3e sv=%11.3e FrLya*/
01507                         fprintf( punch.ipDRout, "change in H2 abundance, d(H2)/H %.2e\n", 
01508                           dH2_abund );
01509                 }
01510 
01511                 else if( fp_equal( radius.drNext, dr_mole_abund ) )
01512                 {
01513                         /* radius_next keys from change in CO mole abundance */
01514                         fprintf( punch.ipDRout, "change in heav ele mole abundance, d(mole)/elem %.2e mole=%i=%s\n", 
01515                           dCO_abund , mole_dr_change , COmole[mole_dr_change]->label);
01516                 }
01517 
01518                 else if( fp_equal( radius.drNext, drSolomon_BigH2 ) )
01519                 {
01520                         /* radius_next keys from change in H2 abundance, old=%11.3e sv=%11.3e FrLya*/
01521                         fprintf( punch.ipDRout, "change in big H2 Solomon rate line opt depth\n");
01522                 }
01523 
01524                 else if( fp_equal( radius.drNext, recom_dr_last_iter ) )
01525                 {
01526                         /* radius_next keys from previous iteration recom logic */
01527                         fprintf( punch.ipDRout, "previous iteration recom logic\n");
01528                 }
01529 
01530                 else if( fp_equal( radius.drNext, drHion ) )
01531                 {
01532                         fprintf( punch.ipDRout, "change in H ionization fm to%11.3e%11.3e\n", 
01533                           SaveOHIIoHI, OHIIoHI );
01534                 }
01535 
01536                 else if( fp_equal( radius.drNext, drEfrac ) )
01537                 {
01538                         fprintf( punch.ipDRout, 
01539                                 "change in elec den, rel chng:%11.3e, cur %11.3e old=%11.3e\n", 
01540                           dEfrac , Efrac_old , Efrac_new );
01541                 }
01542 
01543                 else if( fp_equal( radius.drNext, drdHdStep ) )
01544                 {
01545                         fprintf( punch.ipDRout, 
01546                                 "change in heating; current %10.3e delta=%10.3e\n", 
01547                           thermal.htot, dHdStep );
01548                 }
01549 
01550                 else if( fp_equal( radius.drNext, drLeiden_hack ) )
01551                 {
01552                         fprintf( punch.ipDRout, 
01553                                 "Leiden hack\n" );
01554                 }
01555 
01556                 else if( fp_equal( radius.drNext, drConPres ) )
01557                 {
01558                         fprintf( punch.ipDRout, "change in con accel\n"  );
01559                 }
01560 
01561                 else if( fp_equal( radius.drNext, drdTdStep ) )
01562                 {
01563                         fprintf( punch.ipDRout, 
01564                                 "change in temperature; current= %10.3e, dT/T= %.3f\n", 
01565                           phycon.te, dTdStep );
01566                 }
01567 
01568                 else if( fp_equal( radius.drNext, radius.sdrmin ) )
01569                 {
01570                         fprintf( punch.ipDRout, "sdrmin\n"  );
01571                 }
01572 
01573                 else if( fp_equal( radius.drNext, radius.sdrmax ) )
01574                 {
01575                         fprintf( punch.ipDRout, "sdrmax\n" );
01576                 }
01577 
01578                 else if( fp_equal( radius.drNext, drSphere ) )
01579                 {
01580                         fprintf( punch.ipDRout, "sphericity\n" );
01581                 }
01582 
01583                 else if( fp_equal( radius.drNext, dRTaue ) )
01584                 {
01585                         fprintf( punch.ipDRout, 
01586                                 "optical depth to electron scattering\n" );
01587                 }
01588 
01589                 else if( fp_equal( radius.drNext, drFail ) )
01590                 {
01591                         fprintf( punch.ipDRout, 
01592                                 "temperature failure.\n" );
01593                 }
01594 
01595                 else if( fp_equal( radius.drNext, drmax ) )
01596                 {
01597                         fprintf( punch.ipDRout, 
01598                                 "DRMAX; nu opc dr=%10.2e%10.2e%10.2e\n", 
01599                           freqm, opacm, radius.drChange/
01600                           SDIV(opacm) );
01601                 }
01602 
01603                 else if( fp_equal( radius.drNext, drInter ) )
01604                 {
01605                         fprintf( punch.ipDRout, 
01606                                 "cont inter nu=%10.2e opac=%10.2e\n", 
01607                           freqm, opacm );
01608                 }
01609 
01610                 else if( fp_equal( radius.drNext, DrGrainHeat ) )
01611                 {
01612 
01613                         fprintf( punch.ipDRout, 
01614                                 "grain heating nu=%10.2e opac=%10.2e\n", 
01615                           grfreqm, gropacm );
01616                 }
01617 
01618                 else if( fp_equal( radius.drNext, winddr ) )
01619                 {
01620                         fprintf( punch.ipDRout, 
01621                                 "Wind, dVelRelative=%.3e windv=%.3e\n", 
01622                            dVelRelative ,
01623                            wind.windv);
01624                 }
01625 
01626                 else if( fp_equal( radius.drNext, drFluc ) )
01627                 {
01628                         fprintf( punch.ipDRout, 
01629                                 "density fluctuations\n"  );
01630                 }
01631 
01632                 else if( fp_equal( radius.drNext, GlobDr ) )
01633                 {
01634                         fprintf( punch.ipDRout, 
01635                                 "GLOB law new dr=%10.2e HDEN=%10.2e\n", 
01636                                 GlobDr,
01637                           dense.gas_phase[ipHYDROGEN] );
01638                 }
01639 
01640                 else if( fp_equal( radius.drNext, SpecDr ) )
01641                 {
01642                         fprintf( punch.ipDRout, 
01643                                 "special law new dr=%10.2e HDEN=%10.2e\n", 
01644                                 SpecDr,
01645                           dense.gas_phase[ipHYDROGEN] );
01646                 }
01647 
01648                 else if( fp_equal( radius.drNext, OldDR ) )
01649                 {
01650                         fprintf( punch.ipDRout, "old DR.\n" );
01651                 }
01652 
01653                 else
01654                 {
01655                         fprintf( punch.ipDRout, 
01656                                 " %4ld radius_next keys from insanity %10.2e\n", 
01657                           nzone, radius.drNext );
01658 
01659                         fprintf( ioQQQ, 
01660                                 " %4ld radius_next keys from insanity %10.2e\n", 
01661                           nzone, radius.drNext );
01662                         cdEXIT(EXIT_FAILURE);
01663                 }
01664 
01665                 /*======end active dr sets */
01666         }
01667 
01668         /* this is general code that prevents zone thickness drNext from
01669          * becoming too thin, something that can happen for various bad reasons
01670          * HOWEVER we do not want to do this test for certain density laws,
01671          * for which very small zone thicknesses are unavoidable
01672          * the special cases are:
01673          * special density law,
01674          * globule density law,
01675          * carbon +-0 i front
01676          * the fluctuations command
01677          * drMinimum was set in radius_first to either sdrmin (set drmin) or
01678          * some fraction of the initial radius - it is always set
01679          * to something non-trivial.  
01680          * sdrmin is set with the "set dr" command - is SMALLFLOAT by default */
01681         if( ((strcmp(dense.chDenseLaw,"DLW1") != 0 && 
01682                 strcmp(dense.chDenseLaw ,"GLOB") != 0) )&& 
01683                 /* >>chng 04 feb 19, do not use this test - errors can still happen
01684                  * when all C is atomic! */
01685                 /*dense.xIonDense[ipCARBON][0]/dense.gas_phase[ipCARBON] < 0.05) && */
01686                 (dense.flong == 0.) &&
01687                 /* >>chng 01 aug 11, add check against stopping on depth to go */
01688                 radius.drNext != DepthToGo )
01689         {
01690                 /* don't let dr get smaller than drMinimum, if this resets drNext
01691                  * then code stops with warning that zones got too thin
01692                  * this can happen due to numerical oscillations, which the code
01693                  * tries to damp out by making the zones thinner.
01694                  * scale by radius.Radius/radius.rinner to stop very spherical 
01695                  * simulations from false trigger on dr too small */
01696                 if( radius.drNext* radius.Radius/radius.rinner  * 
01697                         /* drMinimum is drad * hden, to make it proportional to optical depth.
01698                          * This avoids false trigger across thermal fronts. */
01699                         dense.gas_phase[ipHYDROGEN] < 
01700                         radius.drMinimum )
01701                 {
01702                         radius.drNext = radius.drMinimum/ dense.gas_phase[ipHYDROGEN];
01703                         /* leaving this at true will cause the model to stop with a warning
01704                          * that the zone thickness is too small */
01705                         radius.lgDrMinUsed = true;
01706                         /* set abort handler */
01707                         lgAbort = true;
01708                         /* must decrement nzone, since we will not complete this zone, and will not have
01709                          * valid structure data for it */
01710                         --nzone;
01711                         fprintf( ioQQQ, 
01712                                 "\n DISASTER PROBLEM radius_next finds dr too small and aborts.  "
01713                                 "This is zone %ld iteration %ld.  The thickness was %.2e\n", 
01714                                 nzone, 
01715                                 iteration,
01716                                 radius.drNext);
01717                         fprintf( ioQQQ, 
01718                                 " If this simulation seems reasonable you can override this limit "
01719                                 "with the command SET DRMIN %.2e\n\n", 
01720                                 radius.drNext*2);
01721                         /* set abort flag */
01722                         lgAbort = true;
01723                         return(1);
01724                 }
01725         }
01726 
01727         /* factor to allow for slop in floating numbers */
01728 #       define  Z       1.0001
01729 
01730         /* following is to make thickness of model exact
01731          * n.b., on last zone, drNext can be NEGATIVE!!
01732          * DEPTH was incremented at start of zone calc in newrad,
01733          * has been outer edge of zone all throughout */
01734         radius.drNext = MIN2(radius.drNext,(radius.router[iteration-1]-
01735           radius.depth)*Z);
01736 
01737         /* this means outer limit exceeded */
01738         if( radius.drNext < 0. )
01739         {
01740                 radius.lgDrNeg = true;
01741         }
01742         else
01743         {
01744                 radius.lgDrNeg = false;
01745         }
01746 
01747         ASSERT( radius.drNext > 0. );
01748 
01749         if( trace.lgTrace )
01750         {
01751                 fprintf( ioQQQ, " radius_next chooses next drad drNext=%.4e; this drad was%12.4e\n", 
01752                   radius.drNext, radius.drad );
01753         }
01754 
01755         return( 0 );
01756 }
01757 
01758 /*ContRate called by radius_next to find energy of maximum continuum-gas interaction */
01759 STATIC void ContRate(double *freqm, 
01760   double *opacm)
01761 {
01762         long int i, 
01763           ipHi,
01764           iplow, 
01765           limit;
01766         double FreqH, 
01767           Freq_nonIonizing, 
01768           Opac_Hion, 
01769           Opac_nonIonizing, 
01770           Rate_max_Hion, 
01771           Rate_max_nonIonizing;
01772 
01773         DEBUG_ENTRY( "ContRate()" );
01774 
01775         /* 
01776          * find maximum continuum interaction rate,
01777          * these should be reset in following logic without exception,
01778          * if they are still zero at the end we have a logical error 
01779          */
01780         Rate_max_nonIonizing = 0.;
01781         Freq_nonIonizing = 0.;
01782         Opac_nonIonizing = 0.;
01783 
01784         /* this must be reset to val >= 0 */
01785         *opacm = -1.;
01786         *freqm = -1.;
01787 
01788         /* do up to carbon photo edge if carbon is turned on */
01789         /* >>>chng 00 apr 07, add test for whether element is turned on */
01790         if( dense.lgElmtOn[ipCARBON] )
01791         {
01792                 /* carbon is turned on, use carbon 1 edge */
01793                 ipHi = Heavy.ipHeavy[ipCARBON][0] - 1;
01794         }
01795         else
01796         {
01797                 /* carbon turned off, use hydrogen Balmer continuum */
01798                 ipHi = iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH2s]-1;
01799         }
01800 
01801         /* start at plasma frequency */
01802         for( i=rfield.ipPlasma; i < ipHi; i++ )
01803         {
01804                 /* this does not have grain opacity since grains totally passive
01805                  * at energies smaller than CI edge */
01806                 if( rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*(opac.opacity_abs[i] - 
01807                   gv.dstab[i]*dense.gas_phase[ipHYDROGEN]) > Rate_max_nonIonizing )
01808                 {
01809                         Rate_max_nonIonizing = rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*
01810                           (opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN]);
01811                         Freq_nonIonizing = rfield.anu[i];
01812                         Opac_nonIonizing = opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN];
01813                 }
01814         }
01815 
01816         /* not every continuum extends beyond C edge-this whole loop can add to zero
01817          * use total opacity here
01818          * test is to put in fir continuum if free free heating is important */
01819         if( CoolHeavy.brems_heat_total/thermal.htot > 0.05 )
01820         {
01821                 /* this is index for energy where cloud free free optical depth is unity,
01822                  * is zero if no freq are opt thin */
01823                 iplow = MAX2(1 , rfield.ipEnergyBremsThin );
01824         }
01825         else
01826         {
01827                 /* >>>chng 00 apr 07, from Heavy.ipHeavy[0][5] to ipHi defined above, since
01828                  * would crash if element not defined */
01829                 iplow = ipHi;
01830         }
01831         /* do not go below plasma frequency */
01832         iplow = MAX2( rfield.ipPlasma , iplow );
01833 
01834         /* this energy range from carbon edge to hydrogen edge */
01835         limit = MIN2(rfield.nflux,iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1);
01836         for( i=iplow; i < limit; i++ )
01837         {
01838                 if( rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*(opac.opacity_abs[i] - 
01839                   gv.dstab[i]*dense.gas_phase[ipHYDROGEN]) > Rate_max_nonIonizing )
01840                 {
01841                         Rate_max_nonIonizing = rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*
01842                           (opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN]);
01843                         Freq_nonIonizing = rfield.anu[i];
01844                         Opac_nonIonizing = opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN];
01845                 }
01846         }
01847 
01848         /* variables to check continuum interactions over Lyman continuum */
01849         Rate_max_Hion = 0.;
01850         Opac_Hion = 0.;
01851         FreqH = 0.;
01852 
01853         /* not every continuum extends beyond 1 Ryd-this whole loop can add to zero */
01854         for( i=iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1; i < rfield.nflux; i++ )
01855         {
01856                 if( rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*(opac.opacity_abs[i] - 
01857                   gv.dstab[i]*dense.gas_phase[ipHYDROGEN]) > Rate_max_Hion )
01858                 {
01859                         /* Rate_max_Hion = anu(i)*flux(i)/widflx(i)*opac(i) */
01860                         Rate_max_Hion = rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*
01861                           (opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN]);
01862                         FreqH = rfield.anu[i];
01863                         Opac_Hion = opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN];
01864                 }
01865         }
01866 
01867         /* use Lyman continuum if its opacity is larger than non-h ion */
01868         if( Rate_max_nonIonizing < 1e-30  && Opac_Hion > SMALLFLOAT )
01869         {
01870                 /* this happens for laser source - use Lyman continuum */
01871                 *opacm = Opac_Hion;
01872                 *freqm = FreqH;
01873         }
01874         /* >>chng 05 aug 03, add last test on Opac_Hion for case where we go very
01875          * deep and very little radiation is left */
01876         else if( Opac_Hion > Opac_nonIonizing && Rate_max_Hion/Rate_max_nonIonizing > 1e-10 && Opac_Hion > SMALLFLOAT )
01877         {
01878                 /* use Lyman continuum */
01879                 *opacm = Opac_Hion;
01880                 *freqm = FreqH;
01881         }
01882         else
01883         {
01884                 /* not much rate in the Lyman continuum, stick with low energy */
01885                 *opacm = Opac_nonIonizing;
01886                 *freqm = Freq_nonIonizing;
01887         }
01888 
01889 #       if 0
01890         /*>>chng 06 aug 12, do not let zones become optically thick to
01891         * peak of dust emission - one of NA's vastly optically thick dust clouds
01892         * did not conserve total energy - very opticallly thick to ir dust emission
01893         * so ir was absorbed and reemitted many times
01894         * this makes sure the cells remain optically thin to dust emission */
01895         if( gv.lgDustOn && gv.lgGrainPhysicsOn )
01896         {
01897                 double fluxGrainPeak = -1.;
01898                 long int ipGrainPeak = -1;
01899                 long int ipGrainPeak2 = -1;
01900 
01901                 i = 0;
01902                 while( rfield.anu[i] < 0.0912  && i < (rfield.nflux-2) )
01903                 {
01904                         /* >>chng 06 aug 23.  Only want to do this test for the IR dust continuum, therefore only 
01905                          * check on optical depth for wavelengths greater than 1 micron */
01906                         if( gv.GrainEmission[i]*rfield.anu[i]*opac.opacity_abs[i] > fluxGrainPeak )
01907                         {
01908                                 ipGrainPeak = i;
01909                                 fluxGrainPeak = gv.GrainEmission[i]*rfield.anu[i]*opac.opacity_abs[i];
01910                         }
01911                         ++i;
01912                 }
01913                 ASSERT( fluxGrainPeak>=0. && ipGrainPeak >=0 );
01914 
01915                 /* >>chng 06 aug 23.  We have just found the wavelength and flux at the peak of the IR emission.
01916                  * Now we want to find the wavelength, short of the peak, which corresponds to 5% of the 
01917                  * peak emission.  This wavelength will be where the code checks to make sure the zone has
01918                  * not become optically thick.  Since dust opacity generally decreases with wavelength,  this assures that
01919                  * the optical depths are small for all wavelengths where the flux is appreciable.  Tests show that
01920                  * this allows for flux/luminosity conservation to within ~5% for an AV of 1e4 mag and at least 2 iterations*/
01921                 i = ipGrainPeak;
01922                 /* >>chng 06 aug 23.  Only want to do this test for the IR dust continuum, therefore only 
01923                  * check on opacities for wavelengths shortward of the peak and greater than 1 micron 
01924                  * this routine can be called with the dust emission totally zero - it is only evaluated 
01925                  * late to save time.  don't do the following tests if peak dust emission is zero */
01926                 if( fluxGrainPeak > 0. )
01927                 {
01928                         while( rfield.anu[i] < 0.0912 && i < (rfield.nflux-2) )
01929                         {
01930                                 /* find wavelength where flux = 5% of peak flux, shortward of the peak */
01931                                 if( gv.GrainEmission[i]*rfield.anu[i]*opac.opacity_abs[i] > 0.05*fluxGrainPeak )
01932                                 {
01933                                         /* This will be the array number and flux at 10% of the peak */
01934                                         ipGrainPeak2 = i;
01935                                 }
01936                                 ++i;
01937                         }
01938                         ASSERT( ipGrainPeak2>=0 );
01939 
01940                         /* use this as a limit to the dr if opacity is larger */
01941                         if( opac.opacity_abs[ipGrainPeak2] > *opacm )
01942                         {
01943                                 /* scale opacity by factor of 10, which further decreases the zone thickness*/
01944                                 *opacm = opac.opacity_abs[ipGrainPeak2]*10.;
01945                                 *freqm = rfield.anu[ipGrainPeak2];
01946                                 /*fprintf(ioQQQ,"DEBUT opac peak %.2e %.2e \n",
01947                                         *opacm , *freqm );*/
01948                         }
01949                 }
01950         }
01951 #       endif
01952 
01953         {
01954                 /* following should be set true to print contributors */
01955                 enum {DEBUG_LOC=false};
01956                 if( DEBUG_LOC )
01957                 {
01958                         fprintf(ioQQQ,"conratedebug \t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n", 
01959                         Rate_max_nonIonizing,Freq_nonIonizing,Opac_nonIonizing,
01960                         Rate_max_Hion,FreqH ,Opac_Hion,*freqm,*opacm
01961                         );
01962 
01963                 }
01964         }
01965 
01966         /* these were set to -1 at start, and must have been reset in one of the
01967          * two loops.  Logic error if still <0. */
01968         /* >>chng 05 aug 03, change logic to -1 on entry and check at least zero
01969          * here - will be zero if NO radiation field exists, perhaps since totally
01970          * attenuated */
01971         ASSERT( *opacm >= 0. && *freqm >= 0. );
01972         return;
01973 }
01974 
01975 /*GrainRateDr called by radius_next to find grain heating rate dr */
01976 STATIC void GrainRateDr(double *grfreqm, 
01977   double *gropacm)
01978 {
01979         long int i, 
01980           iplow;
01981         double xMax;
01982 
01983         DEBUG_ENTRY( "GrainRateDr()" );
01984 
01985         /* in all following changed from anu2 to anu  july 25 95
01986          *
01987          * find maximum continuum interaction rate */
01988 
01989         /* not every continuum extends beyond C edge-this whole loop can add to zero
01990          * use total opacity here
01991          * test is to put in fir continuum if free free heating is important */
01992         if( CoolHeavy.brems_heat_total/thermal.htot > 0.05 )
01993         {
01994                 /* this is pointer to energy where cloud free free optical depth is unity,
01995                  * is zero if no freq are opt thin */
01996                 iplow = MAX2(1 , rfield.ipEnergyBremsThin );
01997         }
01998         else
01999         {
02000                 /* do up to carbon photo edge if carbon is turned on */
02001                 /* >>>chng 00 apr 07, add test for whether element is turned on */
02002                 if( dense.lgElmtOn[ipCARBON] )
02003                 {
02004                         /* carbon is turned on, use carbon 1 edge */
02005                         iplow = Heavy.ipHeavy[ipCARBON][0];
02006                 }
02007                 else
02008                 {
02009                         /* carbon truned off, use hydrogen balmer continuum */
02010                         iplow = iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH2s];
02011                 }
02012 
02013         }
02014 
02015         xMax = -1.;
02016         /* integrate up to H edge */
02017         for( i=iplow-1; i < Heavy.ipHeavy[ipHYDROGEN][0]; i++ )
02018         {
02019                 if( rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*opac.opacity_abs[i] > xMax )
02020                 {
02021                         xMax = rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*
02022                           opac.opacity_abs[i];
02023                         *grfreqm = rfield.anu[i];
02024                         *gropacm = opac.opacity_abs[i];
02025                 }
02026         }
02027         /* integrate up to heii edge if he is turned on,
02028          * this logic will not make sense if grains on but he off, which in itself makes no sense*/
02029         if( dense.lgElmtOn[ipHELIUM] )
02030         {
02031                 for( i=Heavy.ipHeavy[ipHYDROGEN][0]-1; i < Heavy.ipHeavy[ipHELIUM][1]; i++ )
02032                 {
02033                         if( rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*opac.opacity_abs[i] > xMax )
02034                         {
02035                                 xMax = rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*
02036                                   opac.opacity_abs[i];
02037                                 *grfreqm = rfield.anu[i];
02038                                 *gropacm = opac.opacity_abs[i];
02039                         }
02040                 }
02041         }
02042 
02043         /* possible that there is NO ionizing radiation, in extreme cases,
02044          * if so then xMax will still be negative */
02045         if( xMax <= 0. )
02046         {
02047                 *gropacm = 0.;
02048                 *grfreqm = 0.;
02049         }
02050         return;
02051 }
02052 
02053 #undef WIND_CHNG_VELOCITY_RELATIVE
02054 
02055 #if 0
02056 /*ChkRate called by radius_next to check how rates of destruction of various species changes */
02057 STATIC void ChkRate(
02058           /* element number on C scale */
02059           long int nelem, 
02060           /* change in destruction rate */
02061           double *dDestRate, 
02062           /* old and new destruction rates */
02063           double *DestRateOld,
02064           double *DestRateNew,
02065           /* stage of ionization on the physical scale */
02066           long int *istage)
02067 {
02068         long int i;
02069 
02070         double average, 
02071           dDest;
02072 
02073         static double OldDest[LIMELM][LIMELM];
02074 
02075         DEBUG_ENTRY( "ChkRate()" );
02076 
02077         /* return if this element is not turned on */
02078         if( !dense.lgElmtOn[nelem] )
02079         {
02080                 *dDestRate = 1e-3;
02081                 *DestRateOld = 1e-3;
02082                 *DestRateNew = 1e-3;
02083                 *istage = 0;
02084                 return;
02085         }
02086 
02087         /* for first zone, and during search, we will do nothing but
02088          * still must return finite numbers */
02089         *istage = 1;
02090         *dDestRate = 0.;
02091         *DestRateOld = 0.;
02092         *DestRateNew = 0.;
02093         *dDestRate = 0.;
02094 
02095         if( nzone <= 1 )
02096         {
02097                 for( i=0; i < nelem+1; i++ )
02098                 {
02099                         OldDest[nelem][i] = ionbal.RateIonizTot[nelem][i];
02100                 }
02101         }
02102 
02103         else if( dense.xIonDense[nelem][0]/dense.gas_phase[nelem] <  0.9 )
02104         {
02105                 /* do not use this method if everything is atomic */
02106                 for( i=0; i < (nelem); i++ )
02107                 {
02108                         /* last check below, .5 chosen so that do not key off
02109                          * predominantly neutral species where self-opacity
02110                          * could cause oscillations */
02111                         if( ((dense.xIonDense[nelem][i]/dense.gas_phase[nelem] > 0.01 && 
02112                                 dense.xIonDense[nelem][i]/dense.gas_phase[nelem] < 0.9) && 
02113                                 dense.xIonDense[nelem][i+1]/dense.gas_phase[nelem] > .05) && 
02114                                 OldDest[nelem][i] > 0. )
02115                         {
02116                                 /* last check on old dest in case just lowered ionization
02117                                  * stage, so no history */
02118                                 /* check that rate is positive */
02119                                 if( ionbal.RateIonizTot[nelem][i] <= 0. )
02120                                 {
02121                                         fprintf( ioQQQ, " ChkRate gets insane destruction rate for ion%4ld%4ld%10.2e\n", 
02122                                           nelem+1, i, ionbal.RateIonizTot[nelem][i] );
02123                                         cdEXIT(EXIT_FAILURE);
02124                                 }
02125 
02126                                 /* do not consider unless of middling ionization, and
02127                                  * rate is going down (to prevent dr osciallating)
02128                                  * no absolute value in following since do not want to
02129                                  * consider cases where ionization rate increases */
02130                                 average = (OldDest[nelem][i] + ionbal.RateIonizTot[nelem][i])* 0.5;
02131 
02132                                 dDest = (OldDest[nelem][i] - ionbal.RateIonizTot[nelem][i])/ average;
02133                                 /* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ term + if rate going down */
02134 
02135                                 if( *dDestRate < dDest )
02136                                 {
02137                                         /* biggest rate so far, remember change in rates and ionization stage */
02138                                         *dDestRate = dDest;
02139                                         *istage = i+1;
02140                                         *DestRateOld = OldDest[nelem][i];
02141                                         *DestRateNew = ionbal.RateIonizTot[nelem][i];
02142                                 }
02143 
02144                         }
02145                         OldDest[nelem][i] = ionbal.RateIonizTot[nelem][i];
02146                 }
02147         }
02148         return;
02149 }
02150 #endif
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

Generated for cloudy by doxygen 1.7.6.1