diff --git a/src/enzo/EvolveHierarchy.C b/src/enzo/EvolveHierarchy.C
index 2adee9958..f360a8ccc 100644
--- a/src/enzo/EvolveHierarchy.C
+++ b/src/enzo/EvolveHierarchy.C
@@ -454,7 +454,7 @@ int EvolveHierarchy(HierarchyEntry &TopGrid, TopGridData &MetaData,
 #endif
 
     if (MyProcessorNumber == ROOT_PROCESSOR) {
-      printf("TopGrid dt = %"ESYM"     time = %"GOUTSYM"    cycle = %"ISYM,
+      fprintf(stderr,"TopGrid dt = %"ESYM"     time = %"GOUTSYM"    cycle = %"ISYM,
 	     dt, MetaData.Time, MetaData.CycleNumber);
 
       if (ComovingCoordinates) {
diff --git a/src/enzo/Grid.h b/src/enzo/Grid.h
index 20716bce1..beb8115d8 100644
--- a/src/enzo/Grid.h
+++ b/src/enzo/Grid.h
@@ -2148,6 +2148,7 @@ int zEulerSweep(int j, int NumberOfSubgrids, fluxes *SubgridFluxes[],
 
  int TurbulenceSimulationInitializeGrid(TURBULENCE_INIT_PARAMETERS_DECL);
  int ComputeRandomForcingFields(int mode);
+ int ExtraFunction(char * message); //for debugging.
 
   // The following are private since they should only be called by
   // TurbulenceSimulationInitializeGrid()
diff --git a/src/enzo/Grid_ExtraFunction.C b/src/enzo/Grid_ExtraFunction.C
new file mode 100644
index 000000000..cb0a8e96a
--- /dev/null
+++ b/src/enzo/Grid_ExtraFunction.C
@@ -0,0 +1,21 @@
+//
+// ExtraFunction
+// Generic grid member function for debugging.
+//
+
+#include <stdio.h>
+#include "ErrorExceptions.h"
+#include "performance.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"
+
+int grid::ExtraFunction(char * message){
+
+    return SUCCESS;
+}
diff --git a/src/enzo/Grid_SolveRadiativeCooling.C b/src/enzo/Grid_SolveRadiativeCooling.C
index 8e40b4d95..452785dbe 100644
--- a/src/enzo/Grid_SolveRadiativeCooling.C
+++ b/src/enzo/Grid_SolveRadiativeCooling.C
@@ -86,7 +86,7 @@ extern "C" void FORTRAN_NAME(multi_cool)(
 extern "C" void FORTRAN_NAME(solve_cool)(
 	float *d, float *e, float *ge, float *u, float *v, float *w,
 	int *in, int *jn, int *kn, int *nratec, int *iexpand,
-	hydro_method *imethod, int *idual, int *idim, int *igammah,
+	hydro_method *imethod, int *cool_method, int *idual, int *idim, int *igammah,
 	int *is, int *js, int *ks, int *ie, int *je, int *ke,
 	float *dt, float *aye, float *temstart, float *temend,
 	   float *fh,
@@ -338,7 +338,7 @@ int grid::SolveRadiativeCooling()
        density, totalenergy, gasenergy, velocity1, velocity2, velocity3,
        GridDimension, GridDimension+1, GridDimension+2,
           &CoolData.NumberOfTemperatureBins, &ComovingCoordinates,
-          &HydroMethod,
+          &HydroMethod,&RadiativeCoolingModel,
        &DualEnergyFormalism, &GridRank, &PhotoelectricHeating,
        GridStartIndex, GridStartIndex+1, GridStartIndex+2,
           GridEndIndex, GridEndIndex+1, GridEndIndex+2,
diff --git a/src/enzo/Group_WriteAllData.C b/src/enzo/Group_WriteAllData.C
index 4fef0c88f..5e74b5143 100644
--- a/src/enzo/Group_WriteAllData.C
+++ b/src/enzo/Group_WriteAllData.C
@@ -858,7 +858,7 @@ int Group_WriteAllData(char *basename, int filenumber,
 
   if ( MyProcessorNumber == ROOT_PROCESSOR ){
     sptr = fopen("OutputLog", "a");
-    fprintf(sptr, "DATASET WRITTEN %s \n", name);
+    fprintf(sptr, "DATASET WRITTEN %s %8"ISYM" %18.16e\n", name, MetaData.CycleNumber, MetaData.Time);
     fclose(sptr);
   }
  
diff --git a/src/enzo/InitializeEquilibriumCoolData.C b/src/enzo/InitializeEquilibriumCoolData.C
index 073ed08e7..094b9a8da 100644
--- a/src/enzo/InitializeEquilibriumCoolData.C
+++ b/src/enzo/InitializeEquilibriumCoolData.C
@@ -41,6 +41,8 @@ int InitializeEquilibriumCoolData(FLOAT Time)
  
   /* Open input file for data. */
  
+  if( RadiativeCoolingModel == 1){
+      
   FILE *fptr = fopen("cool_rates.in", "r");
   if (fptr == NULL) {
     ENZO_FAIL("Error opening cool_rates.in\n");
@@ -136,6 +138,12 @@ int InitializeEquilibriumCoolData(FLOAT Time)
 	   float(CoolData.NumberOfTemperatureBins-1),
 	   CoolData.EquilibriumRate[index]);
   fclose(fptr);
+  }else if( RadiativeCoolingModel == 3){
+      CoolData.NumberOfTemperatureBins = 2;
+      CoolData.EquilibriumRate = new float[CoolData.NumberOfTemperatureBins];
+      CoolData.TemperatureEnd = 1e4  ;
+      CoolData.TemperatureStart = 18;
+  }
  
   return SUCCESS;
 }
diff --git a/src/enzo/Make.config.objects b/src/enzo/Make.config.objects
index c21270941..1d059ff96 100644
--- a/src/enzo/Make.config.objects
+++ b/src/enzo/Make.config.objects
@@ -96,6 +96,7 @@ OBJS_CONFIG_LIB = \
         ComputeTable.o \
         Continue.o \
 	cool1d_cloudy.o \
+    cool1d_koyama.o\
         cool1d_multi.o \
         cool1d_sep.o \
         cool1d.o \
@@ -900,6 +901,7 @@ POBJS_CONFIG_LIB = \
         Grid_DeletePhotonPackages.o \
         Grid_DetectIonizationFrontApprox.o \
         Grid_ElectronFractionEstimate.o \
+        Grid_ExtraFunction.o \
         Grid_FinalizeRadiationFields.o \
 	Grid_FindPhotonNewGrid.o \
         Grid_FlagCellsToBeRefinedByOpticalDepth.o \
diff --git a/src/enzo/ReadParameterFile.C b/src/enzo/ReadParameterFile.C
index 2cd9ca8b0..3c0389e22 100644
--- a/src/enzo/ReadParameterFile.C
+++ b/src/enzo/ReadParameterFile.C
@@ -399,6 +399,7 @@ int ReadParameterFile(FILE *fptr, TopGridData &MetaData, float *Initialdt)
     ret += sscanf(line, "RandomForcingMachNumber = %"FSYM, //AK
                   &RandomForcingMachNumber);
     ret += sscanf(line, "RadiativeCooling = %"ISYM, &RadiativeCooling);
+    ret += sscanf(line, "RadiativeCoolingModel = %"ISYM, &RadiativeCoolingModel);
     ret += sscanf(line, "GadgetEquilibriumCooling = %"ISYM, &GadgetEquilibriumCooling);
     ret += sscanf(line, "MultiSpecies = %"ISYM, &MultiSpecies);
     ret += sscanf(line, "PrimordialChemistrySolver = %"ISYM, &PrimordialChemistrySolver);
diff --git a/src/enzo/SetDefaultGlobalValues.C b/src/enzo/SetDefaultGlobalValues.C
index b09f1ce5b..968e924c6 100644
--- a/src/enzo/SetDefaultGlobalValues.C
+++ b/src/enzo/SetDefaultGlobalValues.C
@@ -308,6 +308,8 @@ int SetDefaultGlobalValues(TopGridData &MetaData)
   RandomForcingEdot           = -1.0;              //AK
   RandomForcingMachNumber     = 0.0;               //AK
   RadiativeCooling            = FALSE;             // off
+  RadiativeCoolingModel       = 1;                 //1=cool_rates.in table lookup
+                                                   //3=Koyama&Inutsuka 2002
   GadgetEquilibriumCooling    = FALSE;             // off
   RadiativeTransfer           = 0;                 // off
   RadiativeTransferFLD        = 0;                 // off
diff --git a/src/enzo/WriteAllData.C b/src/enzo/WriteAllData.C
index 3405787ad..df7d0b3e2 100644
--- a/src/enzo/WriteAllData.C
+++ b/src/enzo/WriteAllData.C
@@ -615,8 +615,8 @@ int WriteAllData(char *basename, int filenumber,
   CommunicationBarrier();
 
   if ( MyProcessorNumber == ROOT_PROCESSOR ){
-    sptr = fopen("OutputLog", "a");
-    fprintf(sptr, "DATASET WRITTEN %s \n", name);
+    sptr = fopen("OutputLogA", "a");
+    fprintf(sptr, "DATASET WRITTEN %s %8"ISYM" %18.16e\n", name, MetaData.CycleNumber, MetaData.Time);
     fclose(sptr);
   }
  
diff --git a/src/enzo/WriteParameterFile.C b/src/enzo/WriteParameterFile.C
index 20c886c48..46e03c6b5 100644
--- a/src/enzo/WriteParameterFile.C
+++ b/src/enzo/WriteParameterFile.C
@@ -410,6 +410,7 @@ int WriteParameterFile(FILE *fptr, TopGridData &MetaData)
   fprintf(fptr, "RandomForcing                  = %"ISYM"\n", RandomForcing);
   fprintf(fptr, "RandomForcingEdot              = %"GSYM"\n", RandomForcingEdot);
   fprintf(fptr, "RadiativeCooling               = %"ISYM"\n", RadiativeCooling);
+  fprintf(fptr, "RadiativeCoolingModel          = %"ISYM"\n", RadiativeCoolingModel);
   fprintf(fptr, "GadgetEquilibriumCooling       = %"ISYM"\n", GadgetEquilibriumCooling);
   fprintf(fptr, "MultiSpecies                   = %"ISYM"\n", MultiSpecies);
   fprintf(fptr, "PrimordialChemistrySolver      = %"ISYM"\n", PrimordialChemistrySolver);
diff --git a/src/enzo/cool1d_koyama.F b/src/enzo/cool1d_koyama.F
new file mode 100644
index 000000000..07c615127
--- /dev/null
+++ b/src/enzo/cool1d_koyama.F
@@ -0,0 +1,119 @@
+
+#include "fortran.def"
+c=======================================================================
+c////////////////////////  SUBROUTINE cool1d_koyama  \\\\\\\\\\\\\\\\\\\\\\\\\\
+c
+      subroutine cool1d_koyama(d,e,ge,u,v,w ,
+     &                  in,jn, kn, nratec, idual, idim, imethod, iter,
+     &                  is, ie, j,k, temstart, temend,
+     &                  utem, uxyz, urho, utim, gamma, coola,
+     &                  edot, tgas, tgasold, p2d, cool
+     &                  )
+
+    
+      implicit NONE
+c
+c  From Koyama & Inutsuka 2006, arxiv: asrto-ph 0605528v1
+c
+c  Arguments
+c
+      integer in, jn, kn, is, ie, j, k,
+     &        idual, nratec, idim, imethod, iter
+      real    temstart, temend, utem, uxyz, urho, utim,
+     &        gamma, coola(nratec)
+      real    d(in,jn,kn),   ge(in,jn,kn),     e(in,jn,kn),
+     &        u(in,jn,kn),    v(in,jn,kn),     w(in,jn,kn)
+c
+c  Parameters
+c
+      real pmin
+      parameter (pmin = tiny)
+      double precision mh
+
+      parameter (mh = 1.67262d-24)  ! DPC
+
+c
+c  Locals
+c
+      integer i
+      real thalf
+c
+c  Slice locals
+c 
+      real p2d(in), tgas(in), tgasold(in), cool(in)
+      real edot(in), kboltz, Lambda, Heating
+      parameter ( kboltz = 1.30865d-16 ) !erg / Kelvin
+
+c
+c     Note that the temperature definition here is slightly different
+c     than that in cool1d.src.
+c
+c     Compute Pressure
+c
+      if (imethod .eq. 2) then
+c
+c        Zeus - e() is really gas energy
+c
+         do i = is+1, ie+1
+            p2d(i) = (gamma - 1.0)*d(i,j,k)*e(i,j,k)
+         enddo
+      else
+         if (idual .eq. 1) then
+c
+c           PPM with dual energy -- use gas energy
+c
+            do i = is+1, ie+1
+               p2d(i) = (gamma - 1.0)*d(i,j,k)*ge(i,j,k)
+            enddo
+         else
+c
+c           PPM without dual energy -- use total energy
+c
+            do i = is+1, ie+1
+               p2d(i) = e(i,j,k) - 0.5*u(i,j,k)**2
+               if (idim .gt. 1) p2d(i) = p2d(i) - 0.5*v(i,j,k)**2
+               if (idim .gt. 2) p2d(i) = p2d(i) - 0.5*w(i,j,k)**2
+               p2d(i) = max((gamma - 1.0)*d(i,j,k)*p2d(i), tiny)
+            enddo
+         endif
+      endif
+c
+c     Compute temperature
+c
+      do i = is+1, ie+1
+         tgas(i) = max(p2d(i)*utem/d(i,j,k), temstart)
+      enddo
+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
+            tgasold(i) = tgas(i)
+         enddo
+      endif
+c
+c     Loop over a slice
+c
+      do i = is+1, ie+1
+        thalf = 0.5*(tgas(i) + tgasold(i))
+c        thalf = tgas(i)                              !kludge
+        edot(i) = 0
+        if( thalf .gt. temstart ) then
+            Heating = 2.d-26
+            Lambda = -(1d7 * exp(-1.184d5/(thalf+1000)) +
+     &          14d-3*sqrt(thalf)*exp(-92.0/thalf))
+            Lambda = Lambda * Heating
+            edot(i) = d(i,j,k)*urho/(mh*mh)*Lambda
+            edot(i) = edot(i) + Heating/mh
+            edot(i) = edot(i) * utim /(utem * kboltz)*mh
+
+        endif
+
+
+      enddo
+
+      do i=is+1, ie+1
+         tgasold(i) = tgas(i)
+      enddo
+      end
+c
diff --git a/src/enzo/global_data.h b/src/enzo/global_data.h
index 40f44b49d..f86083fa3 100644
--- a/src/enzo/global_data.h
+++ b/src/enzo/global_data.h
@@ -296,6 +296,7 @@ EXTERN float RootGridCourantSafetyNumber;
 
 EXTERN int RadiativeCooling;
 EXTERN CoolDataType CoolData;
+EXTERN int RadiativeCoolingModel;
 
 /* Cloudy cooling parameters and data. */
 
diff --git a/src/enzo/solve_cool.F b/src/enzo/solve_cool.F
index 5e883a38a..ee0985adb 100644
--- a/src/enzo/solve_cool.F
+++ b/src/enzo/solve_cool.F
@@ -6,7 +6,7 @@
 c
       subroutine solve_cool(
      &                d, e, ge, u, v, w,
-     &                in, jn, kn, nratec, iexpand, imethod,
+     &                in, jn, kn, nratec, iexpand, imethod, cool_method,
      &                idual, idim, igammah,
      &                is, js, ks, ie, je, ke, 
      &                dt, aye, temstart, temend, fh,
@@ -36,7 +36,7 @@ subroutine solve_cool(
 c  Arguments
 c
       integer in, jn, kn, is, js, ks, ie, je, ke,
-     &        idual, iexpand, nratec, idim, imethod, igammah
+     &       idual, iexpand, nratec, idim, imethod, igammah, cool_method
       real    dt, aye, temstart, temend, fh, gammaha,
      &        utem, uxyz, uaye, urho, utim,
      &        eta1, eta2, gamma, coola(nratec)
@@ -91,6 +91,7 @@ subroutine solve_cool(
 c
 c        Compute the cooling rate on a slice
 c
+         if( cool_method .eq. 1 ) then
             call cool1d(d, e, ge, u, v, w,
      &                  in, jn, kn, nratec, idual, idim, imethod, 
      &                  iter, igammah,
@@ -100,6 +101,16 @@ subroutine solve_cool(
      &                  indixe, t1, t2, logtem, tdef, edot,
      &                  tgas, tgasold, p2d, cool
      &                     )
+        else if ( cool_method .eq. 3) then
+            call cool1d_koyama(d,e,ge,u,v,w, 
+     &                  in,jn, kn, nratec, idual, idim, imethod, iter,
+     &                  is, ie, j,k,temstart, temend,
+     &                  utem, uxyz, urho, utim, gamma, coola,
+     &                  edot, tgas, tgasold, p2d, cool
+     &                  )
+
+        endif
+            
 c
 c         Compute maximim timstep that keeps any fractional change to 10%
 c              (maximum timestep is 1/2 hydro step).