diff --git a/climlab/tests/test_domain2D.py b/climlab/tests/test_domain2D.py
index 1edb4df4..7c7bfda5 100644
--- a/climlab/tests/test_domain2D.py
+++ b/climlab/tests/test_domain2D.py
@@ -17,8 +17,6 @@ def test_state():
 @pytest.mark.fast
 def test_2D_EBM():
     '''Can we step forward a 2D lat/lon EBM?'''
-    #state = climlab.surface_state(num_lon=4)
-    #m = climlab.EBM_annual(state=state)
     m = climlab.EBM_annual(num_lon=4)
     m.step_forward()
     assert m.state.Ts.shape == (90, 4, 1)
@@ -28,8 +26,6 @@ def test_2D_EBM():
 @pytest.mark.fast
 def test_2D_EBM_seasonal():
     '''Can we step forward a 2D seasonal lat/lon EBM?'''
-    #state = climlab.surface_state(num_lon=4)
-    #m = climlab.EBM_annual(state=state)
     m = climlab.EBM_seasonal(num_lon=4)
     m.step_forward()
     assert m.state.Ts.shape == (90, 4, 1)
@@ -38,16 +34,62 @@ def test_2D_EBM_seasonal():
 
 @pytest.mark.fast
 def test_2D_insolation():
-    #state = climlab.surface_state(num_lon=4)
-    #m = climlab.EBM_annual(state=state)
     m = climlab.EBM_annual(num_lon=4)
     #  the answers are the mean of 1D insolation arrays
     #  the mean shouldn't change from 1D to 2D...
     #  there are exactly the same amount of each number in 2D array
     assert np.mean(m.subprocess['insolation'].insolation) == pytest.approx(299.30467670961832)
-    #assert np.mean(m.subprocess['insolation'].insolation) == 299.30467670961832
     sfc = m.domains['Ts']
     m.add_subprocess('insolation',
         climlab.radiation.P2Insolation(domains=sfc, **m.param))
     assert np.mean(m.subprocess['insolation'].insolation) == pytest.approx(300.34399999999999)
-    #assert np.mean(m.subprocess['insolation'].insolation) == 300.34399999999999
+
+@pytest.mark.fast
+def test_diurnal_cycle():
+    # We will create a model with a diurnal cycle and a longitude dimension
+    # and another model with no longitude dimension that uses daily average insolation.
+    # The zonal mean insolation should equal the daily average at each latitude 
+    # (same time of year)
+    num_lat = 30
+    num_lon = 1000
+    state = climlab.surface_state(num_lat=num_lat, 
+                              num_lon=num_lon, 
+                              water_depth=0.5)
+    deltat = 3600.  # one hour timestep
+    olr = climlab.radiation.AplusBT(name='OLR', 
+                                    state=state, 
+                                    timestep=deltat)
+    ins = climlab.radiation.InstantInsolation(name='Insolation', 
+                                            domains=state.Ts.domain, 
+                                            timestep=deltat)
+    asr = climlab.radiation.SimpleAbsorbedShortwave(name='ASR', 
+                                                    state=state, 
+                                                    insolation=ins.insolation, 
+                                                    timestep=deltat)
+    ebm = climlab.couple([olr,asr,ins], name='EBM')
+
+    state_daily = climlab.surface_state(num_lat=num_lat, 
+                              water_depth=0.5)
+    deltat_daily = 3600.*24.
+    olr_daily = climlab.radiation.AplusBT(name='OLR (daily)', 
+                                    state=state_daily, 
+                                    timestep=deltat_daily)
+    ins_daily = climlab.radiation.DailyInsolation(name='Insolation (daily)', 
+                                            domains=state_daily.Ts.domain, 
+                                            timestep=deltat_daily)
+    asr_daily = climlab.radiation.SimpleAbsorbedShortwave(name='ASR (daily)', 
+                                            state=state_daily, 
+                                            insolation=ins_daily.insolation, 
+                                            timestep=deltat_daily)
+    ebm_daily = climlab.couple([olr_daily, asr_daily, ins_daily], name='EBM (daily)')
+
+    for model in [ebm, ebm_daily]:
+        model.compute_diagnostics()
+    # zonal average diurnally-varying insolation should match daily mean insolation
+    assert np.allclose(ebm.ASR.mean(axis=1), ebm_daily.ASR)
+
+    # Now we will check to see if the average ASR over one day at each point in longitude
+    # is equal to ASR over a single one-day step in the daily insolation model
+    for model in [ebm, ebm_daily]:
+        model.integrate_days(1)
+    assert np.allclose(ebm.ti)  # how do we access the time mean again?
\ No newline at end of file