diff --git a/src/enzo/CosmologySimulationInitialize.C b/src/enzo/CosmologySimulationInitialize.C
index 95b37ca88..5e6f0ecfc 100644
--- a/src/enzo/CosmologySimulationInitialize.C
+++ b/src/enzo/CosmologySimulationInitialize.C
@@ -111,6 +111,10 @@ int CosmologySimulationInitialize(FILE *fptr, FILE *Outfptr,
   char *HDIName   = "HDI_Density";
   char *MetalName = "Metal_Density";
   char *GPotName  = "Grav_Potential";
+  char *MachName   = "Mach";
+  char *CRName     = "CR_Density";
+  char *PSTempName = "PreShock_Temperature";
+  char *PSDenName  = "PreShock_Density";
   char *ExtraNames[2] = {"Z_Field1", "Z_Field2"};
  
  
@@ -589,7 +593,14 @@ int CosmologySimulationInitialize(FILE *fptr, FILE *Outfptr,
       DataLabel[i++] = HDIName;
     }
   }
- 
+  if (CRModel) {
+    DataLabel[i++] = MachName;
+    if(StorePreShockFields){
+      DataLabel[i++] = PSTempName;
+      DataLabel[i++] = PSDenName;
+    }
+    DataLabel[i++] = CRName;
+  } 
   if (CosmologySimulationUseMetallicityField) {
     DataLabel[i++] = MetalName;
     if(MultiMetals){
diff --git a/src/enzo/EvolveLevel.C b/src/enzo/EvolveLevel.C
index 812a0b0db..9b7219199 100644
--- a/src/enzo/EvolveLevel.C
+++ b/src/enzo/EvolveLevel.C
@@ -34,6 +34,8 @@
 /                Optional StaticSiblingList for root grid
 /  modified8:  April, 2009 by John Wise
 /                Added star particle class and radiative transfer
+/  modified9:  July, 2009 by Sam Skillman
+/                Added shock and cosmic ray analysis
 /
 /  PURPOSE:
 /    This routine is the main grid evolution function.  It assumes that the
@@ -417,6 +419,10 @@ int EvolveLevel(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[],
 
       Grids[grid1]->GridData->StarParticleHandler
 	(Grids[grid1]->NextGridNextLevel, level);
+
+      /* Include shock-finding and cosmic ray acceleration */
+
+      Grids[grid1]->GridData->ShocksHandler();
  
       /* Gravity: clean up AccelerationField. */
 
diff --git a/src/enzo/Grid.h b/src/enzo/Grid.h
index a08643ed6..bcd43ca72 100644
--- a/src/enzo/Grid.h
+++ b/src/enzo/Grid.h
@@ -512,6 +512,10 @@ class grid
 
    int MultiSpeciesHandler();
 
+/* Handle the selection of shock finding algorithm */
+
+   int ShocksHandler();
+
 /* Solve the radiative cooling/heating equations  */
 
    int SolveRadiativeCooling();
@@ -528,6 +532,14 @@ class grid
 
    int RadiationComputeDensities(int level);
 
+// -------------------------------------------------------------------------
+// Functions for shock finding and cosmic ray acceleration
+//
+   int FindShocks();
+   int FindTempSplitShocks();
+   int FindVelShocks();
+   int FindVelSplitShocks();
+
 // -------------------------------------------------------------------------
 // Functions for grid (re)generation.
 //
@@ -1220,6 +1232,10 @@ void SortParticlesByNumber();
 			    int &HMNum, int &H2INum, int &H2IINum,
                             int &DINum, int &DIINum, int &HDINum);
 
+/* Identify shock/cosmic ray fields. */
+  int IdentifyCRSpeciesFields(int &MachNum, int&CRNum, 
+			      int &PSTempNum, int &PSDenNum);
+
 /* Zeus Solver. */
 
   int ZeusSolver(float *gamma, int igamfield, int nhy, 
@@ -1262,6 +1278,9 @@ void SortParticlesByNumber();
 // Functions for Specific problems (usually problem generator functions).
 //
 
+/* Generalized Extra Field Grid Initializer (returns NumberOfBaryonFields) */
+  int InitializeTestProblemGrid(int field_counter);
+
 /* Protostellar Collapse problem: initialize grid (returns SUCCESS or FAIL) */
   int ProtostellarCollapseInitializeGrid(float CoreDensity,
 					 float CoreEnergy,
diff --git a/src/enzo/Grid_CorrectForRefinedFluxes.C b/src/enzo/Grid_CorrectForRefinedFluxes.C
index 8fd30629a..ae5919078 100644
--- a/src/enzo/Grid_CorrectForRefinedFluxes.C
+++ b/src/enzo/Grid_CorrectForRefinedFluxes.C
@@ -235,8 +235,9 @@ int grid::CorrectForRefinedFluxes(fluxes *InitialFluxes,
 	    total number density summed over ionization, etc.) */
  
 	for (field = 0; field < NumberOfBaryonFields; field++)
-	  if ((RadiativeCooling == 0 || (FieldType[field] != TotalEnergy &&
-		 			 FieldType[field] != InternalEnergy))
+ 	  if (FieldTypeNoInterpolate(FieldType[field]) == FALSE &&
+ 	      (RadiativeCooling == 0 || (FieldType[field] != TotalEnergy &&
+					 FieldType[field] != InternalEnergy))
 	      && (FieldType[field] < ElectronDensity)) {
 	  for (k = Start[2]; k <= End[2]; k++)
 	    for (j = Start[1]; j <= End[1]; j++)
diff --git a/src/enzo/Grid_CosmologySimulationInitializeGrid.C b/src/enzo/Grid_CosmologySimulationInitializeGrid.C
index 94b5e6530..eb2e77d58 100644
--- a/src/enzo/Grid_CosmologySimulationInitializeGrid.C
+++ b/src/enzo/Grid_CosmologySimulationInitializeGrid.C
@@ -117,6 +117,8 @@ int grid::CosmologySimulationInitializeGrid(
   int idim, dim, i, j, vel, OneComponentPerFile, ndim, level;
   int DeNum, HINum, HIINum, HeINum, HeIINum, HeIIINum, HMNum, H2INum, H2IINum,
       DINum, DIINum, HDINum, MetalNum;
+
+  int CRNum, MachNum, PSTempNum, PSDenNum;
  
   int ExtraField[2];
  
@@ -250,6 +252,21 @@ int grid::CosmologySimulationInitializeGrid(
     }
     if (WritePotential)
       FieldType[NumberOfBaryonFields++] = GravPotential;
+    if(CRModel){
+      FieldType[MachNum   = NumberOfBaryonFields++] = Mach;
+      if(StorePreShockFields){
+	FieldType[PSTempNum = NumberOfBaryonFields++] = PreShockTemperature;
+	FieldType[PSDenNum = NumberOfBaryonFields++] = PreShockDensity;
+      }
+      FieldType[CRNum     = NumberOfBaryonFields++] = CRDensity;
+    }
+    if (UseMetallicityField) {
+      FieldType[MetalNum = NumberOfBaryonFields++] = Metallicity;
+      if(MultiMetals){
+	FieldType[ExtraField[0] = NumberOfBaryonFields++] = ExtraType0;
+	FieldType[ExtraField[1] = NumberOfBaryonFields++] = ExtraType1;
+      }
+    }
   }
  
   // Set the subgrid static flag
@@ -421,8 +438,20 @@ int grid::CosmologySimulationInitializeGrid(
 	                             BaryonField[H2INum][i];
       }
  
+      //Shock/Cosmic Ray Model
+      if(CRModel && ReadData)
+	for (i = 0; i < size; i++){
+	  BaryonField[MachNum][i] = tiny_number;
+	  BaryonField[CRNum][i] = tiny_number;
+	  if(StorePreShockFields){
+	    BaryonField[PSTempNum][i] = tiny_number;
+	    BaryonField[PSDenNum][i] = tiny_number;
+	  }
+	  
+	}
+      
     }
- 
+  
   // If using metallicity, set the field
  
   if (UseMetallicityField && ReadData)
diff --git a/src/enzo/Grid_InterpolateBoundaryFromParent.C b/src/enzo/Grid_InterpolateBoundaryFromParent.C
index 529958ec3..86efd38cb 100644
--- a/src/enzo/Grid_InterpolateBoundaryFromParent.C
+++ b/src/enzo/Grid_InterpolateBoundaryFromParent.C
@@ -361,7 +361,7 @@ int grid::InterpolateBoundaryFromParent(grid *ParentGrid)
  
       if (FieldType[field] == Density)
 	FieldPointer = TemporaryDensityField;
-      else
+      else if (FieldTypeNoInterpolate(FieldType[field]) == FALSE)
 	FieldPointer = TemporaryField;
  
       /* Copy needed portion of temp field to current grid. */
diff --git a/src/enzo/Grid_InterpolateFieldValues.C b/src/enzo/Grid_InterpolateFieldValues.C
index b9c237134..6a9ac0473 100644
--- a/src/enzo/Grid_InterpolateFieldValues.C
+++ b/src/enzo/Grid_InterpolateFieldValues.C
@@ -327,8 +327,8 @@ int grid::InterpolateFieldValues(grid *ParentGrid)
  
       if (FieldType[field] == Density)
 	FieldPointer = TemporaryDensityField;
-      else
-	FieldPointer = TemporaryField;
+      else if (FieldTypeNoInterpolate(FieldType[field]) == FALSE)
+	  FieldPointer = TemporaryField;
  
       /* Copy needed portion of temp field to current grid. */
  
diff --git a/src/enzo/Grid_ProjectSolutionToParentGrid.C b/src/enzo/Grid_ProjectSolutionToParentGrid.C
index 8f2445089..a02f3d61b 100644
--- a/src/enzo/Grid_ProjectSolutionToParentGrid.C
+++ b/src/enzo/Grid_ProjectSolutionToParentGrid.C
@@ -120,9 +120,10 @@ int grid::ProjectSolutionToParentGrid(grid &ParentGrid)
  
   if (ProcessorNumber == MyProcessorNumber)
     for (field = 0; field < NumberOfBaryonFields; field++)
-      if (FieldTypeIsDensity(FieldType[field]) == FALSE && (
-	  (FieldType[field] < Velocity1 || FieldType[field] > Velocity3)
-          || HydroMethod != Zeus_Hydro      )       )
+      if ((FieldTypeIsDensity(FieldType[field]) == FALSE &&  
+	   ((FieldType[field] < Velocity1 || FieldType[field] > Velocity3)
+	    || HydroMethod != Zeus_Hydro      )       ) &&
+	  (FieldTypeNoInterpolate(FieldType[field]) == FALSE))
       FORTRAN_NAME(mult3d)(BaryonField[DensNum], BaryonField[field],
 			   &Size, &One, &One, &Size, &One, &One,
 			   &Zero, &Zero, &Zero, &Zero, &Zero, &Zero);
@@ -184,8 +185,9 @@ int grid::ProjectSolutionToParentGrid(grid &ParentGrid)
  
 	  for (i = 0; i < Dim[0]; i += skipi) {
 	    i1 = i/Refinement[0];
-	    ParentGrid.BaryonField[field][pindex+i1] +=
-	      BaryonField[field][gindex+i]*weight;
+	    if (FieldTypeNoInterpolate(FieldType[field]) == FALSE)
+	      ParentGrid.BaryonField[field][pindex+i1] +=
+		BaryonField[field][gindex+i]*weight;
 	  }
 	}
       }
@@ -222,9 +224,10 @@ int grid::ProjectSolutionToParentGrid(grid &ParentGrid)
   /* Divide all fields by mass to return to original quantity. */
  
   for (field = 0; field < NumberOfBaryonFields; field++)
-    if (FieldTypeIsDensity(FieldType[field]) == FALSE && (
+    if ((FieldTypeIsDensity(FieldType[field]) == FALSE && (
 	  (FieldType[field] < Velocity1 || FieldType[field] > Velocity3)
-          || HydroMethod != Zeus_Hydro      )       ) {
+          || HydroMethod != Zeus_Hydro      )       ) &&
+	(FieldTypeNoInterpolate(FieldType[field]) == FALSE)){
       if (ProcessorNumber == MyProcessorNumber)
 	FORTRAN_NAME(div3d)(BaryonField[DensNum], BaryonField[field],
 			    &Size, &One, &One, &Size, &One, &One,
diff --git a/src/enzo/Grid_SolveHydroEquations.C b/src/enzo/Grid_SolveHydroEquations.C
index d2776a9e8..fe84806ec 100644
--- a/src/enzo/Grid_SolveHydroEquations.C
+++ b/src/enzo/Grid_SolveHydroEquations.C
@@ -121,8 +121,44 @@ int grid::SolveHydroEquations(int CycleNumber, int NumberOfSubgrids,
       colnum[NumberOfColours++] = SNColourNum;
     }
 
-    /* Determine if Gamma should be a scalar or a field. */
+    /* Add shock/cosmic ray variables as a colour variable. */
 
+    if(CRModel){
+      int MachNum, CRNum, PSTempNum,PSDenNum;
+      
+      if (IdentifyCRSpeciesFields(MachNum,CRNum,PSTempNum,PSDenNum) == FAIL) {
+	ENZO_FAIL("Error in IdentifyCRSpeciesFields.")
+      }
+      
+      colnum[NumberOfColours++] = MachNum;
+      if(StorePreShockFields){
+	colnum[NumberOfColours++] = PSTempNum;
+	colnum[NumberOfColours++] = PSDenNum;
+      }
+      colnum[NumberOfColours++] = CRNum;
+      
+      /*  Zero out Mach number array so that we don't get strange advection */
+      
+      float *mach = BaryonField[MachNum];
+      if(StorePreShockFields){
+	float *pstemp = BaryonField[PSTempNum];
+	float *psden = BaryonField[PSDenNum];
+	for (i = 0; i < size; i++){
+	  pstemp[i] = tiny_number;
+	  psden[i] = tiny_number;
+	}
+      }
+      for (i = 0; i < size; i++)
+	mach[i] = tiny_number;
+      if(CRModel == 2){  //zero out CR to get instantanous injection
+	float *cr = BaryonField[CRNum];
+	for (i = 0; i < size; i++)
+	  cr[i] = tiny_number;
+      }
+    }
+    
+    /* Determine if Gamma should be a scalar or a field. */
+    
     int UseGammaField = FALSE;
     float *GammaField;
     if (HydroMethod == Zeus_Hydro && MultiSpecies > 1) {
diff --git a/src/enzo/Make.config.objects b/src/enzo/Make.config.objects
index 9270ead35..abb1331bf 100644
--- a/src/enzo/Make.config.objects
+++ b/src/enzo/Make.config.objects
@@ -258,6 +258,7 @@ OBJS_CONFIG_LIB = \
         Grid_destructor.o \
         Grid_DetachForcingFromBaryonFields.o \
         Grid_DoubleMachInitializeGrid.o \
+	Grid_FindShocks.o \
         Grid_FastSiblingLocatorAddGrid.o \
         Grid_FastSiblingLocatorFindSiblings.o \
         Grid_FindAllStarParticles.o \
@@ -282,6 +283,7 @@ OBJS_CONFIG_LIB = \
         Grid_Group_ReadGrid.o \
         Grid_Group_WriteGridInterpolate.o \
         Grid_Group_WriteGrid.o \
+	Grid_IdentifyCRSpeciesFields.o \
         Grid_IdentifyNewSubgrids.o \
         Grid_IdentifyNewSubgridsSmall.o \
         Grid_IdentifyPhysicalQuantities.o \
@@ -350,6 +352,7 @@ OBJS_CONFIG_LIB = \
         Grid_SetParticleMassFlaggingField.o \
         Grid_ShearingBoxInitializeGrid.o \
         Grid_ShockTubeInitializeGrid.o \
+	Grid_ShocksHandler.o \
         Grid_SolveForPotential.o \
         Grid_SolveHydroEquations.o \
         Grid_SolveRadiativeCooling.o \
@@ -369,6 +372,7 @@ OBJS_CONFIG_LIB = \
         Grid_TracerParticleCreateParticles.o \
         Grid_TracerParticleOutputData.o \
         Grid_TracerParticleSetVelocity.o \
+	Grid_TestProblemInitializeGrid.o \
         Grid_TransferSubgridParticles.o \
         Grid_TransferSubgridStars.o \
         Grid_TurbulenceSimulationInitialize.o \
@@ -399,6 +403,7 @@ OBJS_CONFIG_LIB = \
         IdentifyNewSubgridsBySignature.o \
         ImplosionInitialize.o \
         InitializeCloudyCooling.o \
+	InitializeCosmicRayData.o \
         InitializeEquilibriumCoolData.o \
         InitializeLocal.o \
         InitializeMovieFile.o \
@@ -558,6 +563,7 @@ OBJS_CONFIG_LIB = \
         TestGravitySphereInitialize.o \
         TestOrbitInitialize.o \
         TracerParticleCreation.o \
+	TestProblemDataInitialize.o \		
         tscint1d.o \
         tscint2d.o \
         tscint3d.o \
diff --git a/src/enzo/Make.mach.darwin b/src/enzo/Make.mach.darwin
index 658f6b54a..c618bbad0 100644
--- a/src/enzo/Make.mach.darwin
+++ b/src/enzo/Make.mach.darwin
@@ -1,6 +1,6 @@
 #=======================================================================
 #
-# FILE:        Make.mach.rpwagner-cable
+# FILE:        Make.mach.darwin
 #
 # DESCRIPTION: Makefile settings for Leopard OSX
 #  This was written to use:
@@ -28,15 +28,16 @@ MACH_FILE  = Make.mach.darwin
 # Install paths (local variables)
 #-----------------------------------------------------------------------
 
-LOCAL_PACKAGES       = /sw
+LOCAL_PACKAGES       = /usr/local
 
-LOCAL_MPI_INSTALL    = /usr/local
-LOCAL_FC_INSTALL     = /sw/lib/gcc4.4
-LOCAL_HDF5_INSTALL   = /opt/local
+LOCAL_MPI_INSTALL    = /usr
+LOCAL_FC_INSTALL     = /usr/local
+LOCAL_HDF5_INSTALL   = $(LOCAL_PACKAGES)
 LOCAL_SZIP_INSTALL   = $(LOCAL_PACKAGES)
 LOCAL_HYPRE_INSTALL  = $(HOME)
 LOCAL_JBPERF_INSTALL = 
 LOCAL_PNG_INSTALL    = $(LOCAL_PACKAGES)
+LOCAL_PYTHON_INSTALL = /Library/Frameworks/Python.framework/Versions/Current/
 
 #-----------------------------------------------------------------------
 # Compiler settings
@@ -51,6 +52,7 @@ MACH_CXX_MPI   = mpic++
 MACH_FC_MPI    = gfortran
 MACH_F90_MPI   = gfortran
 MACH_LD_MPI    = mpic++
+MACH_CUDACOMPILER = /usr/local/cuda/bin/nvcc
 
 # Without MPI
 
@@ -68,7 +70,8 @@ MACH_LD_NOMPI  = g++ # Linker when not using MPI
 # compile HDF5 with --with-default-api-version=v16, or Enzo with
 # -DH5_USE_16_API.
 
-MACH_DEFINES   = -DLINUX -DH5_USE_16_API -DOPTIMIZED_CTP
+MACH_DEFINES   = -DLINUX -DH5_USE_16_API 
+#-DOPTIMIZED_CTP -DSTATIC_SIBLING_LIST
 
 #-----------------------------------------------------------------------
 # Compiler flag settings
@@ -77,7 +80,7 @@ MACH_DEFINES   = -DLINUX -DH5_USE_16_API -DOPTIMIZED_CTP
 MACH_CPPFLAGS = -P -traditional 
 MACH_CFLAGS   = 
 MACH_CXXFLAGS =
-MACH_FFLAGS   = -fno-second-underscore -ffixed-line-length-132
+MACH_FFLAGS   = -fno-second-underscore
 MACH_F90FLAGS = -fno-second-underscore
 MACH_LDFLAGS  = 
 
@@ -97,8 +100,7 @@ MACH_FFLAGS_REAL_64    = -fdefault-real-8
 MACH_OPT_WARN        = -Wall -g
 MACH_OPT_DEBUG       = -g
 MACH_OPT_HIGH        = -O2
-MACH_OPT_AGGRESSIVE  = -O3 -fomit-frame-pointer -fstrict-aliasing \
-		     -momit-leaf-frame-pointer -fno-tree-pre -falign-loops -g
+MACH_OPT_AGGRESSIVE  = -O3
 
 #-----------------------------------------------------------------------
 # Includes
@@ -106,28 +108,33 @@ MACH_OPT_AGGRESSIVE  = -O3 -fomit-frame-pointer -fstrict-aliasing \
 
 LOCAL_INCLUDES_MPI    = -I$(LOCAL_MPI_INSTALL)/include 
 LOCAL_INCLUDES_HDF5   = -I$(LOCAL_HDF5_INSTALL)/include
+LOCAL_INCLUDES_PYTHON = -I$(LOCAL_PYTHON_INSTALL)/include/python2.5/ \
+                        -I$(LOCAL_PYTHON_INSTALL)/lib/python2.5/site-packages/numpy/core/include
 LOCAL_INCLUDES_HYPRE  = 
 LOCAL_INCLUDES_JBPERF = 
 LOCAL_INCLUDES_PAPI   = # PAPI includes
+LOCAL_INCLUDES_CUDA   = -I/Developer/CUDA/common/inc
 
-MACH_INCLUDES         = $(LOCAL_INCLUDES_HDF5)
-
+MACH_INCLUDES         = $(LOCAL_INCLUDES_HDF5) $(LOCAL_INCLUDES_CUDA)
 MACH_INCLUDES_MPI     = $(LOCAL_INCLUDES_MPI)
 MACH_INCLUDES_HYPRE   = $(LOCAL_INCLUDES_HYPRE)
 MACH_INCLUDES_JBPERF  = $(LOCAL_INCLUDES_JBPERF)
-
+MACH_INCLUDES_PYTHON  = $(LOCAL_INCLUDES_PYTHON)
 
 #-----------------------------------------------------------------------
 # Libraries
 #-----------------------------------------------------------------------
 
-LOCAL_LIBS_MACH  = -L$(LOCAL_FC_INSTALL)/lib -lgfortran
-#LOCAL_LIBS_HDF5  =  -L$(LOCAL_HDF5_INSTALL)/lib -lhdf5 -lsz
-LOCAL_LIBS_HDF5  =  -L$(LOCAL_HDF5_INSTALL)/lib -lhdf5 -lz
+LOCAL_LIBS_MACH   = -L$(LOCAL_FC_INSTALL)/lib -lgfortran
+LOCAL_LIBS_HDF5   = -L$(LOCAL_HDF5_INSTALL)/lib -lhdf5 
+LOCAL_LIBS_PYTHON = -framework Python 
+LOCAL_LIBS_CUDA   = -L/usr/local/cuda/lib/ -lcudart
+
+MACH_LIBS        = $(LOCAL_LIBS_HDF5) $(LOCAL_LIBS_MACH) 
 
-MACH_LIBS        = $(LOCAL_LIBS_HDF5) $(LOCAL_LIBS_MACH)
 MACH_LIBS_MPI    = $(LOCAL_LIBS_MPI)
 MACH_LIBS_HYPRE  = $(LOCAL_LIBS_HYPRE)
 MACH_LIBS_JBPERF = $(LOCAL_LIBS_JBPERF)
-MACH_LIBS_PAPI    = $(LOCAL_LIBS_PAPI)
-
+MACH_LIBS_PAPI   = $(LOCAL_LIBS_PAPI)
+MACH_LIBS_PYTHON = $(LOCAL_LIBS_PYTHON)
+MACH_LIBS_CUDA   = $(LOCAL_LIBS_CUDA)
diff --git a/src/enzo/ReadParameterFile.C b/src/enzo/ReadParameterFile.C
index 40b5443aa..f16826144 100644
--- a/src/enzo/ReadParameterFile.C
+++ b/src/enzo/ReadParameterFile.C
@@ -48,6 +48,7 @@ int ReadListOfInts(FILE *fptr, int N, int nums[]);
 int CosmologyReadParameters(FILE *fptr, FLOAT *StopTime, FLOAT *InitTime);
 int ReadUnits(FILE *fptr);
 int InitializeCloudyCooling(FLOAT Time);
+int InitializeCosmicRayData();
 int InitializeRateData(FLOAT Time);
 int InitializeEquilibriumCoolData(FLOAT Time);
 int InitializeRadiationFieldData(FLOAT Time);
@@ -335,6 +336,12 @@ int ReadParameterFile(FILE *fptr, TopGridData &MetaData, float *Initialdt)
     ret += sscanf(line, "MetalCooling = %d", &MetalCooling);
     if (sscanf(line, "MetalCoolingTable = %s", dummy) == 1) 
       MetalCoolingTable = dummy;
+
+    ret += sscanf(line, "CRModel = %"ISYM, &CRModel);
+    ret += sscanf(line, "ShockMethod = %"ISYM, &ShockMethod);
+    ret += sscanf(line, "ShockTemperatureFloor = %"FSYM, &ShockTemperatureFloor);
+    ret += sscanf(line, "StorePreShockFields = %"ISYM, &StorePreShockFields);
+
     ret += sscanf(line, "RadiationFieldType = %"ISYM, &RadiationFieldType);
     ret += sscanf(line, "AdjustUVBackground = %"ISYM, &AdjustUVBackground);
     ret += sscanf(line, "SetUVBAmplitude = %"FSYM, &SetUVBAmplitude);
@@ -817,6 +824,15 @@ int ReadParameterFile(FILE *fptr, TopGridData &MetaData, float *Initialdt)
     }
   }
 
+  /* If set, initialze Cosmic Ray Efficiency Models */
+
+  if(CRModel){
+    if(InitializeCosmicRayData() == FAIL){
+      ENZO_FAIL("Error in Initialize CosmicRayData.");
+      return FAIL;
+    }
+  }
+
   /* If set, initialize CloudyCooling. */
 
   if (MetalCooling == CLOUDY_METAL_COOLING) {
diff --git a/src/enzo/SetDefaultGlobalValues.C b/src/enzo/SetDefaultGlobalValues.C
index 1d4a6f326..12609f630 100644
--- a/src/enzo/SetDefaultGlobalValues.C
+++ b/src/enzo/SetDefaultGlobalValues.C
@@ -249,6 +249,10 @@ int SetDefaultGlobalValues(TopGridData &MetaData)
   RandomForcingMachNumber     = 0.0;               //AK
   RadiativeCooling            = FALSE;             // off
   MultiSpecies                = FALSE;             // off
+  CRModel                     = 0;                 // off
+  ShockMethod                 = 0;                 // temperature unsplit
+  ShockTemperatureFloor       = 1.0;               // Set to 1K
+  StorePreShockFields         = 0;
   RadiationFieldType          = 0;
   RadiationFieldLevelRecompute = 0;
   AdjustUVBackground          = 1;
diff --git a/src/enzo/TestProblemData.h b/src/enzo/TestProblemData.h
index f7bbae272..996070435 100644
--- a/src/enzo/TestProblemData.h
+++ b/src/enzo/TestProblemData.h
@@ -42,4 +42,9 @@ struct TestProblemDataType
   float MaximumTemperature;
   int ResetEnergies;
 
+  /* Shock/Cosmic Ray Fields */
+  int CRModel;
+  int ShockMethod;
+  int StorePreShockFields;
+  float ShockTemperatureFloor;
 };
diff --git a/src/enzo/WriteParameterFile.C b/src/enzo/WriteParameterFile.C
index ee31acb97..45777c89f 100644
--- a/src/enzo/WriteParameterFile.C
+++ b/src/enzo/WriteParameterFile.C
@@ -316,6 +316,10 @@ int WriteParameterFile(FILE *fptr, TopGridData &MetaData)
   fprintf(fptr, "RadiativeTransfer              = %"ISYM"\n", RadiativeTransfer);
   fprintf(fptr, "RadiationXRaySecondaryIon      = %"ISYM"\n", 
 	  RadiationXRaySecondaryIon);
+  fprintf(fptr, "CRModel                        = %"ISYM"\n", CRModel);
+  fprintf(fptr, "ShockMethod                    = %"ISYM"\n", ShockMethod);
+  fprintf(fptr, "ShockTemperatureFloor          = %"FSYM"\n", ShockTemperatureFloor);
+  fprintf(fptr, "StorePreShockFields            = %"ISYM"\n", StorePreShockFields);
   fprintf(fptr, "RadiationFieldType             = %"ISYM"\n", RadiationFieldType);
   fprintf(fptr, "AdjustUVBackground             = %"ISYM"\n", AdjustUVBackground);
   fprintf(fptr, "SetUVBAmplitude                = %"GSYM"\n", SetUVBAmplitude);
@@ -329,6 +333,7 @@ int WriteParameterFile(FILE *fptr, TopGridData &MetaData)
 
   fprintf(fptr, "OutputCoolingTime              = %"ISYM"\n", OutputCoolingTime);
   fprintf(fptr, "OutputTemperature              = %"ISYM"\n", OutputTemperature);
+
   fprintf(fptr, "OutputSmoothedDarkMatter       = %"ISYM"\n", 
 	  OutputSmoothedDarkMatter);
   fprintf(fptr, "SmoothedDarkMatterNeighbors    = %"ISYM"\n", 
diff --git a/src/enzo/euler.src b/src/enzo/euler.src
index 3b1a0fae2..7002aa00f 100644
--- a/src/enzo/euler.src
+++ b/src/enzo/euler.src
@@ -526,7 +526,7 @@ c
 c             colslice(i,j,n) = (colslice(i,j,n)*dslice(i,j) +
 c     &         (colf(i,j,n)-colf(i+1,j,n))/dx(i)      )/dnu(i) ! c*d conserved
              if (colslice(i,j,n)              +
-     &         (colf(i,j,n)-colf(i+1,j,n))/dx(i) .le. 0.0) write(6,*)
+     &         (colf(i,j,n)-colf(i+1,j,n))/dx(i) .lt. 0.0) write(6,*)
      &         'euler_c:',i,j,n,dx(i),
      &          colf(i-1,j,n),colf(i,j,n),colf(i+1,j,n),
      &          colslice(i-1,j,n),colslice(i,j,n),colslice(i+1,j,n),
diff --git a/src/enzo/global_data.h b/src/enzo/global_data.h
index 5a31672c3..6ce43c31c 100644
--- a/src/enzo/global_data.h
+++ b/src/enzo/global_data.h
@@ -258,6 +258,24 @@ EXTERN RateDataType RateData;
 
 EXTERN int MultiMetals;
 
+/* Cosmic Ray Model 
+ * 0: Off - default
+ * 1: On, Let CRs accululate on Grid
+ * 2: On, Zero out CRs each step to only look at instantaneous acceleration
+ * 3: Highly experimental, takes energy out of gas.  Unstable.
+ */
+EXTERN int CRModel; 
+/* Shock Finding Method: Always on when CRModel nonzero
+ * 0: temperature unsplit - default
+ * 1: temperature split 
+ * 2: velocity unsplit
+ * 3: velocity split
+ */
+EXTERN int ShockMethod; 
+EXTERN CosmicRayDataType CosmicRayData;
+EXTERN float ShockTemperatureFloor;
+EXTERN int StorePreShockFields;
+
 /* Type of radiation field. 
    0 - none,                    1 - Haardt & Madau alpha=-1.5
    2 - H&M alpha = -1.8       
diff --git a/src/enzo/typedefs.h b/src/enzo/typedefs.h
index ab2156c9c..2efbb5085 100644
--- a/src/enzo/typedefs.h
+++ b/src/enzo/typedefs.h
@@ -17,6 +17,7 @@
 #include "RateData.h"
 #include "RadiationFieldData.h"
 #include "TestProblemData.h"
+#include "CosmicRayData.h"
 
 /* These are the different types of baryon fields. */
 
@@ -108,7 +109,13 @@ const field_type
   AccelerationField2         = 59, 
   AccelerationField3         = 60,
 
-  FieldUndefined  = 61;
+/* these are required for Sam Skillman's Shock/Cosmic ray models. */
+  Mach            = 61,
+  PreShockTemperature = 62,
+  PreShockDensity = 63,  
+  CRDensity       = 64,
+
+  FieldUndefined  = 65;
    
 /*
 enum field_type {Density, TotalEnergy, InternalEnergy, Pressure,
@@ -122,7 +129,7 @@ enum field_type {Density, TotalEnergy, InternalEnergy, Pressure,
 */
 
 #define FieldTypeIsDensity(A) (((A) >= TotalEnergy && (A) <= Velocity3) ? FALSE : TRUE)
-
+#define FieldTypeNoInterpolate(A) ((((A) >= Mach) && ((A) <= Mach + 1 + CRModel)) ? TRUE : FALSE) 
 /* These are the different types of fluid boundary conditions. */
 
 const boundary_type