Commits

Matthew Turk committed cfafd76

Removing:

* PrimordialChemistrySolver, which has been unused for some time and has been
discouraged for a similar amount of time.
* Experiments with CVODE, which were not finished and should have not been
pushed.

  • Participants
  • Parent commits 75fe3d2

Comments (0)

Files changed (22)

 
    int SolveRateAndCoolEquations(int RTCoupledSolverIntermediateStep);
 
-/* Solve the joint rate and radiative cooling/heating equations using MTurk's Solver */
-
-   int SolveHighDensityPrimordialChemistry();
-#ifdef USE_CVODE
-   int SolvePrimordialChemistryCVODE();
-#endif
-
 /* Compute densities of various species for RadiationFieldUpdate. */
 
    int RadiationComputeDensities(int level);

src/enzo/Grid_MultiSpeciesHandler.C

   LCAPERF_START("grid_MultiSpeciesHandler");
 
   if (MultiSpecies && RadiativeCooling ) {
-    if((MultiSpecies == 3) && (PrimordialChemistrySolver > 0)) {
-      if (PrimordialChemistrySolver == 1) {
-	this->SolveHighDensityPrimordialChemistry();
-      } else if (PrimordialChemistrySolver == 2) {
-#ifdef USE_CVODE
-	this->SolvePrimordialChemistryCVODE();
-#else
-	ENZO_FAIL("You have asked for the CVODE solver but CVODE not enabled!");
-#endif
-      }
-    } // ENDIF PrimordialChemistrySolver
-    else {
-      int RTCoupledSolverIntermediateStep = FALSE;
-      this->SolveRateAndCoolEquations(RTCoupledSolverIntermediateStep);
-    }
+    int RTCoupledSolverIntermediateStep = FALSE;
+    this->SolveRateAndCoolEquations(RTCoupledSolverIntermediateStep);
   } else {
     if (MultiSpecies)
       this->SolveRateEquations();

src/enzo/Grid_SolveHighDensityPrimordialChemistry.C

-/***********************************************************************
-/
-/  GRID CLASS (SOLVE THE COOLING/HEATING AND RATE EQUATIONS)
-/
-/  written by: Greg Bryan
-/  date:       October, 1996
-/  modified1:  July, 2005 to solve cool and rate equations simultaneously
-/  modified2:  2007-2008 to use new solver
-/
-/  PURPOSE:
-/
-/  RETURNS:
-/    SUCCESS or FAIL
-/
-************************************************************************/
-
-#include <stdio.h>
-#include <math.h>
-#include "ErrorExceptions.h"
-#include "macros_and_parameters.h"
-#include "typedefs.h"
-#include "global_data.h"
-#include "Fluxes.h"
-#include "GridList.h"
-#include "ExternalBoundary.h"
-#include "Grid.h"
-#include "fortran.def"
-#include "CosmologyParameters.h"
-
-/* This parameter controls whether the cooling function recomputes
-   the metal cooling rates.  It is reset by RadiationFieldUpdate. */
-
-/* function prototypes */
-
-int CosmologyComputeExpansionFactor(FLOAT time, FLOAT *a, FLOAT *dadt);
-int GetUnits(float *DensityUnits, float *LengthUnits,
-	     float *TemperatureUnits, float *TimeUnits,
-	     float *VelocityUnits, FLOAT Time);
-int RadiationFieldCalculateRates(FLOAT Time);
-int FindField(int field, int farray[], int numfields);
-double ReturnWallTime();
-extern "C" void FORTRAN_NAME(primordial_solver)(
-	float *d, float *e, float *ge, float *u, float *v, float *w, float *de,
-	float *HI, float *HII, float *HeI, float *HeII, float *HeIII,
-	int *in, int *jn, int *kn, int *nratec, int *iexpand, 
-           hydro_method *imethod,
-        int *idual, int *ispecies, int *idim,
-	int *is, int *js, int *ks, int *ie, int *je, int *ke, int *ih2co, 
-	   int *ipiht,
-	float *dt, float *aye, float *temstart, float *temend,
-	float *utem, float *uxyz, float *uaye, float *urho, float *utim, float *uvel,
-	float *eta1, float *eta2, float *gamma, float *fh, float *dtoh,
-	float *k1a, float *k2a, float *k3a, float *k4a, float *k5a, 
-	   float *k6a, float *k7a, float *k8a, float *k9a, float *k10a,
-	float *k11a, float *k12a, float *k13a, float *k13dda, float *k14a, 
-           float *k15a,
-        float *k16a, float *k17a, float *k18a, float *k19a, float *k21a,
-            float *k22a, float *k23a,
-	float *k24, float *k25, float *k26, float *k27, float *k28, float *k29,
-	   float *k30, float *k31,
-	float *k50a, float *k51a, float *k52a, float *k53a, float *k54a,
-	   float *k55a, float *k56a,
-	float *ceHIa, float *ceHeIa, float *ceHeIIa, float *ciHIa, 
-	   float *ciHeIa, 
-	float *ciHeISa, float *ciHeIIa, float *reHIIa, float *reHeII1a, 
-	float *reHeII2a, float *reHeIIIa, float *brema, float *compa,
-	float *comp_xraya, float *comp_temp, 
-           float *piHI, float *piHeI, float *piHeII,
-	float *HM, float *H2I, float *H2II, float *DI, float *DII, float *HDI,
-	float *hyd01ka, float *h2k01a, float *vibha, float *rotha, float *rotla,
-	float *gpldl, float *gphdl, float *HDltea, float *HDlowa, float *HDcool, float *ciecoa,
-	float *gaHIa, float *gaH2a, float *gaHea, float *gaHpa, float *gaela,
-	float *inutot, int *iradtype, int *nfreq, 
-	int *iradshield, float *avgsighp, float *avgsighep, float *avgsighe2p,
-    int *iciecool, int *ih2optical, int *errcode, int *omaskflag, int *threebody, int *subgridmask
-#ifdef UNUSED_TABULATED_EQ
-    ,float *HIeqtable, float *HIIeqtable, float *H2Ieqtable, 
-        int *nrhobins, int *nebins, float *rhostart, float *rhoend,
-    float *estart, float *eend
-#endif
-    );
-
-int grid::SolveHighDensityPrimordialChemistry()
-{
-  /* Return if this doesn't concern us. */
-  if (!( (MultiSpecies == 3)
-      &&  RadiativeCooling 
-      && (PrimordialChemistrySolver == 1))) return SUCCESS;
-
-  /* Return if this doesn't concern us. */
-  
-  if (ProcessorNumber != MyProcessorNumber)
-    return SUCCESS;
-
-  if (NumberOfBaryonFields == 0)
-    return SUCCESS;
-
-  this->DebugCheck("SolveHighDensityPrimordialChemistry");
-
-  /* Declarations */
-
-  int DensNum, GENum, TENum, Vel1Num, Vel2Num, Vel3Num, B1Num, B2Num, B3Num;
-  int DeNum, HINum, HIINum, HeINum, HeIINum, HeIIINum, HMNum, H2INum, H2IINum,
-      DINum, DIINum, HDINum;
-  FLOAT a = 1.0, dadt;
-    
-  /* Find fields: density, total energy, velocity1-3. */
-
-  if (this->IdentifyPhysicalQuantities(DensNum, GENum, Vel1Num, Vel2Num, 
-				       Vel3Num, TENum, B1Num, B2Num, B3Num) == FAIL) {
-        ENZO_FAIL("Error in IdentifyPhysicalQuantities.");
-  }
-
-  /* Find Multi-species fields. */
-
-  if (MultiSpecies)
-    if (IdentifySpeciesFields(DeNum, HINum, HIINum, HeINum, HeIINum, HeIIINum, 
-                      HMNum, H2INum, H2IINum, DINum, DIINum, HDINum) == FAIL) {
-            ENZO_FAIL("Error in grid->IdentifySpeciesFields.");
-    }
-
-  /* Find photo-ionization fields */
-
-  int kphHINum, kphHeINum, kphHeIINum, kdissH2INum;
-  int gammaNum;
-  IdentifyRadiativeTransferFields(kphHINum, gammaNum, kphHeINum, 
-				  kphHeIINum, kdissH2INum);
-
-  /* Compute size of the current grid. */
-
-  int size = 1;
-  for (int dim = 0; dim < GridRank; dim++) {
-    size *= GridDimension[dim];
-  }
-
-  /* Get easy to handle pointers for each variable. */
-
-  float *density     = BaryonField[DensNum];
-  float *totalenergy = BaryonField[TENum];
-  float *gasenergy   = BaryonField[GENum];
-  float *velocity1   = BaryonField[Vel1Num];
-  float *velocity2   = BaryonField[Vel2Num];
-  float *velocity3   = BaryonField[Vel3Num];
-
-  /* Compute total gas energy if using MHD */
-  if (HydroMethod == MHD_RK) {
-    totalenergy = new float[size];
-    float B2;
-    for (int n=0; n<size; n++) {
-      B2 = pow(BaryonField[B1Num][n],2) + pow(BaryonField[B2Num][n],2) + pow(BaryonField[B3Num][n],2);
-      totalenergy[n] = BaryonField[TENum][n] - 0.5*B2/BaryonField[DensNum][n];
-    }
-  }
-  else {
-    totalenergy = BaryonField[TENum];
-  }
-
-
-  /* If using cosmology, compute the expansion factor and get units. */
-
-  float TemperatureUnits = 1, DensityUnits = 1, LengthUnits = 1, 
-    VelocityUnits = 1, TimeUnits = 1, MassUnits = 1, aUnits = 1;
-
-  if (ComovingCoordinates) {
-
-    if (CosmologyComputeExpansionFactor(Time+0.5*dtFixed, &a, &dadt) 
-	== FAIL) {
-            ENZO_FAIL("Error in CosmologyComputeExpansionFactors.");
-    }
-
-    aUnits = 1.0/(1.0 + InitialRedshift);
-
-  }
-
-  if (GetUnits(&DensityUnits, &LengthUnits, &TemperatureUnits,
-	       &TimeUnits, &VelocityUnits, Time) == FAIL) {
-        ENZO_FAIL("Error in GetUnits.");
-  }
-
-  float afloat = float(a);
-
-  /* Calculate the rates due to the radiation field. */
-
-  if (RadiationFieldCalculateRates(Time+0.5*dtFixed) == FAIL) {
-        ENZO_FAIL("Error in RadiationFieldCalculateRates.");
-  }
-
-  /* Set up information for rates which depend on the radiation field. 
-     Precompute factors for self shielding (this is the cross section * dx). */
-
-  float HIShieldFactor = RadiationData.HIAveragePhotoHeatingCrossSection * 
-                         double(LengthUnits) * CellWidth[0][0];
-  float HeIShieldFactor = RadiationData.HeIAveragePhotoHeatingCrossSection * 
-                          double(LengthUnits) * CellWidth[0][0];
-  float HeIIShieldFactor = RadiationData.HeIIAveragePhotoHeatingCrossSection * 
-                           double(LengthUnits) * CellWidth[0][0];
-
-  /* Call the fortran routine to solve cooling equations. */
-
-  int RTCoupledSolverIntermediateStep = FALSE;
-  int mask = 1;
-  int ErrCode = 0, OutputFlag = 0;
-
-  FORTRAN_NAME(primordial_solver)(
-    density, totalenergy, gasenergy, velocity1, velocity2, velocity3,
-    BaryonField[DeNum], BaryonField[HINum], BaryonField[HIINum], 
-       BaryonField[HeINum], BaryonField[HeIINum], BaryonField[HeIIINum], 
-    GridDimension, GridDimension+1, GridDimension+2, 
-       &CoolData.NumberOfTemperatureBins, &ComovingCoordinates, &HydroMethod, 
-    &DualEnergyFormalism, &MultiSpecies, &GridRank,
-    GridStartIndex, GridStartIndex+1, GridStartIndex+2, 
-       GridEndIndex, GridEndIndex+1, GridEndIndex+2,
-       &CoolData.ih2co, &CoolData.ipiht,
-    &dtFixed, &afloat, &CoolData.TemperatureStart, &CoolData.TemperatureEnd,
-    &TemperatureUnits, &LengthUnits, &aUnits, &DensityUnits, &TimeUnits, &VelocityUnits,
-    &DualEnergyFormalismEta1, &DualEnergyFormalismEta2, &Gamma,
-       &CoolData.HydrogenFractionByMass, &CoolData.DeuteriumToHydrogenRatio,
-    RateData.k1, RateData.k2, RateData.k3, RateData.k4, RateData.k5, 
-       RateData.k6, RateData.k7, RateData.k8, RateData.k9, RateData.k10,
-    RateData.k11, RateData.k12, RateData.k13, RateData.k13dd, RateData.k14, 
-       RateData.k15, RateData.k16,
-    RateData.k17, RateData.k18, RateData.k19, RateData.k21, RateData.k22,
-       RateData.k23,
-    &RateData.k24, &RateData.k25, &RateData.k26, &RateData.k27,
-       &RateData.k28, &RateData.k29, &RateData.k30, &RateData.k31,
-    RateData.k50, RateData.k51, RateData.k52, RateData.k53,
-       RateData.k54, RateData.k55, RateData.k56,
-    CoolData.ceHI, CoolData.ceHeI, CoolData.ceHeII, CoolData.ciHI,
-       CoolData.ciHeI, 
-    CoolData.ciHeIS, CoolData.ciHeII, CoolData.reHII, CoolData.reHeII1, 
-    CoolData.reHeII2, CoolData.reHeIII, CoolData.brem, &CoolData.comp,
-    &CoolData.comp_xray, &CoolData.temp_xray,
-       &CoolData.piHI, &CoolData.piHeI, &CoolData.piHeII,
-    BaryonField[HMNum], BaryonField[H2INum], BaryonField[H2IINum],
-       BaryonField[DINum], BaryonField[DIINum], BaryonField[HDINum],
-    CoolData.hyd01k, CoolData.h2k01, CoolData.vibh, CoolData.roth,CoolData.rotl,
-    CoolData.GP99LowDensityLimit, CoolData.GP99HighDensityLimit, 
-       CoolData.HDlte, CoolData.HDlow, CoolData.HDcool, CoolData.cieco,
-    CoolData.GAHI, CoolData.GAH2, CoolData.GAHe, CoolData.GAHp, CoolData.GAel,
-    RadiationData.Spectrum[0], &RadiationFieldType, 
-          &RadiationData.NumberOfFrequencyBins, 
-    &RadiationData.RadiationShield, &HIShieldFactor, &HeIShieldFactor, &HeIIShieldFactor,
-    &CIECooling, &H2OpticalDepthApproximation, &ErrCode, &OutputFlag,
-           &ThreeBodyRate, &mask
-#ifdef UNUSED_TABULATED_EQ
-        ,RateData.HighDensityEquilibriumRate[0], 
-        RateData.HighDensityEquilibriumRate[1], 
-        RateData.HighDensityEquilibriumRate[2], 
-    &RateData.HighDensityNumberOfDensityBins, &RateData.HighDensityNumberOfEnergyBins,
-    &RateData.HighDensityDensityStart, &RateData.HighDensityDensityStop,
-    &RateData.HighDensityEnergyStart, &RateData.HighDensityEnergyStop
-#endif
-    );
-
-  if (ErrCode) {
-      fprintf(stdout, "GridLeftEdge = %"FSYM" %"FSYM" %"FSYM"\n",
-	      GridLeftEdge[0], GridLeftEdge[1], GridLeftEdge[2]);
-      fprintf(stdout, "GridRightEdge = %"FSYM" %"FSYM" %"FSYM"\n",
-	      GridRightEdge[0], GridRightEdge[1], GridRightEdge[2]);
-      fprintf(stdout, "GridDimension = %"ISYM" %"ISYM" %"ISYM"\n",
-	      GridDimension[0], GridDimension[1], GridDimension[2]);
-      ENZO_FAIL("Error in FORTRAN rate/cool solver!\n");
-  }
-
-  if (HydroMethod == MHD_RK) {
-    float B2, v2;
-    for (int n = 0; n < size; n++) {
-      B2 = pow(BaryonField[B1Num][n],2) + pow(BaryonField[B2Num][n],2) + pow(BaryonField[B3Num][n],2);
-
-      /* Always trust gas energy in cooling routine */
-      if (DualEnergyFormalism) {
-
-	v2 = pow(BaryonField[Vel1Num][n],2) + 
-	  pow(BaryonField[Vel2Num][n],2) + pow(BaryonField[Vel3Num][n],2);
-	BaryonField[TENum][n] = gasenergy[n] + 0.5*v2 + 0.5*B2/BaryonField[DensNum][n];
-      }
-      else {
-	BaryonField[TENum][n] = totalenergy[n] + 0.5*B2/BaryonField[DensNum][n];
-      }
-      
-    }
-    
-    delete totalenergy;
-  }
-
-  return SUCCESS;
-
-}

src/enzo/Make.config.assemble

     endif
 
 #-----------------------------------------------------------------------
-# DETERMINE CVODE SETTINGS
-#-----------------------------------------------------------------------
-
-    ERROR_CVODE = 1
-
-    # Settings to turn CVODE ON
-
-    ifeq ($(CONFIG_CVODE),yes)
-        ERROR_CVODE = 0
-        ASSEMBLE_CVODE_DEFINES = -DUSE_CVODE
-        ASSEMBLE_CVODE_INCLUDES = $(MACH_INCLUDES_CVODE)
-        ASSEMBLE_CVODE_LIBS     = $(MACH_LIBS_CVODE)
-    endif
-
-    # Settings to turn CVODE OFF
-
-    ifeq ($(CONFIG_CVODE),no)
-        ERROR_CVODE = 0
-    endif
-
-    # error if CONFIG_CVODE is incorrect
-
-    ifeq ($(ERROR_CVODE),1)
-       .PHONY: error_CVODE
-       error_CVODE:
-	$(error Illegal value '$(CONFIG_CVODE)' for $$(CONFIG_CVODE))
-    endif
-
-#-----------------------------------------------------------------------
 # DETERMINE PAPI SETTINGS
 #-----------------------------------------------------------------------
 
               $(ASSEMBLE_LCAPERF_DEFINES) \
               $(ASSEMBLE_PYTHON_DEFINES) \
               $(ASSEMBLE_NEW_PROBLEM_TYPES_DEFINES) \
-              $(ASSEMBLE_CVODE_DEFINES) \
               $(ASSEMBLE_MPI_DEFINES) \
               $(ASSEMBLE_OOC_BOUNDARY_DEFINES) \
               $(ASSEMBLE_PAPI_DEFINES) \
     	       $(ASSEMBLE_MPI_INCLUDES) \
                $(ASSEMBLE_HYPRE_INCLUDES) \
                $(ASSEMBLE_LCAPERF_INCLUDES) \
-               $(ASSEMBLE_CVODE_INCLUDES) \
                $(ASSEMBLE_PYTHON_INCLUDES) \
                $(ASSEMBLE_NEW_PROBLEM_TYPES_INCLUDES) \
                $(ASSEMBLE_PAPI_INCLUDES) \
            $(ASSEMBLE_HYPRE_LIBS) \
            $(ASSEMBLE_LCAPERF_LIBS) \
            $(ASSEMBLE_PAPI_LIBS) \
-           $(ASSEMBLE_CVODE_LIBS) \
            $(ASSEMBLE_PYTHON_LIBS) \
            $(ASSEMBLE_NEW_PROBLEM_TYPES_LIBS) \
            $(ASSEMBLE_CUDA_LIBS)

src/enzo/Make.config.objects

         Grid_ShearingBoxStratifiedInitializeGrid.o \
 	Grid_ShocksHandler.o \
 	Grid_SolveForPotential.o \
-        Grid_SolveHighDensityPrimordialChemistry.o \
-        Grid_SolvePrimordialChemistryCVODE.o \
         Grid_SolveHydroEquations.o \
 	Grid_SolveOneZoneFreefall.o \
 	Grid_SolvePPM_DE.o \
         pop3_properties.o \
         power_of_2.o \
         prefort2.o \
-        primordial_calc_derivs.o \
-        primordial_cool.o \
-        primordial_lookup_rates.o \
-        primordial_solver.o \
 	PrepareDensityField.o \
         PrepareGravitatingMassField.o \
         PrepareIsolatedGreensFunction.o \

src/enzo/Make.config.settings

 #    CONFIG_LCAPERF
 #    CONFIG_PYTHON
 #    CONFIG_NEW_PROBLEM_TYPES
-#    CONFIG_CVODE
 #    CONFIG_PAPI
 #    CONFIG_OOC_BOUNDARY
 #    CONFIG_ACCELERATION_BOUNDARY
      CONFIG_NEW_PROBLEM_TYPES = no
 
 #=======================================================================
-# CONFIG_CVODE
-#=======================================================================
-#    yes           use CVODE solver
-#    no            don't use CVODE solver
-#-----------------------------------------------------------------------
-
-     CONFIG_CVODE = no
-
-#=======================================================================
 # CONFIG_OOC_BOUNDARY
 #=======================================================================
 #    yes           use out-of-core top-grid boundary conditions

src/enzo/Make.config.targets

 	@echo "      gmake new-problem-types-yes"
 	@echo "      gmake new-problem-types-no"
 	@echo
-	@echo "   Set whether to use embedded CVODE solver"
-	@echo
-	@echo "      gmake cvode-yes"
-	@echo "      gmake cvode-no"
-	@echo
 	@echo "   Set whether to perform unigrid communication transpose"
 	@echo
 	@echo "      gmake unigrid-transpose-yes"
 	@echo "   CONFIG_PAPI  [papi-{yes,no}]                              : $(CONFIG_PAPI)"
 	@echo "   CONFIG_PYTHON  [python-{yes,no}]                          : $(CONFIG_PYTHON)"
 	@echo "   CONFIG_NEW_PROBLEM_TYPES  [new-problem-types-{yes,no}]    : $(CONFIG_NEW_PROBLEM_TYPES)"
-	@echo "   CONFIG_CVODE  [cvode-{yes,no}]                            : $(CONFIG_CVODE)"
 	@echo "   CONFIG_ECUDA  [cuda-{yes,no}]                             : $(CONFIG_ECUDA)"
 	@echo "   CONFIG_OOC_BOUNDARY  [ooc-boundary-{yes,no}]              : $(CONFIG_OOC_BOUNDARY)"
 	@echo "   CONFIG_ACCELERATION_BOUNDARY  [acceleration-boundary-{yes,no}]              : $(CONFIG_ACCELERATION_BOUNDARY)"
 
 #-----------------------------------------------------------------------
 
-VALID_CVODE = cvode-yes cvode-no
-.PHONY: $(VALID_CVODE)
-
-cvode-yes: CONFIG_CVODE-yes
-cvode-no: CONFIG_CVODE-no
-cvode-%:
-	@printf "\n\tInvalid target: $@\n\n\tValid targets: [$(VALID_CVODE)]\n\n"
-CONFIG_CVODE-%: suggest-clean
-	@tmp=.config.temp; \
-        grep -v CONFIG_CVODE $(MAKE_CONFIG_OVERRIDE) > $${tmp}; \
-        mv $${tmp} $(MAKE_CONFIG_OVERRIDE); \
-        echo "CONFIG_CVODE = $*" >> $(MAKE_CONFIG_OVERRIDE); \
-	$(MAKE)  show-config | grep CONFIG_CVODE; \
-	echo
-
-#-----------------------------------------------------------------------
-
 VALID_ECUDA = cuda-yes cuda-no
 .PHONY: $(VALID_ECUDA)
 

src/enzo/Make.mach.darwin

 LOCAL_INCLUDES_HYPRE  = 
 LOCAL_INCLUDES_PAPI   = # PAPI includes
 LOCAL_INCLUDES_CUDA   = -I/Developer/CUDA/common/inc
-LOCAL_INCLUDES_CVODE  = 
-#-I$(LOCAL_CVODE_INSTALL)/include/cvode/
 
-MACH_INCLUDES         = $(LOCAL_INCLUDES_HDF5) $(LOCAL_INCLUDES_CUDA) \
-                        $(LOCAL_INCLUDES_CVODE)
+MACH_INCLUDES         = $(LOCAL_INCLUDES_HDF5) $(LOCAL_INCLUDES_CUDA)
 MACH_INCLUDES_PYTHON  = $(LOCAL_INCLUDES_PYTHON)
 MACH_INCLUDES_MPI     = $(LOCAL_INCLUDES_MPI)
 MACH_INCLUDES_HYPRE   = $(LOCAL_INCLUDES_HYPRE)
 LOCAL_LIBS_PYTHON = -lpython2.6
 LOCAL_LIBS_PYTHON = -framework Python 
 LOCAL_LIBS_CUDA   = -L/usr/local/cuda/lib/ -lcudart
-LOCAL_LIBS_CVODE  = 
-#-L$(LOCAL_CVODE_INSTALL)/lib -lsundials_cvode \
-#                    -lsundials_nvecserial -lhdf5_hl
 
-MACH_LIBS        = $(LOCAL_LIBS_HDF5) $(LOCAL_LIBS_MACH) $(LOCAL_LIBS_CVODE)
+MACH_LIBS        = $(LOCAL_LIBS_HDF5) $(LOCAL_LIBS_MACH) 
 MACH_LIBS_PYTHON  = $(LOCAL_LIBS_PYTHON)
 MACH_LIBS_MPI    = $(LOCAL_LIBS_MPI)
 MACH_LIBS_HYPRE  = $(LOCAL_LIBS_HYPRE)

src/enzo/Make.mach.trestles

 #LOCAL_COMPILER_DIR = /opt/pgi/linux86-64/10.5
 LOCAL_COMPILER_DIR = /opt/intel/Compiler/11.1/072
 LOCAL_HYPRE_INSTALL = 
-LOCAL_CVODE_INSTALL = 
 
 # With MPI
 
 LOCAL_INCLUDES_PAPI   = # PAPI includes
 LOCAL_INCLUDES_PYTHON = -I$(LOCAL_PYTHON_INSTALL)/include/python2.7 \
                         -I$(LOCAL_PYTHON_INSTALL)/lib/python2.7/site-packages/numpy/core/include
-LOCAL_INCLUDES_CVODE  = 
 
 MACH_INCLUDES         = $(LOCAL_INCLUDES_HDF5)
 MACH_INCLUDES_PYTHON  = $(LOCAL_INCLUDES_PYTHON)
 MACH_INCLUDES_MPI     = $(LOCAL_INCLUDES_MPI)
 MACH_INCLUDES_HYPRE   = $(LOCAL_INCLUDES_HYPRE)
 MACH_INCLUDES_PAPI    = $(LOCAL_INCLUDES_PAPI)
-MACH_INCLUDES_CVODE   = $(LOCAL_INCLUDES_CVODE)
 
 #-----------------------------------------------------------------------
 # Libraries
 LOCAL_LIBS_PAPI   = # PAPI libraries
 LOCAL_LIBS_PYTHON  = -L$(LOCAL_PYTHON_INSTALL)/lib/ -lpython2.7 \
                      -lreadline -ltermcap -lutil
-LOCAL_LIBS_CVODE  = 
 
 #LOCAL_LIBS_MACH   = -L$(LOCAL_COMPILER_DIR)/lib \
 #			-lpgf90 -lpgf90_rpm1 -lpgf902 -lpgf90rtl -lpgftnrtl -lrt
 MACH_LIBS_HYPRE   = $(LOCAL_LIBS_HYPRE)
 MACH_LIBS_PAPI    = $(LOCAL_LIBS_PAPI)
 MACH_LIBS_PYTHON  = $(LOCAL_LIBS_PYTHON)
-MACH_LIBS_CVODE   = $(LOCAL_LIBS_CVODE)

src/enzo/Make.mach.triton-gnu

 LOCAL_PYTHON_INSTALL = /home/mjturk/yt-x86_64-libenzo/
 #LOCAL_COMPILER_DIR = /opt/intel/Compiler/11.1/072/bin/intel64/
 #LOCAL_HYPRE_INSTALL = /home/mjturk/intel_mx/
-#LOCAL_CVODE_INSTALL = /projects/lca-group/local-dev/
 
 # With MPI
 
 LOCAL_INCLUDES_PAPI   = # PAPI includes
 LOCAL_INCLUDES_PYTHON = -I$(LOCAL_PYTHON_INSTALL)/include/python2.6/ \
                         -I$(LOCAL_PYTHON_INSTALL)/lib/python2.6/site-packages/numpy/core/include
-LOCAL_INCLUDES_CVODE  = -I$(LOCAL_CVODE_INSTALL)/include/cvode/
 
 MACH_INCLUDES         = $(LOCAL_INCLUDES_HDF5)
 MACH_INCLUDES_PYTHON  = $(LOCAL_INCLUDES_PYTHON)
 MACH_INCLUDES_MPI     = $(LOCAL_INCLUDES_MPI)
 MACH_INCLUDES_HYPRE   = $(LOCAL_INCLUDES_HYPRE)
 MACH_INCLUDES_PAPI    = $(LOCAL_INCLUDES_PAPI)
-MACH_INCLUDES_CVODE   = $(LOCAL_INCLUDES_CVODE)
 
 #-----------------------------------------------------------------------
 # Libraries
 LOCAL_LIBS_PAPI   = # PAPI libraries
 LOCAL_LIBS_PYTHON  = $(LOCAL_PYTHON_INSTALL)/lib/python2.6/config/libpython2.6.a \
                      -lreadline -ltermcap -lutil -lpython2.6
-LOCAL_LIBS_CVODE  = -L$(LOCAL_CVODE_INSTALL)/lib -lsundials_cvode \
-                    -lsundials_nvecserial -lhdf5_hl \
-                    /projects/lca-group/local-dev/lib/liblapack.a \
-                    /projects/lca-group/local-dev//lib/libblas.a
 
 LOCAL_LIBS_MACH   = -lgfortran
 
 MACH_LIBS_HYPRE   = $(LOCAL_LIBS_HYPRE)
 MACH_LIBS_PAPI    = $(LOCAL_LIBS_PAPI)
 MACH_LIBS_PYTHON  = $(LOCAL_LIBS_PYTHON)
-MACH_LIBS_CVODE   = $(LOCAL_LIBS_CVODE)

src/enzo/Make.mach.triton-intel

 LOCAL_PYTHON_INSTALL = /home/mjturk/yt-x86_64-libenzo-intel/
 LOCAL_COMPILER_DIR = /opt/intel/Compiler/11.1/072/bin/intel64/
 LOCAL_HYPRE_INSTALL = /home/mjturk/intel_mx/
-LOCAL_CVODE_INSTALL = /projects/lca-group/local-dev/
 
 # With MPI
 
 LOCAL_INCLUDES_PAPI   = # PAPI includes
 LOCAL_INCLUDES_PYTHON = -I$(LOCAL_PYTHON_INSTALL)/include/python2.6/ \
                         -I$(LOCAL_PYTHON_INSTALL)/lib/python2.6/site-packages/numpy/core/include
-LOCAL_INCLUDES_CVODE  = -I$(LOCAL_CVODE_INSTALL)/include/cvode/
 
 MACH_INCLUDES         = $(LOCAL_INCLUDES_HDF5)
 MACH_INCLUDES_PYTHON  = $(LOCAL_INCLUDES_PYTHON)
 MACH_INCLUDES_MPI     = $(LOCAL_INCLUDES_MPI)
 MACH_INCLUDES_HYPRE   = $(LOCAL_INCLUDES_HYPRE)
 MACH_INCLUDES_PAPI    = $(LOCAL_INCLUDES_PAPI)
-MACH_INCLUDES_CVODE   = $(LOCAL_INCLUDES_CVODE)
 
 #-----------------------------------------------------------------------
 # Libraries
 #                     -lreadline -ltermcap -lutil -lpython2.6 -ldl
 LOCAL_LIBS_PYTHON = -L$(LOCAL_PYTHON_INSTALL)/lib/ -lpython2.6 \
                      -lreadline -ltermcap -lutil -ldl
-LOCAL_LIBS_CVODE  = -L$(LOCAL_CVODE_INSTALL)/lib -lsundials_cvode \
-                    -lsundials_nvecserial -lhdf5_hl \
-                    /projects/lca-group/local-dev/lib/liblapack.a \
-                    /projects/lca-group/local-dev//lib/libblas.a
 
 LOCAL_LIBS_MACH   = -limf -lifcore -lifport 
 
 MACH_LIBS_HYPRE   = $(LOCAL_LIBS_HYPRE)
 MACH_LIBS_PAPI    = $(LOCAL_LIBS_PAPI)
 MACH_LIBS_PYTHON  = $(LOCAL_LIBS_PYTHON)
-MACH_LIBS_CVODE   = $(LOCAL_LIBS_CVODE)

src/enzo/ReadParameterFile.C

     ret += sscanf(line, "RadiativeCoolingModel = %"ISYM, &RadiativeCoolingModel);
     ret += sscanf(line, "GadgetEquilibriumCooling = %"ISYM, &GadgetEquilibriumCooling);
     ret += sscanf(line, "MultiSpecies = %"ISYM, &MultiSpecies);
-    ret += sscanf(line, "PrimordialChemistrySolver = %"ISYM, &PrimordialChemistrySolver);
     ret += sscanf(line, "CIECooling = %"ISYM, &CIECooling);
     ret += sscanf(line, "H2OpticalDepthApproximation = %"ISYM, &H2OpticalDepthApproximation);
     ret += sscanf(line, "ThreeBodyRate = %"ISYM, &ThreeBodyRate);

src/enzo/SetDefaultGlobalValues.C

 
   MultiSpecies                = FALSE;             // off
   NoMultiSpeciesButColors     = FALSE;             // off
-  PrimordialChemistrySolver   = 0;
   ThreeBodyRate               = 0;                 // ABN02
   CIECooling                  = 1;
   H2OpticalDepthApproximation = 1;

src/enzo/WriteParameterFile.C

   fprintf(fptr, "RadiativeCoolingModel          = %"ISYM"\n", RadiativeCoolingModel);
   fprintf(fptr, "GadgetEquilibriumCooling       = %"ISYM"\n", GadgetEquilibriumCooling);
   fprintf(fptr, "MultiSpecies                   = %"ISYM"\n", MultiSpecies);
-  fprintf(fptr, "PrimordialChemistrySolver      = %"ISYM"\n", PrimordialChemistrySolver);
   fprintf(fptr, "CIECooling                     = %"ISYM"\n", CIECooling);
   fprintf(fptr, "H2OpticalDepthApproximation    = %"ISYM"\n", H2OpticalDepthApproximation);
   fprintf(fptr, "ThreeBodyRate                  = %"ISYM"\n", ThreeBodyRate);

src/enzo/global_data.h

 
 EXTERN int MultiSpecies;
 EXTERN int NoMultiSpeciesButColors;
-EXTERN int PrimordialChemistrySolver;
 EXTERN int ThreeBodyRate;
 EXTERN RateDataType RateData;
 EXTERN int H2FormationOnDust;

src/enzo/primordial_calc_derivs.F

-#include "fortran.def"
-#define TOLERANCE 1.d-20
-#define RTOL 5.d-1
-#define ITERTOL 1.d-4
-c -----------------------------------------------------------
-c  This routine is very simple.  It calculates the derivatives
-c  of the species.
-c
-      subroutine primordial_calc_derivs(
-     &                 in,  is,  ie,   s,  sp, sm1, dsq, dsp,  fd,
-     &                k01, k02, k03, k04, k05, k06, k07, k08, k09, k10,
-     &                ar01,ar02,ar03,ar04,ar05,ar06,ar07,ar08,ar09,ar10,
-     &                k11, k12, k13, k14, k15, k16, k17, k18, k19, k21,
-     &                ar11,ar12,ar13,ar14,ar15,ar16,ar17,ar18,ar19,ar21,
-     &                k22, k23, k24, k25, k26, k27, k28, k29, k30, 
-     &                ar22,ar23,ar24,ar25,ar26,ar27,ar28,ar29,ar30,
-     &                k31, k50, k51, k52, k53, k54, k55, k56,
-     &                ar31,ar50,ar51,ar52,ar53,ar54,ar55,ar56,
-     &                k24shield, k25shield, k26shield,
-     &                dmask, chunit, ncspecies, doupdate,
-     &                tsgamma, dtit, tscapy, checker, mins, iter,
-     &                dtcoef, dtot, rtot,
-     &                h2heatp, h2heatm, hiheatp, heiheatp, heiiheatp,
-     &                tsc, esterr, tgas,
-     &                rejected, itererror, eqmask, dom, dt, ttot )
-c -------------------------------------------------------------------
-c
-      implicit NONE
-c
-c     arguments
-c
-      integer in, is, ie, ijk, ncspecies, n
-      parameter (ijk = MAX_ANY_SINGLE_DIRECTION)
-      integer mins(ijk), iter
-      real*8  chunit, checker, dtcoef(NSPECIES)
-      real*8 s(ijk,NSPECIES), sp(ijk,NSPECIES), sm1(ijk,NSPECIES),
-     &       dsq(ijk,NSPECIES), dsp(ijk,NSPECIES), tempde, peqHq, peqHp,
-     &       temp1, temp2, temp3, temp4, tsgamma(ijk),
-     &       tscapy(ijk,NSPECIES), dtot(ijk), rtot(ijk), h2heatp(ijk),
-     &       h2heatm(ijk), hiheatp(ijk), heiheatp(ijk), heiiheatp(ijk),
-     &       tsc(ijk), esterr(ijk), qbar, qdot, tbar, prevval,
-     &       itererror, fd(ijk)
-      logical dmask(ijk), doupdate, rejected(ijk), eqmask(ijk)
-      real dtit(ijk), dom, dt, ttot(ijk), tgas(ijk)
-c
-c     Rate values
-c
-      real*8 r01, r02, r03, r04, r05,
-     &       r06, r07, r08, r09, r10,
-     &       r11, r12, r13, r14, r15,
-     &       r16, r17, r18, r19, r21,
-     &       r22, r23, r50, r51, r52,
-     &       r53, r54, r55, r56, r61, r62
-      real ar01(in), ar02(in), ar03(in), ar04(in), ar05(in),
-     &     ar06(in), ar07(in), ar08(in), ar09(in), ar10(in),
-     &     ar11(in), ar12(in), ar13(in), ar14(in), ar15(in),
-     &     ar16(in), ar17(in), ar18(in), ar19(in), ar21(in),
-     &     ar22(in), ar23(in), ar50(in), ar51(in), ar52(in),
-     &     ar53(in), ar54(in), ar55(in), ar56(in),
-     &     ar24, ar25, ar26, ar27, ar28, ar29, ar30, ar31
-      real k01(in), k02(in), k03(in), k04(in), k05(in),
-     &     k06(in), k07(in), k08(in), k09(in), k10(in),
-     &     k11(in), k12(in), k13(in), k14(in), k15(in),
-     &     k16(in), k17(in), k18(in), k19(in), k21(in),
-     &     k22(in), k23(in), k50(in), k51(in), k52(in),
-     &     k53(in), k54(in), k55(in), k56(in), 
-     &     k24shield(in), k25shield(in), k26shield(in),
-     &     k24, k25, k26, k27, k28, k29, k30, k31
-c
-c     locals
-c
-      integer i, retloc
-      logical stopthepresses, flipdmask
-
-      integer HI,HII,de,HeI,HeII,HeIII,H2I,H2II,HM,ge,DI,DII,HDI
-      parameter(HI=1,HII=2,HeI=3,HeII=4,HeIII=5,
-     &          H2I=6,ge=7,DI=8,DII=9,HDI=10,
-     &          de=11,HM=12,H2II=13 ) ! eq species
-
-        stopthepresses = .false.
-        flipdmask=.false.
-        if(itererror.gt.1)flipdmask=.true.
-        itererror = -1.d0
-
-        do i = is+1, ie+1
-        if(dmask(i).eqv..true.) then
-          rejected(i) = .false.
-          temp3 = tiny
-          mins(i) = -1
-
-#include "primordial_reactions.inc"
-c
-c-------------HM-----------------------------------------------------
-c
-        n = HM
-        prevval = s(i,n)
-        dsq(i,n) = 
-     &       + r07
-        dsp(i,n) = 
-     &       + r08
-     &       + r14
-     &       + r15
-     &       + r16
-     &       + r17
-     &       + r19
-        if(doupdate.eqv..true.)then
-          s(i,n) = s(i,n) * (dsq(i,n)/dsp(i,n))
-          s(i,n) = max(s(i,n),tiny*rtot(i))
-#include "primordial_reactions.inc"
-        endif
-
-c
-c-------------H2II---------------------------------------------------
-c
-        n = H2II
-        prevval = s(i,n)
-        dsq(i,n) = (
-     &       + r09
-     &       + r11
-     &       + r17 ) * 2.0d0
-        dsp(i,n) = (
-     &       + r10
-     &       + r18
-     &       + r19 ) * 2.0d0
-        if(doupdate.eqv..true.)then
-          s(i,n) = s(i,n) * (dsq(i,n)/dsp(i,n))
-          s(i,n) = max(s(i,n),tiny*rtot(i))
-#include "primordial_reactions.inc"
-        endif
-
-c Notes: r10 and r11 included in HI and HII cause problems with the H2 and
-c recombination rates.
-c
-c-------------HII----------------------------------------------------
-c
-        n = HII
-        prevval = s(i,n)
-        dsq(i,n) = 
-     &       + r01
-     &       + r10 !
-     &       + r51 !
-        dsp(i,n) = 
-     &       + r02
-     &       + r09 !
-     &       + r11 !
-     &       + r16 !
-     &       + r17 !
-     &       + r50 !
-        if(doupdate.eqv..true..and.eqmask(i).eqv..false.)then
-          s(i,n) = (dsq(i,n)*tsgamma(i)*dtit(i)+tscapy(i,n))
-     &            / (1.0d0+tsgamma(i)*dtit(i)*dsp(i,n)/s(i,n))
-          temp2 = abs((2.0d0/(tsc(i)+1.0d0)) *
-     &         (tsc(i)*s(i,n) - (1.0d0+tsc(i))*sp(i,n) + sm1(i,n)))
-     &      /  (RTOL*dtcoef(n)*sp(i,n))
-          itererror = max(itererror, abs(s(i,n)-prevval)/prevval)
-          if(temp2.gt.temp3)then
-            temp3 = temp2
-            mins(i) = n
-          endif
-#include "primordial_reactions.inc"
-        endif
-
-c
-c-------------electrons----------------------------------------------
-c
-        n = de
-        prevval = s(i,n)
-        dsq(i,n) = 
-     &       + r01
-     &       + r03
-     &       + r05
-     &       + r08 !
-     &       + r14 !
-     &       + r15 !
-     &       + r17 !
-     &       + r61 !
-        dsp(i,n) = 
-     &       + r02
-     &       + r04
-     &       + r06
-     &       + r07 !
-     &       + r18 !
-     &       + r62 !
-        if(doupdate.eqv..true.)then
-c        s(i,n) = (dsq(i,n)*dtit(i)+s(i,n))
-c     &         / (1.0d0+dtit(i)*dsp(i,n)/s(i,n))
-        s(i,de) = s(i,HII) + s(i,HeII)/4.0d0 + s(i,HeIII)/2.0d0
-     &           + s(i,H2II)/2.0d0 - s(i,HM) + s(i,DII)/2.0d0
-        s(i,de) = max(tiny,s(i,de))
-#include "primordial_reactions.inc"
-        endif
-
-c
-c-------------HI-----------------------------------------------------
-c
-        n = HI
-        prevval = s(i,n)
-        dsq(i,n) = (
-     &       + r02
-     &       + r11 !
-     &       + r12  * 2.0d0
-     &       + r14 !
-     &       + r15 !
-     &       + r16  * 2.0d0 !
-     &       + r18  * 2.0d0 !
-     &       + r19 !
-     &       + r13  * 2.0d0
-     &       + r23  * 2.0d0
-     &       + r50 !
-     &       + r54 !
-     &       )
-        dsp(i,n) = (
-     &       + r01
-     &       + r07 !
-     &       + r08 !
-     &       + r09 !
-     &       + r10 !
-     &       + r21  * 2.0d0
-     &       + r22  * 2.0d0
-     &       + r51 !
-     &       + r55 !
-     &       )
-        if(doupdate.eqv..true..and.eqmask(i).eqv..false.)then
-          s(i,n) = (dsq(i,n)*tsgamma(i)*dtit(i)+tscapy(i,n))
-     &            / (1.0d0+tsgamma(i)*dtit(i)*dsp(i,n)/s(i,n))
-          temp2 = abs((2.0d0/(tsc(i)+1.0d0)) *
-     &         (tsc(i)*s(i,n) - (1.0d0+tsc(i))*sp(i,n) + sm1(i,n)))
-     &      /  (RTOL*dtcoef(n)*sp(i,n))
-          itererror = max(itererror, abs(s(i,n)-prevval)/prevval)
-          if(temp2.gt.temp3)then
-            temp3 = temp2
-            mins(i) = n
-          endif
-#include "primordial_reactions.inc"
-        endif
-        hiheatp(i) = 1.0d0 * (r01)
-
-c
-c-------------H2I----------------------------------------------------
-c
-        n = H2I
-        prevval = s(i,n)
-        dsq(i,n) = (
-     &       + r08
-     &       + r10
-     &       + r19
-     &       + r21
-     &       + r22
-     &       + r53 !
-     &       + r55 !
-     &       ) * 2.0d0
-        dsp(i,n) = (
-     &       + r11
-     &       + r12
-     &       + r13
-     &       + r23
-     &       + r52 !
-     &       + r54 !
-     &       ) * 2.0d0
-! With quantitiyes in number density:
-! H2_eq = (sqrt(k13)*sqrt(4*k22*rho+k13)-k13)/(2*k22)
-        if(doupdate.eqv..true..and.eqmask(i).eqv..false.)then
-          s(i,n) = (dsq(i,n)*tsgamma(i)*dtit(i)+tscapy(i,n))
-     &            / (1.0d0+tsgamma(i)*dtit(i)*dsp(i,n)/s(i,n))
-          temp2 = abs((2.0d0/(tsc(i)+1.0d0)) *
-     &         (tsc(i)*s(i,n) - (1.0d0+tsc(i))*sp(i,n) + sm1(i,n)))
-     &      /  (RTOL*dtcoef(n)*s(i,n))
-          itererror = max(itererror, abs(s(i,n)-prevval)/prevval)
-          if(temp2.gt.temp3)then
-            temp3 = temp2
-            mins(i) = n
-          endif
-        endif
-#include "primordial_reactions.inc"
-        if(s(i,H2I)/fd(i).gt.1e-7)then
-          h2heatp(i) = 2.0d0 * (r21 + r22)
-          h2heatm(i) = 2.0d0 * (r13 + r23)
-        endif
-c        if((tgas(i).gt.8000).and.(rtot(i)*dom.gt.1e8))then
-c          h2heatp(i) = max((s(i,H2I)-sp(i,H2I))/dtit(i),0.0d0)
-c          h2heatm(i) = max((sp(i,H2I)-s(i,H2I))/dtit(i),0.0d0)
-c        endif
-
-c
-c-------------HeI----------------------------------------------------
-c
-        n = HeI
-        prevval = s(i,n)
-        dsq(i,n) = (
-     &       + r04
-     &       ) * 4.0d0
-        dsp(i,n) = (
-     &       + r03
-     &       ) * 4.0d0
-        if(doupdate.eqv..true.)then
-        s(i,n) = (dsq(i,n)*tsgamma(i)*dtit(i)+tscapy(i,n))
-     &          / (1.0d0+tsgamma(i)*dtit(i)*dsp(i,n)/s(i,n))
-        s(i,n) = max(s(i,n), tiny)
-#include "primordial_reactions.inc"
-        temp2 = abs((2.0d0/(tsc(i)+1.0d0)) *
-     &       (tsc(i)*s(i,n) - (1.0d0+tsc(i))*sp(i,n) + sm1(i,n)))
-     &    /  (RTOL*dtcoef(n)*sp(i,n))
-c        itererror = max(itererror, abs(s(i,n)-prevval)/prevval)
-        if(temp2.gt.temp3)then
-          temp3 = temp2
-          mins(i) = n
-        endif
-        endif
-        heiheatp(i)  = 4.0d0 * (r03)
-c        write(0,*) heiheatp(i)
-
-c
-c-------------HeII---------------------------------------------------
-c
-        n = HeII
-        prevval = s(i,n)
-        dsq(i,n) = (
-     &       + r03
-     &       + r06
-     &       ) * 4.0d0
-        dsp(i,n) = (
-     &       + r04
-     &       + r05
-     &       ) * 4.0d0
-        if(doupdate.eqv..true.)then
-        s(i,n) = (dsq(i,n)*tsgamma(i)*dtit(i)+tscapy(i,n))
-     &          / (1.0d0+tsgamma(i)*dtit(i)*dsp(i,n)/s(i,n))
-        s(i,n) = max(tiny,s(i,n))
-#include "primordial_reactions.inc"
-        temp2 = abs((2.0d0/(tsc(i)+1.0d0)) *
-     &       (tsc(i)*s(i,n) - (1.0d0+tsc(i))*sp(i,n) + sm1(i,n)))
-     &    /  (RTOL*dtcoef(n)*sp(i,n))
-c        itererror = max(itererror, abs(s(i,n)-prevval)/prevval)
-        if(temp2.gt.temp3)then
-c          temp3 = temp2
-c          mins(i) = n
-        endif
-        endif
-        heiiheatp(i) = 4.0d0 * (r05)
-
-c
-c-------------HeIII--------------------------------------------------
-c
-        n = HeIII
-        prevval = s(i,n)
-        dsq(i,n) = (
-     &       + r05
-     &       ) * 4.0d0
-        dsp(i,n) = (
-     &       + r06
-     &       ) * 4.0d0
-        if(doupdate.eqv..true.)then
-        s(i,n) = (dsq(i,n)*tsgamma(i)*dtit(i)+tscapy(i,n))
-     &          / (1.0d0+tsgamma(i)*dtit(i)*dsp(i,n)/s(i,n))
-        s(i,n) = max(tiny,s(i,n))
-#include "primordial_reactions.inc"
-        temp2 = abs((2.0d0/(tsc(i)+1.0d0)) *
-     &       (tsc(i)*s(i,n) - (1.0d0+tsc(i))*sp(i,n) + sm1(i,n)))
-     &    /  (RTOL*dtcoef(n)*sp(i,n))
-c        itererror = max(itererror, abs(s(i,n)-prevval)/prevval)
-        if(temp2.gt.temp3)then
-c          temp3 = temp2
-c          mins(i) = n
-        endif
-        endif
-
-c NB:
-c   We assume that the charge exchange with the H/H+ species will completely
-c   regulate the abundance of D/D+
-c Recall:
-c   nD = (nDTot - nHD) * (nH * k51) / (nH*k51 + nH*k50)
-        temp1 = max(tiny, dtot(i) - (2.d0/3.d0) * s(i,HDI))
-c        s(i,DI)  = max(tiny*rtot(i), temp1 * 
-c     &     (s(i,HI) * k51(i)) / (s(i,HI)*k51(i) + s(i,HII)*k50(i)))
-c        s(i,DI) = s(i,DI) * (r62 + r51 + r55) / (r61 + r50 + r54)
-c        s(i,DII) = s(i,DII) * (r61 + r50) / (r62 + r51)
-c        s(i,DI) = max(temp1, tiny*rtot(i))
-c        s(i,DII) = max(temp1 - s(i,DI),tiny*rtot(i))
-c#include "primordial_reactions.inc"
-        
-c
-c-------------DII----------------------------------------------------
-c
-        n = DII
-        dsq(i,DII) = (  ! DII
-     &       + r61
-     &       + r50
-     &       + r53
-     & ) * 2.0d0
-        dsp(i,DII) = (  ! DII
-     &       + r62
-     &       + r51
-     &       + r52
-     & ) * 2.0d0
-        if(doupdate.eqv..true.)then
-        s(i,n) = (dsq(i,n)*dtit(i)+sp(i,n))
-     &         / (1.0d0+dtit(i)*dsp(i,n)/s(i,n))
-        s(i,n) = max(tiny, s(i,n))
-#include "primordial_reactions.inc"
-        endif
-
-c
-c-------------DI-----------------------------------------------------
-c
-        n = DI
-        dsq(i,DI) = (  ! DI
-     &       + r62
-     &       + r51
-     &       + r55
-     & ) * 2.0d0
-        dsp(i,DI) = (  ! DI
-     &       + r61
-     &       + r50
-     &       + r54
-     & ) * 2.0d0
-        if(doupdate.eqv..true.)then
-        s(i,n) = (dsq(i,n)*dtit(i)+sp(i,n))
-     &         / (1.0d0+dtit(i)*dsp(i,n)/s(i,n))
-        s(i,n) = max(tiny, s(i,n))
-#include "primordial_reactions.inc"
-        endif
-
-c
-c-------------HDI----------------------------------------------------
-c
-        n = HDI
-        dsq(i,HDI) = (  ! HDI
-     &       + r52
-     &       + r54
-     & ) * 3.0d0
-        dsp(i,HDI) = (  ! HDI
-     &       + r53
-     &       + r55
-     & ) * 3.0d0
-        if(doupdate.eqv..true.)then
-        s(i,n) = (dsq(i,n)*dtit(i)+sp(i,n))
-     &         / (1.0d0+dtit(i)*dsp(i,n)/s(i,n))
-        s(i,n) = max(tiny, s(i,n))
-#include "primordial_reactions.inc"
-        endif
-
-c Charge conservation constraint step only gets applied when we actually
-c update, except, we probably only want to do it once both corrections are
-c done.
-        if(doupdate.eqv..true.)then
-          s(i,de) = s(i,HII) + s(i,HeII)/4.0d0 + s(i,HeIII)/2.0d0
-     &            + s(i,H2II)/2.0d0 - s(i,HM)
-          s(i,de) = max(s(i,de),tiny)
-        endif
-
-        if(flipdmask.and.(itererror.lt.ITERTOL))dmask(i)=.false.
-        esterr(i) = temp3
-        stopthepresses = .false.
-        do n=1,NSPECIES
-          if((dsp(i,n).ne.dsp(i,n)).or.
-     &       (dsq(i,n).ne.dsq(i,n)).or.
-     &       (s(i,n).ne.s(i,n)).or.
-     &       (s(i,n).lt.0.0))then
-            write(0,*) "CALC_DERIV PROBLEM",i,n,iter,eqmask(i)
-            write(0,*)'a',dsp(i,n),dsq(i,n),s(i,n)
-            write(0,*)'b',sp(i,n),tsgamma(i),tscapy(i,n)
-            write(0,*)'c',mins(i),doupdate,dtit(i)
-            write(0,*)'d',sm1(i,n),esterr(i)
-            write(0,*)'e',rejected(i)
-            stopthepresses = .true.
-          else
-            s(i,n) = max(s(i,n),tiny*rtot(i))
-          endif
-        enddo
-        if(stopthepresses.eqv..true.) then
-          do n=1,NSPECIES
-            write(0,*) n, s(i,n), dsp(i,n), dsq(i,n), sp(i,n)
-          enddo
-          write(0,*)'done.'
-          stop
-          dmask(i) = .false.
-        endif
-
-#ifdef NOTIMPLEMENTED
-        ar01(i) = r01
-        ar02(i) = r02
-        ar03(i) = r03
-        ar04(i) = r04
-        ar05(i) = r05
-        ar06(i) = r06
-        ar07(i) = r07
-        ar08(i) = r08
-        ar09(i) = r09
-        ar10(i) = r10
-        ar11(i) = r11
-        ar12(i) = r12
-        ar13(i) = r13
-        ar14(i) = r14
-        ar15(i) = r15
-        ar16(i) = r16
-        ar17(i) = r17
-        ar18(i) = r18
-        ar19(i) = r19
-        ar20(i) = r20
-        ar21(i) = r21
-        ar22(i) = r22
-        ar23(i) = r23
-        ar50(i) = r50
-        ar51(i) = r51
-        ar52(i) = r52
-        ar53(i) = r53
-        ar54(i) = r54
-        ar55(i) = r55
-        ar56(i) = r56
-#endif
-
-        endif
-        enddo
-
-        return
-        end

src/enzo/primordial_cool.F

-#include "fortran.def"
-#include "phys_const.def"
-c=======================================================================
-c//////////////////////  SUBROUTINE COOL1D_MULTI  \\\\\\\\\\\\\\\\\\\\\\
-c
-      subroutine primordial_cool(
-     &                d, e, ge, u, v, w, de, HI, HII, HeI, HeII, HeIII,
-     &                in, jn, kn, nratec, idual, imethod,
-     &                iexpand, ispecies, idim,
-     &                is, ie, j, k, ih2co, ipiht, iter,
-     &                aye, temstart, temend,
-     &                utem, uxyz, uaye, urho, utim,
-     &                eta1, eta2, gamma,
-     &                ceHIa, ceHeIa, ceHeIIa, ciHIa, ciHeIa, 
-     &                ciHeISa, ciHeIIa, reHIIa, reHeII1a, 
-     &                reHeII2a, reHeIIIa, brema, compa, 
-     &                comp_xraya, comp_temp,
-     &                piHI, piHeI, piHeII, comp1, comp2,
-     &                HM, H2I, H2II, DI, DII, HDI, 
-     &                hyd01ka, h2k01a, vibha, rotha, rotla,
-     &                hyd01k, h2k01, vibh, roth, rotl,
-     &                gpldla, gphdla, gpldl, gphdl,
-     &                hdltea, hdlowa, hdlte, hdlow, 
-     &                hdcoola, hdcool, ciecoa, cieco, 
-     &                gaHIa, gaH2a, gaHea, gaHpa, gaela,
-     &                ceHI, ceHeI, ceHeII, ciHI, ciHeI, ciHeIS, ciHeII,
-     &                reHII, reHeII1, reHeII2, reHeIII, brem,
-     &                indixe, t1, t2, logtem, tdef, edotplus,
-     &                edotminus, tgas, tgasold, p2d,
-     &                inutot, iradtype, nfreq, 
-     &                iradshield, avgsighp, avgsighep, 
-     &                avgsighe2p, itmask, iciecool,ih2optical,ifedtgas,
-     &                gamma2, omask)!, 
-!     &                sobnstart, sobnend, sobtstart, sobtend, 
-!     &                sobnn, sobnt, sobratio, neff )
-c
-c  SOLVE RADIATIVE COOLING/HEATING EQUATIONS
-c
-c  written by: Yu Zhang, Peter Anninos and Tom Abel
-c  date:       
-c  modified1: January, 1996 by Greg Bryan; adapted to KRONOS
-c  modified2: October, 1996 by GB; moved to AMR
-c
-c  PURPOSE:
-c    Solve the energy cooling equations.
-c
-c  INPUTS:
-c    is,ie   - start and end indicies of active region (zero-based!)
-c
-c  PARAMETERS:
-c
-c-----------------------------------------------------------------------
-c
-      implicit NONE
-c
-c  Arguments
-c
-      integer in, jn, kn, is, ie, j, k, nratec, imethod, idim,
-     &        idual, iexpand, ih2co, ipiht, ispecies,
-     &        nfreq, iradshield, iradtype, iciecool,
-     &        ih2optical, ifedtgas, ijk
-      parameter (ijk = MAX_ANY_SINGLE_DIRECTION)
-
-      real    aye, temstart, temend,
-     &        utem, uxyz, uaye, urho, utim,
-     &        eta1, eta2, gamma
-      real    d(in,jn,kn), e(in,jn,kn)
-      real    u(ijk),    v(ijk),     w(ijk)
-      real*8
-     &       de(ijk),   HI(ijk),   HII(ijk),
-     &      HeI(ijk), HeII(ijk), HeIII(ijk),
-     &       ge(ijk),
-     &        HM(ijk),  H2I(ijk), H2II(ijk),
-     &        DI(ijk),  DII(ijk), HDI(ijk)
-      logical omask(ijk)
-      real    hyd01ka(nratec), h2k01a(nratec), vibha(nratec), 
-     &        rotha(nratec), rotla(nratec), gpldla(nratec),
-     &        gphdla(nratec), hdltea(nratec), hdlowa(nratec),
-     $        hdcoola(nratec, 5), ciecoa(nratec)
-      real    gaHIa(nratec), gaH2a(nratec), gaHea(nratec),
-     &        gaHpa(nratec), gaela(nratec)
-      real    ceHIa(nratec), ceHeIa(nratec), ceHeIIa(nratec),
-     &        ciHIa(nratec), ciHeIa(nratec), ciHeISa(nratec), 
-     &        ciHeIIa(nratec), reHIIa(nratec), reHeII1a(nratec), 
-     &        reHeII2a(nratec), reHeIIIa(nratec), brema(nratec)
-      real    compa, piHI, piHeI, piHeII, comp_xraya, comp_temp,
-     &        inutot(nfreq), avgsighp, avgsighep, avgsighe2p,
-     &        ciecont, mciecont, medot, gamma2(in)
-c
-c  Parameters
-c
-      double precision mh
-      real    ZSOLAR
-      parameter (mh = mass_h, ZSOLAR = 0.02041d0)
-c
-c  Locals
-c
-      integer i, m, j1, iter, iradfield, mciei
-      real dom, qq, vibl, logtem0, logtem9, dlogtem, energy, zr
-      real dt2, ttmin, comp1, comp2, scoef, acoef, HIdot,
-     &     hdlte1, hdlow1, x, nH2, nother, fudge, fH2,
-     &     gphdl1, factor, ciefudge, tau, comp_edot, comp_xray_edot
-      double precision coolunit, dbase1, tbase1, xbase1, 
-     &                 tm, lt, t3, HDLR, HDLV, cierate, ttt, turnoff
-c
-c  Slice locals
-c 
-      integer indixe(in)
-      real t1(in), t2(in), logtem(in), tdef(in), p2d(in),
-     &     tgas(in), tgasold(in), tdef_inv(in)
-      real*8 edotplus(ijk), edotminus(ijk)
-      double precision ttgas
-      real temp_table(5)
-c
-c  Cooling/heating slice locals
-c
-      double precision ceHI(in), ceHeI(in), ceHeII(in),
-     &     ciHI(in), ciHeI(in), ciHeIS(in), ciHeII(in),
-     &     reHII(in), reHeII1(in), reHeII2(in), reHeIII(in),
-     &     brem(in)
-      real hyd01k(in), h2k01(in), vibh(in), roth(in), rotl(in),
-     &     gpldl(in), gphdl(in), hdlte(in), hdlow(in), 
-     &     hdcool(in, 5), hdcoolval(in), cieco(in)
-      real gaHI(in), gaH2(in), gaHe(in), gaHp(in), gael(in),
-     &     galdl(in)
-c
-c     Iteration mask
-      
-      logical itmask(in)
-      integer inneriter
-      real oldgamma2(in)
-
-      real g2p(10)
-      data g2p/0.0493452d0, 2.28021d0, 0.115357d0,
-     &         0.114188d0,  2.90778d0, 0.689708d0,
-     &         64.2416d0, -9.36334d0, -0.377266d0, 69.8091d0/
-c
-c\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\/////////////////////////////////
-c=======================================================================
-c
-c     Set log values of start and end of lookup tables
-c
-      mciecont = 0
-      logtem0 = log(temstart)
-      logtem9 = log(temend)
-      dlogtem= (log(temend) - log(temstart))/real(nratec-1)
-c
-c     Set units
-c
-      dom      = urho*(aye**3)/mh
-      tbase1   = utim
-      xbase1   = uxyz/(aye*uaye)    ! uxyz is [x]*a      = [x]*[a]*a'        '
-      dbase1   = urho*(aye*uaye)**3 ! urho is [dens]/a^3 = [dens]/([a]*a')^3 '
-      coolunit = (uaye**5 * xbase1**2 * mh**2) / (tbase1**3 * dbase1)
-c      write(6,*) "coolunit: ", coolunit
-c      write(6,*) "dom:      ", dom
-      zr       = 1.d0/(aye*uaye) - 1.d0
-      fudge    = 1.d0
-      ciefudge = 1.d0
-
-      do i = is+1, ie+1
-        edotminus(i) = 0.d0
-        edotplus(i) = 0.d0
-      enddo
-
-c
-c     Set compton cooling coefficients (and temperature)
-c
-      if (iexpand .eq. 1) then
-         comp1 = compa * (1.d0 + zr)**4
-         comp2 = 2.73  * (1.d0 + zr)
-      else
-         comp1 = tiny
-         comp2 = tiny
-      endif
-c
-c     Compute Pressure
-c
-      do i = is+1, ie+1
-      if (itmask(i)) then
-        p2d(i) = (gamma - 1.d0)*d(i,j,k)*ge(i)
-      endif
-      enddo
-c
-c     Compute temperature
-c
-      do i = is+1, ie+1
-         if (itmask(i)) then
-            tgas(i) = 
-     &           (HeI(i) + HeII(i) + HeIII(i))/4.d0 +
-     &           HI(i) + HII(i) + de(i)
-         endif
-      enddo
-c
-c          (include molecular hydrogen, but ignore deuterium)
-c
-      if (ispecies .gt. 1) then
-         do i = is+1, ie+1
-            if (itmask(i)) then
-               tgas(i) = tgas(i) +
-     &              HM(i) + (H2I(i) + H2II(i))/2.d0
-            endif
-         enddo
-      endif
-c
-      do i = is+1, ie+1
-        if (itmask(i)) then
-           tgasold(i) = tgas(i)
-           tgas(i) = max(p2d(i)*(utem/tgas(i)), temstart)
-c           write(6,*), 'tgas ', tgas(i)
-           if(tgas(i).ne.tgas(i))then
-              write(6,*) 'cool1d_multi tgas0a:',p2d(i),tgas(i)
-              write(6,*) 'cool1d_multi tgas0b:',HI(i),HII(i),HM(i)
-              write(6,*) 'cool1d_multi tgas0c:',H2I(i)
-           endif
-        endif
-      enddo
-c   Omukai
-c
-c     Correct temperature for gamma from H2
-c
-      if (ispecies .gt. 1) then
-         do i = is+1, ie+1
-         if (itmask(i).eqv..true.) then
-           nH2 = 0.5*(H2I(i) + H2II(i))
-           nother = (HeI(i) + HeII(i) + HeIII(i))/4.d0 +
-     &          HI(i) + HII(i) + de(i) + HM(i)
-           if(ispecies.gt.2)then
-             nother = nother + (DI(i)+DII(i))/2.d0
-     &                       + (HDI(i)/3.d0)
-           endif
-#define OLDGAMMA2
-#ifdef OLDGAMMA2
-           oldgamma2(i) = gamma
-           do inneriter=1,1
-             if (nH2/nother .gt. 1.0d-3) then
-               x = 6100/tgas(i) ! not quite self-consistent
-               if (x .gt. 10.0) then
-                 gamma2(i) = 0.5d0*5.d0
-               else
-                 gamma2(i) = 0.5d0*(5.d0 + 2.d0*x**2 * exp(x)
-     &                     /(exp(x)-1.d0)**2)
-               endif
-             else
-               gamma2(i) = 2.5d0
-             endif
-             gamma2(i) = 1.d0 + (nH2 + nother)/
-     &            (nH2*gamma2(i) + nother/(gamma-1.d0))
-             tgas(i) = tgas(i) * (gamma2(i) - 1.d0)/(oldgamma2(i) -1.d0)
-             oldgamma2(i) = gamma2(i)
-           enddo
-#else
-           oldgamma2(i) = gamma
-           do inneriter=1,1
-             x = log10(tgas(i))
-             if(nH2/nother.lt.-1.0d-3.or.tgas(i).lt.-50.d0)then
-               gamma2(i) = 2.5 !gamma
-             else
-               gamma2(i) = g2p(1)*exp(-(x-g2p(2))**2.d0/g2p(3))
-     &                   + g2p(4)*exp(-(x-g2p(5))**2.d0/g2p(6))
-     &                   + exp(-g2p(7)*(x**g2p(8)))
-     &                     *(g2p(9)+x**(-1.0d0*g2p(10)))
-     &                   + 5.d0/3.d0
-             endif
-             gamma2(i) = 1.d0 + (nH2 + nother)/
-     &            (nH2/(gamma2(i)-1.d0) + nother/(gamma-1.d0))
-             x = tgas(i)
-             tgas(i) = tgas(i) * (gamma2(i) - 1.d0)/(oldgamma2(i) - 1.d0)
-             oldgamma2(i) = gamma2(i)
-             if(tgas(i).ne.tgas(i).or.tgas(i).lt.0.0d0)then
-               write(0,*)'tgas issue', tgas(i), gamma2(i), x, nH2,
-     &                      nother, gamma
-               stop
-             endif
-           enddo
-#endif
-c         write(0,*) tgas(i)
-         endif
-         enddo
-      endif
-c      write(0,*) tgas(8)
-c
-c     If this is the first time through, just set tgasold to tgas
-c
-      if (iter .eq. 1) then
-         do i = is+1, ie+1
-            if (itmask(i)) then 
-c               write(6,*) i, tgas(i)
-               tgasold(i) = tgas(i)
-            endif
-         enddo
-      endif
-c
-c     --- 6 species cooling ---
-c
-      do i = is+1, ie+1
-
-         if (itmask(i)) then
-
-c
-c        Compute log temperature and truncate if above/below table max/min
-c
-c         logtem(i) = log(0.5*(tgas(i)+tgasold(i)))
-c Changed to be non-time-centered
-         logtem(i) = log(tgas(i))
-         logtem(i) = max(logtem(i), logtem0)
-         logtem(i) = min(logtem(i), logtem9)
-c
-c        Compute index into the table and precompute parts of linear interp
-c
-         indixe(i) = min(nratec-1,
-     .                  max(1,int((logtem(i)-logtem0)/dlogtem)+1))
-         t1(i) = (logtem0 + (indixe(i) - 1)*dlogtem)
-         t2(i) = (logtem0 + (indixe(i)    )*dlogtem)
-         tdef(i) = t2(i) - t1(i)
-         tdef_inv(i) = 1.d0 / tdef(i)
-c
-c        Lookup cooling values and do a linear temperature in log(T)
-c
-         turnoff = 1.d0
-         if(d(i,j,k)*dom.gt.1e18)turnoff=0.d0
-         ceHI(i) = turnoff*(ceHIa(indixe(i)) + (logtem(i) - t1(i))
-     .         *(ceHIa(indixe(i)+1) -ceHIa(indixe(i)))*tdef_inv(i))
-         ceHeI(i) = turnoff*(ceHeIa(indixe(i)) + (logtem(i) - t1(i))
-     .         *(ceHeIa(indixe(i)+1) -ceHeIa(indixe(i)))*tdef_inv(i))
-         ceHeII(i) = turnoff*(ceHeIIa(indixe(i)) + (logtem(i) - t1(i))
-     .         *(ceHeIIa(indixe(i)+1) -ceHeIIa(indixe(i)))*tdef_inv(i))
-         ciHI(i) = (ciHIa(indixe(i)) + (logtem(i) - t1(i))
-     .         *(ciHIa(indixe(i)+1) -ciHIa(indixe(i)))*tdef_inv(i))!*turnoff
-         ciHeI(i) = (ciHeIa(indixe(i)) + (logtem(i) - t1(i))
-     .         *(ciHeIa(indixe(i)+1) -ciHeIa(indixe(i)))*tdef_inv(i))!*turnoff
-         ciHeIS(i) = (ciHeISa(indixe(i)) + (logtem(i) - t1(i))
-     .         *(ciHeISa(indixe(i)+1) -ciHeISa(indixe(i)))*tdef_inv(i))!*turnoff
-         ciHeII(i) = (ciHeIIa(indixe(i)) + (logtem(i) - t1(i))
-     .         *(ciHeIIa(indixe(i)+1) -ciHeIIa(indixe(i)))*tdef_inv(i))!*turnoff
-         reHII(i) = (reHIIa(indixe(i)) + (logtem(i) - t1(i))
-     .         *(reHIIa(indixe(i)+1) -reHIIa(indixe(i)))*tdef_inv(i))!*turnoff
-         reHeII1(i)=(reHeII1a(indixe(i)) + (logtem(i) - t1(i))
-     .        *(reHeII1a(indixe(i)+1)-reHeII1a(indixe(i)))*tdef_inv(i))!*turnoff
-         reHeII2(i)=(reHeII2a(indixe(i)) + (logtem(i) - t1(i))
-     .        *(reHeII2a(indixe(i)+1)-reHeII2a(indixe(i)))*tdef_inv(i))!*turnoff
-         reHeIII(i)=(reHeIIIa(indixe(i)) + (logtem(i) - t1(i))
-     .        *(reHeIIIa(indixe(i)+1)-reHeIIIa(indixe(i)))*tdef_inv(i))!*turnoff
-         brem(i) = turnoff*(brema(indixe(i)) + (logtem(i) - t1(i))
-     .         *(brema(indixe(i)+1) -brema(indixe(i)))*tdef_inv(i))
-         cieco(i) = ciecoa(indixe(i)) + (logtem(i) - t1(i))
-     .      *(ciecoa(indixe(i)+1) - ciecoa(indixe(i)))*tdef_inv(i)
-      endif
-      enddo
-c
-c     Compute the cooling function
-c
-      do i = is+1, ie+1
-         if (itmask(i)) then
-         
-!         if(omask(i))write(0,*)i,j,k,'a0',edotminus(i),edotplus(i),tgas(i)
-         edotminus(i) = edotminus(i) + abs(
-c
-c                    Collisional excitations
-c
-     .             - (ceHI  (i)*HI  (i))*de(i)           ! ce of HI
-     .             - ((ceHeI (i)*HeII(i))*de(i))
-     .                  *de(i)*dom/4.d0                       ! ce of HeI
-     .             - (ceHeII(i)*HeII(i))*de(i)/4.d0         ! ce of HeII
-     .                 ) 
-!         if(omask(i))write(0,*)i,j,k,'a1',edotminus(i),edotplus(i),tgas(i)
-c
-c                    Collisional ionizations
-c
-         edotminus(i) = edotminus(i) + abs(
-     .             - (ciHI  (i)*HI  (i))*de(i)             ! ci of HI
-     .             - (ciHeI (i)*HeI (i))*de(i)/4.d0         ! ci of HeI
-     .             - (ciHeII(i)*HeII(i))*de(i)/4.d0         ! ci of HeII
-     .             - ((ciHeIS(i)*HeII(i))*de(i))
-     .                  *de(i)*dom/4.d0  ! ci of HeIS
-     .                 ) 
-!         if(omask(i))write(0,*)i,j,k,'a2',edotminus(i),edotplus(i),tgas(i)
-c
-c                    Recombinations
-c
-         edotminus(i) = edotminus(i) + abs(
-     .             - (reHII  (i)*HII  (i))*de(i)           ! re of HII
-     .             - (reHeII1(i)*HeII (i))*de(i)/4.d0      ! re of HeII
-     .             - (reHeII2(i)*HeII (i))*de(i)/4.d0      ! re of HeII
-     .             - (reHeIII(i)*HeIII(i))*de(i)/4.d0      ! re of HeIII
-     .                 ) 
-!         if(omask(i))write(0,*)i,j,k,'a3',edotminus(i),edotplus(i),tgas(i)
-!         if(omask(i))write(0,*)i,j,k,'re1',reHII(i),HII(i),de(i)
-!         if(omask(i))write(0,*)i,j,k,'re2',reHeII1(i),HeII(i),de(i)/4.d0
-!         if(omask(i))write(0,*)i,j,k,'re3',reHeII2(i),HeII(i),de(i)/4.d0
-!         if(omask(i))write(0,*)i,j,k,'re4',reHeIII(i),HeIII(i),de(i)/4.d0
-c
-c                    Bremsstrahlung
-c
-         edotminus(i) = edotminus(i) + abs(
-     .             - (brem(i)*(HII(i)+HeII(i)/4.d0+HeIII(i)))
-     .                        *de(i)
-     .                 ) 
-!         if(omask(i))write(0,*)i,j,k,'a4',edotminus(i),edotplus(i),tgas(i)
-          
-c
-c                    Compton cooling or heating
-c
-           if(omask(i))write(0,*)'stats',i,j,k,d(i,j,k),
-     &                 HI(i),HII(i),de(i),HeI(i),HeII(i),HeIII(i)
-           comp_edot = - comp1*(tgas(i)-comp2)*de(i)/dom
-           if (comp_edot.lt.0.0d0) then
-!         if(omask(i))write(0,*)i,j,k,'b',edotminus(i),edotplus(i),tgas(i)
-             edotminus(i) = edotminus(i)
-     .          + abs(comp_edot)
-           else
-!         if(omask(i))write(0,*)i,j,k,'c',edotminus(i),edotplus(i),tgas(i)
-             edotplus(i) = edotplus(i)
-     .          + abs(comp_edot)
-           endif
-c
-c                    X-ray compton heating
-c
-           comp_xray_edot=-comp_xraya*(tgas(i)-comp_temp)*de(i)/dom
-           if (comp_xray_edot.lt.0.0d0) then
-!         if(omask(i))write(0,*)i,j,k,'d1',edotminus(i),edotplus(i),tgas(i)
-             edotminus(i) = edotminus(i) 
-     .             + abs(comp_xray_edot)
-           else
-!         if(omask(i))write(0,*)i,j,k,'d2',edotminus(i),edotplus(i),tgas(i)
-             edotplus(i) = edotplus(i)
-     .            +  abs(comp_xray_edot)
-           endif
-         endif
-      enddo
-c     
-c     --- H2 cooling ---
-c
-      if (ispecies .gt. 1) then
-c
-
-#define USE_GLOVER_ABEL2008
-#ifdef USE_GLOVER_ABEL2008
-         do i = is+1, ie+1
-            if ( itmask(i) ) then
-            gaHI(i) = gaHIa(indixe(i)) + (logtem(i) - t1(i))
-     &         *(gaHIa(indixe(i)+1) - gaHIa(indixe(i)))*tdef_inv(i)
-            gaH2(i) = gaH2a(indixe(i)) + (logtem(i) - t1(i))
-     &         *(gaH2a(indixe(i)+1) - gaH2a(indixe(i)))*tdef_inv(i)
-            gaHe(i) = gaHea(indixe(i)) + (logtem(i) - t1(i))
-     &         *(gaHea(indixe(i)+1) - gaHea(indixe(i)))*tdef_inv(i)
-            gaHp(i) = gaHpa(indixe(i)) + (logtem(i) - t1(i))
-     &         *(gaHpa(indixe(i)+1) - gaHpa(indixe(i)))*tdef_inv(i)
-            gael(i) = gaela(indixe(i)) + (logtem(i) - t1(i))
-     &         *(gaela(indixe(i)+1) - gaHIa(indixe(i)))*tdef_inv(i)
-            gphdl(i) = gphdla(indixe(i)) + (logtem(i) - t1(i))
-     &         *(gphdla(indixe(i)+1) - gphdla(indixe(i)))*tdef_inv(i)
-            end if
-         enddo
-
-         do i = is+1, ie+1
-            if ( itmask(i) ) then
-            if (ih2optical.eq.1) then
-c
-c            nH2 = 0.5*H2I(i)
-c            nother = (HeI(i) + HeII(i) + HeIII(i))/4.0 +
-c     .                HI(i) + HII(i) + de(i)
-c            fH2 = nH2/(nH2 + nother)
-c            fudge = sqrt((40.d0 * 10**(4.8d0 * 
-c     .               sqrt(max(log10(tgas(i)),2.d0)-2.d0)) / fH2**2)/
-c     .               ((nH2 + nother)*dom) )
-              fudge = (0.76d0*d(i,j,k)*dom/8.0d9)**(-0.45d0)
-              fudge = min(fudge, 1.d0)
-            endif
-
-            galdl(i) = gaHI(i) * HI(i)  + gaH2(i) * H2I(i)
-     &               + gaHe(i) * HeI(i) + gaHp(i) * HII(i)
-     &               + gael(i) * de(i)
-            gphdl1 = gphdl(i)/dom
-            edotminus(i) = edotminus(i) + float(ih2co)*fudge*H2I(i)*
-     &           gphdl(i)/(1.d0 + gphdl1/galdl(i)) / (2.d0*dom)
-
-            end if
-         enddo
-#else
-
-#define USE_GALLI_PALLA1999
-
-c        Use the Galli and Palla (1999) cooling rates for molecular H.
-c
-#ifdef USE_GALLI_PALLA1999
-c
-c
-         do i = is+1, ie+1
-            if ( itmask(i) ) then
-            gpldl(i) = gpldla(indixe(i)) + (logtem(i) - t1(i))
-     &         *(gpldla(indixe(i)+1) - gpldla(indixe(i)))*tdef_inv(i)
-            gphdl(i) = gphdla(indixe(i)) + (logtem(i) - t1(i))
-     &         *(gphdla(indixe(i)+1) - gphdla(indixe(i)))*tdef_inv(i)
-            cieco(i) = ciecoa(indixe(i)) + (logtem(i) - t1(i))
-     .         *(ciecoa(indixe(i)+1) - ciecoa(indixe(i)))*tdef_inv(i)
-            endif
-         enddo
-
-         do i = is+1, ie+1
-
-            if (itmask(i)) then
-            fudge = 1.0d0
-            if (ih2optical.eq.1) then
-c
-c            nH2 = 0.5*H2I(i)
-c            nother = (HeI(i) + HeII(i) + HeIII(i))/4.0 +
-c     .                HI(i) + HII(i) + de(i)
-c            fH2 = nH2/(nH2 + nother)
-c            fudge = sqrt((40.0 * 10**(4.8 * 
-c     .               sqrt(max(log10(tgas(i)),2.0)-2.0)) / fH2**2)/
-c     .               ((nH2 + nother)*dom) )
-              fudge = (0.76d0*d(i,j,k)*dom/8.0d9)**(-0.45d0)
-              fudge = min(fudge, 1.d0)
-              ciefudge = 1.d0
-              tau = ((d(i,j,k)/2d16)*dom)**2.8d0  ! 2e16 is in units of cm^-3
-              tau = max(tau, 1.d-5)
-              ciefudge = min((1.d0-exp(-tau))/tau,1.d0)
-              fudge = fudge * ciefudge
-c              ciefudge = 1.d0
-c              tau = ((d(i,j,k)/2.d16)*dom)**2.8d0  ! 2e16 is in units of cm^-3
-c              tau = max(tau, 1.d-5)
-c              ciefudge = min((1.d0-exp(-tau))/tau,1.d0)
-c              fudge = fudge! * ciefudge
-c              if (fudge.lt.1.d-2)
-c     .          write(6,*) 'fudge:', fudge
-            endif
-c
-            gphdl1 = gphdl(i)/(HI(i)*dom)
-!         if(omask(i))write(0,*)i,j,k,'e',edotminus(i),edotplus(i),tgas(i)
-            edotminus(i) = edotminus(i) + abs(
-     .           float(ih2co)*fudge*H2I(i)*
-     .                 gphdl(i)/(1.d0 + gphdl1/gpldl(i)) / (2.d0*dom) )
-!         if(omask(i))write(0,*)i,j,k,'e2',edotminus(i),edotplus(i),tgas(i),gphdl(i),gphdl1,gpldl(i),H2I(i)
-            endif
-         enddo
-c
-#else /* USE_GALLI_PALLA1999 */
-cc
-         do i = is+1, ie+1
-            if (itmask(i)) then
-            hyd01k(i) = hyd01ka(indixe(i)) + (logtem(i) - t1(i))
-     .         *(hyd01ka(indixe(i)+1)-hyd01ka(indixe(i)))*tdef_inv(i)
-            h2k01(i) = h2k01a(indixe(i)) + (logtem(i) - t1(i))
-     .         *(h2k01a(indixe(i)+1) - h2k01a(indixe(i)))*tdef_inv(i)
-            vibh(i) = vibha(indixe(i)) + (logtem(i) - t1(i))
-     .         *(vibha(indixe(i)+1) - vibha(indixe(i)))*tdef_inv(i)
-            roth(i) = rotha(indixe(i)) + (logtem(i) - t1(i))
-     .         *(rotha(indixe(i)+1) - rotha(indixe(i)))*tdef_inv(i)
-            rotl(i) = rotla(indixe(i)) + (logtem(i) - t1(i))
-     .         *(rotla(indixe(i)+1) - rotla(indixe(i)))*tdef_inv(i)
-            cieco(i) = ciecoa(indixe(i)) + (logtem(i) - t1(i))
-     .         *(ciecoa(indixe(i)+1) - ciecoa(indixe(i)))*tdef_inv(i)
-            endif
-         enddo
-c
-         do i = is+1, ie+1
-            if (itmask(i)) then
-            qq   = 1.2d0*(HI(i)*dom)**0.77d0 + 
-     .                (H2I(i)*dom/2.d0)**0.77d0
-            vibl = (HI(i)*hyd01k(i) + 
-     .             H2I(i)/2.*h2k01(i))
-     .             *dom*8.18d-13
-c
-            if (ih2optical .eq. 1) then
-c            nH2 = 0.5d0*H2I(i)
-c            nother = (HeI(i) + HeII(i) + HeIII(i))/4.d0 +
-c     .                HI(i) + HII(i) + de(i)
-c            fH2 = nH2/(nH2 + nother)
-c            fudge = sqrt((40.d0 * 10.d0**(4.8d0 * 
-c     .               sqrt(max(log10(tgas(i)),2.d0)-2.d0)) / fH2**2)/
-c     .               ((nH2 + nother)*dom) )
-              fudge = (0.76d0*d(i,j,k)*dom/8.0d9)**(-0.45d0)
-              fudge = min(fudge, 1.d0)
-            endif
-c
-!         if(omask(i))write(0,*)i,j,k,'f',edotminus(i),edotplus(i),tgas(i)
-            edotminus(i)=edotminus(i) + float(ih2co)*fudge*H2I(i)*(
-     .          vibh(i)/(1.d0+vibh(i)/max(   vibl     ,tiny)) +
-     .          roth(i)/(1.d0+roth(i)/max(qq*rotl(i),tiny))     
-     .                                                      )/2.d0/dom
-            endif
-         enddo
-c
-#endif /* USE_GALLI_PALLA1999 */
-#endif
-c
-
-c     CIE
-c     cooling from H2-H2 and He-H2 collisional induced emission comes
-C     with its own radiative transfer correction as discussed in
-C     Ripamonti & Abel 2003
-         if (iciecool.eq.1) then
-            do i = is+1, ie+1
-            if (itmask(i)) then
-c     Only calculate if H2I(i) is a substantial fraction
-              if (d(i,j,k)*dom.gt.1.d10) then
-                ciefudge = 1.d0
-                tau = ((d(i,j,k)/2e16)*dom)**2.8d0  ! 2e16 is in units of cm^-3
-                tau = max(tau, 1.d-5)
-                ciefudge = min((1.d0-exp(-tau))/tau,1.d0)
-c               Matt's attempt at a second exponentialier cutoff
-                tau = ((d(i,j,k)/2.d18)*dom)**8.d0  ! 2e16 is in units of cm^-3
-                tau = max(tau, 1.d-5)
-                ciefudge = ciefudge*min((1.-exp(-tau))/tau,1.d0)
-c We apply the ciefudge factor at the end, because it should apply to the
-c entire continuum
-                edotminus(i) = ciefudge*edotminus(i) 
-     &                 + ciefudge*H2I(i)*(d(i,j,k)*cieco(i))
-              endif
-            endif
-            enddo
-         endif
-
-      endif                     !  ispecies = 2
-c
-c     --- Cooling from HD ---
-c
-      if (ispecies .gt. 2) then
-         do i = is+1, ie+1
-            if (itmask(i)) then
-c            hdlte(i) = hdltea(indixe(i)) + (logtem(i) - t1(i))
-c     .         *(hdltea(indixe(i)+1) - hdltea(indixe(i)))*tdef_inv(i)
-c            hdlow(i) = hdlowa(indixe(i)) + (logtem(i) - t1(i))
-c     .         *(hdlowa(indixe(i)+1) - hdlowa(indixe(i)))*tdef_inv(i)
-#ifdef PYFORTFIEJ
-            ttgas = tgas(i)
-            call calc_hdcool(ttgas, temp_table)
-            do m = 1, 5
-              hdcool(i,m) = temp_table(m)
-            enddo
-#else
-            do m = 1, 5
-              hdcool(i,m) = hdcoola(indixe(i),m) + (logtem(i) - t1(i))
-     .        *(hdcoola(indixe(i)+1,m)-hdcoola(indixe(i),m))*tdef_inv(i)
-            enddo
-#endif
-            endif
-         enddo
-c
-         do i = is+1, ie+1
-            if (itmask(i)) then
-c               hdlte1 = hdlte(i)/(HDI(i)*dom/2.0)
-c               hdlow1 = max(hdlow(i), tiny)
-c               edotminus(i) = edotminus(i) + HDI(i)*
-c     .                     (hdlte1/(1.0 + hdlte1/hdlow1)/(2.0*dom))
-c   Do we want to consider HII as well?
-c   Note that, as per Flower 2000, we do not include H2 collisions
-                hdcoolval(i) = 0.d0
-                do m = 1,5
-                  hdcoolval(i) = hdcoolval(i) + (hdcool(i,m)
-     .                          *(log10(HI(i))*dom)**(m-1))
-                enddo
-                ciefudge = 1.
-                tau = ((d(i,j,k)/2.d16)*dom)**2.8d0  ! 2e16 is in units of cm^-3
-                tau = max(tau, 1.d-5)
-                ciefudge = min((1.d0-exp(-tau))/tau,1.d0)
-c                omask(i) = .true.
-!         if(omask(i))write(0,*)i,'h',edotminus(i),edotplus(i),tgas(i),
-!     .                               (HDI(i)*10**hdcoolval(i)/coolunit)*ciefudge
-c                omask(i) = .false.
-                edotminus(i) = edotminus(i)
-     .                   + 0*HDI(i)*(10**hdcoolval(i)/coolunit)*ciefudge
-            endif
-            omask(i) = .false.
-         enddo
-      endif
-c
-c     --- Compute (external) radiative heating terms ---
-c
-c                       Photoionization heating
-c
-#ifdef RADIATION
-      if (iradshield .eq. 0) then
-c
-c        regular version
-c
-         if (iradtype .eq. 8) then
-c
-c           1) heating assuming high energy photons produces secondary
-c              electrons which do the heating (Shull & Steenberg, 1985).
-c
-            do i = is+1, ie+1
-               if (itmask(i)) then
-               x = max(HII(i)/(HI(i)+HII(i)), 1.d-4)
-               factor = 0.9971d0*(1.d0-(1.d0-x**0.2663d0)**1.3163d0)
-               edotplus(i) = edotplus(i) + float(ipiht)*factor*(
-     .                + piHI  *HI  (i)            ! pi of HI
-     .                + piHeI *HeI (i)*0.25d0     ! pi of HeI
-     .                + piHeII*HeII(i)*0.25d0     ! pi of HeII
-     .              )/dom
-               endif
-            enddo
-         else
-c
-c           2) standard heating
-c
-            do i = is+1, ie+1
-               if (itmask(i)) then
-               edotplus(i) = edotplus(i) + float(ipiht)*(
-     .                + piHI  *HI  (i)            ! pi of HI
-     .                + piHeI *HeI (i)*0.25d0     ! pi of HeI
-     .                + piHeII*HeII(i)*0.25d0     ! pi of HeII
-     .              )/dom
-               endif
-            enddo
-         endif
-      else
-c
-c        version with approximate self-shielding
-c
-         do i = is+1, ie+1
-            if (itmask(i)) then
-            edotplus(i) = edotplus(i) + float(ipiht)*(
-     .                + piHI  *HI  (i)*
-     .                   exp(-avgsighp*HI(i)*dom)
-     .                + piHeI *HeI (i)*0.25d0*
-     .                   exp(-avgsighep*0.25d0*HeI(i)*dom)
-     .                + piHeII*HeII(i)*0.25d0*
-     .                   exp(-avgsighe2p*0.25d0*HeII(i)*dom)
-     .           )/dom
-            endif
-         enddo
-      endif
-c
-#endif /* RADIATION */
-      return
-      end
-

src/enzo/primordial_derivstep.inc

-c
-c           Calculate all of our derivatives
-c
-            call primordial_cool(
-     &                d, e, s(1,ige), u, v, w,
-     &                s(1,ide), s(1,iHI), s(1,iHII),
-     &                s(1,iHeI), s(1,iHeII), s(1,iHeIII),
-     &                in, jn, kn, nratec, idual, imethod,
-     &                iexpand, ispecies, idim,
-     &                is, ie, j, k, ih2co, ipiht, iter,
-     &                aye, temstart, temend,
-     &                utem, uxyz, uaye, urho, utim,
-     &                eta1, eta2, gamma,
-     &                ceHIa, ceHeIa, ceHeIIa, ciHIa, ciHeIa, 
-     &                ciHeISa, ciHeIIa, reHIIa, reHeII1a, 
-     &                reHeII2a, reHeIIIa, brema, compa, 
-     &                comp_xraya, comp_temp,
-     &                piHI, piHeI, piHeII, comp1, comp2,
-     &                s(1,iHM), s(1,iH2I), s(1,iH2II),
-     &                s(1,iDI), s(1,iDII), s(1,iHDI),
-     &                hyd01ka, h2k01a, vibha, rotha, rotla,
-     &                hyd01k, h2k01, vibh, roth, rotl,
-     &                gpldla, gphdla, gpldl, gphdl,
-     &                hdltea, hdlowa, hdlte, hdlow, 
-     &                hdcoola, hdcool, ciecoa, cieco, 
-     &                gaHIa, gaH2a, gaHea, gaHpa, gaela,
-     &                ceHI, ceHeI, ceHeII, ciHI, ciHeI, ciHeIS, ciHeII,
-     &                reHII, reHeII1, reHeII2, reHeIII, brem,
-     &                indixe, t1, t2, logtem, tdef, dsq(1,ige),
-     &                dsp(1,ige), tgas, tgasold, p2d,
-     &                inutot, iradtype, nfreq, 
-     &                iradshield, avgsighp, avgsighep, 
-     &                avgsighe2p,itmask, iciecool,ih2optical,ifedtgas,
-     &                gamma2, omask )
-#ifdef UNUSED_TABULATED_EQ
-            do i=is+1,ie+1
-            if(itmask(i).eqv..true..and.eqmask(i).eqv..true.)then
-              geprime(i) = log10((ge(i,j,k) +
-     &                   (fh - H2I(i,j,k)/d(i,j,k))*chunit)*(uvel**2.0))
-            endif
-            enddo
-#endif
-c
-c        Look-up rates as a function of temperature for 1D set of zones
-c
-            call primordial_lookup_rates(temstart, temend, nratec, j, k,
-     &                is, ie, ijk, iradtype, iradshield, in, jn, kn,
-     &                ispecies, tgas, tgasold, d, 
-     &                s(1,iHI), s(1,iHII), s(1,iHeI), s(1,iHeII),
-     &                k1a, k2a, k3a, k4a, k5a, k6a, k7a, k8a, k9a, k10a,
-     &                k11a, k12a, k13a, k13dda, k14a, k15a, k16a,
-     &                k17a, k18a, k19a, k21a, k22a, k23a,
-     &                k50a, k51a, k52a, k53a, k54a, k55a, k56a,
-     &                avgsighp, avgsighep, avgsighe2p, piHI, piHeI,
-     &                k1, k2, k3, k4, k5, k6, k7, k8, k9, k10,
-     &                k11, k12, k13, k14, k15, k16, k17, k18,
-     &                k19, k21, k22, k23, k24, k25, k26,
-     &                k50, k51, k52, k53, k54, k55,
-     &                k56, k13dd, k24shield, k25shield, k26shield,
-     &                t1, t2, tdef, logtem, indixe, 
-     &                dom, coolunit, tbase1, itmask,
-     &                eqmask, threebody )
-            do i=is+1,ie+1
-              dmask(i) = itmask(i)
-              eulermask(i) = .false.
-            enddo
-
-