diff --git a/doc/source/User_Guide/newtonian_cooling.txt b/doc/source/User_Guide/newtonian_cooling.txt
new file mode 100644
index 00000000..88dbc143
--- /dev/null
+++ b/doc/source/User_Guide/newtonian_cooling.txt
@@ -0,0 +1,40 @@
+.. _newtonian_cooling:
+
+Newtonian Cooling
+-----------------------
+
+We have added an initial implementation of Newtonian cooling to the code that adds a term to the thermal energy equation of the form
+
+.. math::
+   :label: newton
+
+   \frac{\partial \Theta}{\partial t}  = \frac{\Delta\Theta_{\mathrm{eq}}-\Theta}{\tau_\mathrm{eq}}.
+
+Here, :math:`\Delta\Theta_\mathrm{eq}` is a target temperature/entropy variation about the background state and :math:`\tau_\mathrm{eq}` is the cooling timescale.  Newtonian cooling can be turned on by setting ``newtonian_cooling=.true.`` in the ``physical_controls_namelist``, and the cooling time is similarly controlled by specifying the value of ``newtonian_cooling_time``.  
+
+At present, :math:`\Delta\Theta_\mathrm{eq}` is allowed to take one of two forms.  These are controlled by setting the ``newtonian_cooling_type`` variable in ``physical_controls_namelist``.  A value of 1 yields
+
+.. math::
+   :label: ncteq1
+   
+   \Delta\Theta_\mathrm{eq} = A,
+
+and a value of 2 yields
+
+.. math::
+   :label: ncteq2
+
+   \Delta\Theta_\mathrm{eq} = A\mathrm{cos}(\theta)\mathrm{sin}(\phi).
+
+The amplitude :math:`A` is controlled by setting the variable  ``newtonian_cooling_tvar_amp``.   As an example, to use Newtonian cooling, one might add the following four lines to the ``physical_controls_namelist``.
+
+::
+
+   &physical_controls_namelist
+    newtonian_cooling = .true.
+    newtonian_cooling_type = 1
+    newtonian_cooling_time = 1.0d0
+    newtonian_cooling_tvar_amp = 1.0d0
+   /
+
+
diff --git a/doc/source/User_Guide/under_development.rst b/doc/source/User_Guide/under_development.rst
index 3114f778..4617c1bc 100644
--- a/doc/source/User_Guide/under_development.rst
+++ b/doc/source/User_Guide/under_development.rst
@@ -34,8 +34,6 @@ scalars respectively.
 
 See `tests/chi_scalar` for example input files.
 
-
-
-
+.. include:: newtonian_cooling.txt
 
 
diff --git a/src/Physics/Controls.F90 b/src/Physics/Controls.F90
index b8f4cb9c..1de01030 100755
--- a/src/Physics/Controls.F90
+++ b/src/Physics/Controls.F90
@@ -80,14 +80,14 @@ Module Controls
     Logical :: newtonian_cooling = .false.  ! Turn newtonian_cooling on/off
     Integer :: newtonian_cooling_type = 1
     Real*8  :: newtonian_cooling_time = 1.0d22
-    Real*8  :: newtonian_delta_tvar_Eq = 0.0d0
+    Real*8  :: newtonian_cooling_tvar_amp = 0.0d0
 
     Namelist /Physical_Controls_Namelist/ magnetism, nonlinear, rotation, lorentz_forces, &
                 & viscous_heating, ohmic_heating, advect_reference_state, benchmark_mode, &
                 & benchmark_integration_interval, benchmark_report_interval, &
                 & momentum_advection, inertia, n_active_scalars, n_passive_scalars, &
                 & newtonian_cooling, newtonian_cooling_type, newtonian_cooling_time, &
-                & newtonian_delta_tvar_eq
+                & newtonian_cooling_tvar_amp
 
     !///////////////////////////////////////////////////////////////////////////
     !   Temporal Controls
diff --git a/src/Physics/Sphere_Physical_Space.F90 b/src/Physics/Sphere_Physical_Space.F90
index dc5e3ae6..b6b9b9e9 100755
--- a/src/Physics/Sphere_Physical_Space.F90
+++ b/src/Physics/Sphere_Physical_Space.F90
@@ -53,10 +53,11 @@ Subroutine Physical_Space_Init()
 
             If (newtonian_cooling_type .eq. 1) Then
                 ! No angular variation
+                If (my_rank .eq. 0) Write(6,*) 'Newtonian cooling is active.  Type = 1'
                 Do t = my_theta%min, my_theta%max
                     Do r = my_r%min, my_r%max
                         Do k =1, n_phi
-                            tvar_eq(k,r,t) = newtonian_delta_tvar_eq
+                            tvar_eq(k,r,t) = newtonian_cooling_tvar_amp
                         Enddo
                     Enddo
                 Enddo
@@ -65,10 +66,11 @@ Subroutine Physical_Space_Init()
 
             If (newtonian_cooling_type .eq. 2) Then
                 ! Angular variation (ell=1,m=1, motivated by hot Jupiters)
+                If (my_rank .eq. 0) Write(6,*) 'Newtonian cooling is active.  Type = 2'
                 Do t = my_theta%min, my_theta%max
                     Do r = my_r%min, my_r%max
                         Do k =1, n_phi
-                            tvar_eq(k,r,t) = newtonian_delta_tvar_eq*costheta(t)*sinphi(k)
+                            tvar_eq(k,r,t) = newtonian_cooling_tvar_amp*costheta(t)*sinphi(k)
                         Enddo
                     Enddo
                 Enddo