diff --git a/dapper/mods/Stommel/clima.py b/dapper/mods/Stommel/clima.py
index 9027e307..94351ab7 100644
--- a/dapper/mods/Stommel/clima.py
+++ b/dapper/mods/Stommel/clima.py
@@ -25,22 +25,28 @@
                             BurnIn=0)  # 1 observation/year
 # Create default Stommel model
 model = stommel.StommelModel()
-#Heat air fluxes 
+#Switch on heat exchange with atmosphere. 
+#Start with default stationary atm. temperature. 
 functions = stommel.default_air_temp(N)
+#Add linear warming trend. Increase taking place over period T_warming and warming over this 
+#period is 3C over equator and 6C over pole. 
 trend = interp1d(np.array([0.,T_warming]), np.array([[0.,6.],[0.,3.]]), 
                  fill_value='extrapolate', axis=1)
 trended = [stommel.add_functions(func, trend) for func in functions]
+#Switch on heat flux. 
 model.fluxes.append(stommel.TempAirFlux(trended))
-#Salinity air fluxes 
+#Default stationary salinity flux with atm. 
 functions = stommel.default_air_salt(N)
+#Switch on atm. fluxes. 
 model.fluxes.append(stommel.SaltAirFlux(functions))
-#Melt flux 
+#Add fresh water from melting as constant rate over the period T_warming/
 melt_rate = -stommel.V_ice * np.array([1.0/(model.dx[0,0]*model.dy[0,0]), 0.0]) / T_warming #ms-1
 functions = stommel.default_air_ep(N)
+#Add melting rate to the default rate (=0) and switch Evaporation-Percipitation fluxe (EP). 
 functions = [stommel.merge_functions(T_warming, lambda t:func(t)+melt_rate, func)
              for func in functions]
 model.fluxes.append(stommel.EPFlux(functions))
-# Initial conditions
+# Default initial conditions
 x0 = model.x0
 #Variance Ocean temp[2], ocean salinity[2], temp diffusion parameter,
 #salt diffusion parameter, transport parameter
@@ -53,7 +59,7 @@
        }
 # Observational variances
 R = np.array([1.0, 1.0, 0.1, 0.1])**2  # C2,C2,ppt2,ppt2
-# Observation
+# Default bservation
 Obs = model.obs_ocean()
 # Create model.
 HMM = modelling.HiddenMarkovModel(Dyn, Obs, tseq, X0)
@@ -67,6 +73,7 @@
 fig, ax = stommel.time_figure(tseq)    
 stommel.plot_truth(ax, xx, yy)
 
+#Add equilibrium based on unperturbed initial conditions. 
 model.ens_member=0
 stommel.plot_eq(ax, tseq, model)
 
diff --git a/dapper/mods/Stommel/clima_forcing.py b/dapper/mods/Stommel/clima_forcing.py
index 65bcf5db..99454771 100644
--- a/dapper/mods/Stommel/clima_forcing.py
+++ b/dapper/mods/Stommel/clima_forcing.py
@@ -17,36 +17,47 @@ def exp_clima_forcing(N=100, seed=1000):
     T_warming=100*stommel.year
     #Data assimilation period in this or other experiments.
     Tda = 20 * stommel.year
-    # Timestepping. Timesteps of 1 day, running for 10 year.
+    # Timestepping. Timesteps of 1 day, running for 200 year.
     tseq = modelling.Chronology(stommel.year, 
                                 kko=np.array([],dtype=int), 
                                 T=200*stommel.year, BurnIn=0)  # 1 observation/year
     #Create model
     model = stommel.StommelModel()
     #Heat air fluxes 
+    #Start with default stationary atm. temperature. 
     functions = stommel.default_air_temp(N)
+    #Add linear warming with 6C/T_warming over the pole and 3C/T_warming over the equatior. 
     trend = interp1d(np.array([0.,T_warming]), np.array([[0.,6.],[0.,3.]]), 
                      fill_value='extrapolate', axis=1)
     trended = [stommel.add_functions(func, trend) for func in functions]
+    #Add random temperature perturbations with std dev. of 2C
     noised = [stommel.add_noise(func, seed=seed+n*20+1, sig=np.array([2.,2.])) 
               for n,func in enumerate(trended)]
+    #For time<Tda all ensemble member n uses noised[0] after that noised[n]
     functions = [stommel.merge_functions(Tda, noised[0], func) 
                  for func in noised]
+    #Activate surface heat flux. 
     model.fluxes.append(stommel.TempAirFlux(functions))
     #Salinity air fluxes 
     functions = stommel.default_air_salt(N)
+    #Add random salinity perturbations with std dev. of 0.2ppt
     noised = [stommel.add_noise(func, seed=seed+n*20+2, sig=np.array([.2,.2])) 
               for n,func in enumerate(functions)]
+    #For time<Tda all ensemble member n uses noised[0] after that noised[n]
     functions = [stommel.merge_functions(Tda, noised[0], func) 
                  for func in noised]
+    #Activate surface salinity flux. 
     model.fluxes.append(stommel.SaltAirFlux(functions))
     #Melt flux 
     melt_rate = -stommel.V_ice * np.array([1.0/(model.dx[0,0]*model.dy[0,0]), 0.0]) / T_warming #ms-1
+    #Default evaporation-percipitation flux (=0)
     functions = stommel.default_air_ep(N)
+    #Add effect Greenland melt with annual rate melt_rate
     functions = [stommel.merge_functions(T_warming, lambda t:func(t)+melt_rate, func)
                  for func in functions]
+    #Activate EP flux. 
     model.fluxes.append(stommel.EPFlux(functions))
-    # Initial conditions
+    #Default nitial conditions
     x0 = model.x0
     #Variance in initial conditions and parameters.
     B = stommel.State().zero()
@@ -61,7 +72,7 @@ def exp_clima_forcing(N=100, seed=1000):
            'model': model.step,
            'noise': 0
            }
-    # Observation
+    # Default observations. 
     Obs = model.obs_ocean()
     # Create model.
     HMM = modelling.HiddenMarkovModel(Dyn, Obs, tseq, X0)
@@ -83,6 +94,7 @@ def exp_clima_forcing(N=100, seed=1000):
     for n in range(np.size(Efor,1)):
         stommel.plot_truth(ax, Efor[:,n,:], yy)
         
+    #Add equilibrium based on unperturbed initial conditions. 
     model.ens_member=0
     stommel.plot_eq(ax, HMM.tseq, model, stommel.prob_change(Efor) * 100.)
     
diff --git a/dapper/mods/Stommel/clima_forcing_da.py b/dapper/mods/Stommel/clima_forcing_da.py
index ef303651..7d4d0d75 100644
--- a/dapper/mods/Stommel/clima_forcing_da.py
+++ b/dapper/mods/Stommel/clima_forcing_da.py
@@ -15,7 +15,7 @@ def exp_clima_forcing_da(N=10, seed=1000):
     Tda = 20*stommel.year
     #Time period over which climate change takes place
     T_warming=100*stommel.year
-    # Timestepping. Timesteps of 1 day, running for 10 year.
+    # Timestepping. Timesteps of 1 day, running for 200 year.
     kko = np.arange(1, int(Tda/stommel.year)+1)
     tseq = modelling.Chronology(stommel.year, kko=kko, 
                                 T=200*stommel.year, BurnIn=0)  # 1 observation/year
@@ -23,28 +23,38 @@ def exp_clima_forcing_da(N=10, seed=1000):
     model = stommel.StommelModel()
     #Heat air fluxes 
     functions = stommel.default_air_temp(N)
+    #Add linear warming with 6C/T_warming over the pole and 3C/T_warming over the equatior. 
     trend = interp1d(np.array([0.,T_warming]), np.array([[0.,6.],[0.,3.]]), 
                      fill_value='extrapolate', axis=1)
     trended = [stommel.add_functions(func, trend) for func in functions]
+    #Add random temperature perturbations with std dev. of 2C
     noised = [stommel.add_noise(func, seed=seed+n*20+1, sig=np.array([2.,2.])) 
               for n,func in enumerate(trended)]
+    #For time<Tda all ensemble member n uses noised[0] after that noised[n]
     functions = [stommel.merge_functions(Tda, noised[0], func) 
                  for func in noised]
+    #Activate surface heat flux. functions[n] contains atm. temperature for ensemble member n. 
     model.fluxes.append(stommel.TempAirFlux(functions))
     #Salinity air fluxes 
     functions = stommel.default_air_salt(N)
+    #Add random salinity perturbations with std dev. of 0.2ppt
     noised = [stommel.add_noise(func, seed=seed+n*20+2, sig=np.array([.2,.2])) 
               for n,func in enumerate(functions)]
+    #For time<Tda all ensemble member n uses noised[0] after that noised[n]
     functions = [stommel.merge_functions(Tda, noised[0], func) 
                  for func in noised]
+    #Activate surface salinity flux. 
     model.fluxes.append(stommel.SaltAirFlux(functions))
     #Melt flux 
     melt_rate = -stommel.V_ice * np.array([1.0/(model.dx[0,0]*model.dy[0,0]), 0.0]) / T_warming #ms-1
+    #Default evaporation-percipitation flux (=0)
     functions = stommel.default_air_ep(N)
+    #Add effect Greenland melt with annual rate melt_rate
     functions = [stommel.merge_functions(T_warming, lambda t:func(t)+melt_rate, func)
                  for func in functions]
+    #Activate EP flux. 
     model.fluxes.append(stommel.EPFlux(functions))
-    # Initial conditions
+    #Default initial conditions
     x0 = model.x0
     #Variance in initial conditions and parameters.
     B = stommel.State().zero()
@@ -59,7 +69,7 @@ def exp_clima_forcing_da(N=10, seed=1000):
            'model': model.step,
            'noise': 0
            }
-    # Observation
+    #Default observations. 
     Obs = model.obs_ocean()
     # Create model.
     HMM = modelling.HiddenMarkovModel(Dyn, Obs, tseq, X0)
@@ -80,6 +90,7 @@ def exp_clima_forcing_da(N=10, seed=1000):
     for n in range(np.size(Efor,1)):
         stommel.plot_truth(ax, Efor[:,n,:], yy)
         
+    #Add equilibrium based on unperturbed initial conditions. 
     model.ens_member=0
     stommel.plot_eq(ax, HMM.tseq, model, stommel.prob_change(Efor) * 100.)
     
diff --git a/dapper/mods/Stommel/clima_uncertainty.py b/dapper/mods/Stommel/clima_uncertainty.py
index fc642d3a..3060570f 100644
--- a/dapper/mods/Stommel/clima_uncertainty.py
+++ b/dapper/mods/Stommel/clima_uncertainty.py
@@ -19,51 +19,56 @@
 Tda = 20 * stommel.year
 #Time over which climate change will occur
 T_warming = 100 * stommel.year
-# Timestepping. Timesteps of 1 day, running for 10 year.
+# Timestepping. Timesteps of 1 day, running for 200 year.
 tseq = modelling.Chronology(stommel.year, 
                             kko=np.array([],dtype=int), 
                             T=200*stommel.year, BurnIn=0)  # 1 observation/year
 #Create model
 model = stommel.StommelModel()
-#Heat air fluxes 
+#Default stationary atm. temperature. 
 functions = stommel.default_air_temp(N)
 
 for n, func in enumerate(functions):
-    #Same heating rate for all members for t<=Tda, after that differences.
+    #Set same heating rate for all members for t<=Tda, after that differences.
     trend = interp1d(np.array([0.,Tda,T_warming]), 
                  np.array([[0.,1.2,np.random.normal(loc=6,scale=1.5)],
                            [0.,1.2,np.random.normal(loc=3,scale=0.5)]]),
                  fill_value='extrapolate', axis=1)
+    #Set function for ensemble member n. 
     functions[n]=stommel.add_functions(functions[n], trend)
-    
+
+#Add white noise with std dev of 2C over both pole and equator basin separately.   
 noised = [stommel.add_noise(func, seed=n*20+1, sig=np.array([2.,2.])) 
           for n,func in enumerate(functions)]
 functions = [stommel.merge_functions(Tda, noised[0], func) 
              for func in noised]
+#Switch on the atm. heat fluxes. 
 model.fluxes.append(stommel.TempAirFlux(functions))
 #Salinity air fluxes 
 functions = stommel.default_air_salt(N)
+#Add white with std dev. of .2 ppt. 
 noised = [stommel.add_noise(func, seed=n*20+2, sig=np.array([.2,.2])) 
           for n,func in enumerate(functions)]
 functions = [stommel.merge_functions(Tda, noised[0], func) 
              for func in noised]
+#Switch on salinity fluxes. 
 model.fluxes.append(stommel.SaltAirFlux(functions))
 #Melt flux 
 A = np.array([1,0])/(model.dx[0]*model.dy[0]) #area scaling for flux
     
 functions = []
 for func in range(N+1):
-    #Period over which ice melts.
+    #Period over which ice melts. Period varies with std. dev. of 10 years between ensemble members. 
     T1 = np.random.normal() * 10 * stommel.year + T_warming
     #Same melting rate for t<=Tda (rate0), after that differences (rate1)
     rate0 = -stommel.V_ice / T_warming * A
     rate1 = (-stommel.V_ice - rate0 * Tda) / (T1-Tda) * A
-    
+    #Set function for ensemble member n. 
     functions.append(interp1d(np.array([0,Tda,T1]), np.array([rate0, rate1, 0*A]), 
                               kind='previous', fill_value='extrapolate', axis=0))
-    
-#model.fluxes.append(stommel.EPFlux(functions))
-# Initial conditions
+#Switch on evaporation-percipitation flux. 
+model.fluxes.append(stommel.EPFlux(functions))
+# Default initial conditions
 x0 = model.x0
 #Variance Ocean temp[2], ocean salinity[2], temp diffusion parameter,
 #salt diffusion parameter, transport parameter
@@ -81,7 +86,7 @@
        }
 # Observational variances
 R = np.array([1.0, 1.0, 0.1, 0.1])**2  # C2,C2,ppt2,ppt2
-# Observation
+# Default observations
 Obs = model.obs_ocean()
 Obs['noise'] = modelling.GaussRV(C=R, mu=np.zeros_like(R))
 # Create model.
@@ -104,6 +109,7 @@
 for n in range(np.size(Efor,1)):
     stommel.plot_truth(ax, Efor[:,n,:], yy)
     
+#Add equilibrium based on unperturbed initial conditions. 
 model.ens_member=0
 stommel.plot_eq(ax, tseq, model, stommel.prob_change(Efor) * 100.)
 
diff --git a/dapper/mods/Stommel/ref_forcing.py b/dapper/mods/Stommel/ref_forcing.py
index 5a8859b2..35a5d14c 100644
--- a/dapper/mods/Stommel/ref_forcing.py
+++ b/dapper/mods/Stommel/ref_forcing.py
@@ -15,26 +15,32 @@
 
 def exp_ref_forcing(N=100, seed=1000):
     Tda = 20 * stommel.year #time period over which DA takes place. 
-    # Timestepping. Timesteps of 1 day, running for 10 year.
+    # Timestepping. Timesteps of 1 day, running for 200 year.
     tseq = modelling.Chronology(stommel.year, kko=np.array([], dtype=int), 
                                 T=200*stommel.year, BurnIn=0)  # 1 observation/year
     # Create default Stommel model
     model = stommel.StommelModel()
-    #Heat air fluxes 
+    #Switch on heat exchange with atmosphere. 
+    #Start with default stationary atm. temperature. 
     functions = stommel.default_air_temp(N)
+    #Add white noise with std dev of 2C over both pole and equator basin separately. 
     noised = [stommel.add_noise(func, seed=seed+n*20+1, sig=np.array([2.,2.])) 
               for n,func in enumerate(functions)]
     functions = [stommel.merge_functions(Tda, noised[0], func2) 
                  for func2 in noised]
+    #Switch on the atm. heat fluxes. 
     model.fluxes.append(stommel.TempAirFlux(functions))
-    #Salinity air fluxes 
+    #Switch on salinity exchange with atmosphere. 
+    #Start with default stationary atm. salinity. 
     functions = stommel.default_air_salt(N)
+    #Add white with std dev. of .2 ppt. 
     noised = [stommel.add_noise(func, seed=seed+n*20+2, sig=np.array([.2,.2])) 
               for n,func in enumerate(functions)]
     functions = [stommel.merge_functions(Tda, noised[0], func2) 
                  for func2 in noised]
+    #Switch on salinity fluxes. 
     model.fluxes.append(stommel.SaltAirFlux(functions))
-    # Initial conditions
+    #Default initial conditions
     x0 = model.x0
     #Variance Ocean temp[2], ocean salinity[2], temp diffusion parameter,
     #salt diffusion parameter, transport parameter
@@ -50,7 +56,7 @@ def exp_ref_forcing(N=100, seed=1000):
            'model': model.step,
            'noise': 0
            }
-    # Observation
+    #Default observation/
     Obs = model.obs_ocean()
     # Create model.
     HMM = modelling.HiddenMarkovModel(Dyn, Obs, tseq, X0)
@@ -71,6 +77,7 @@ def exp_ref_forcing(N=100, seed=1000):
     for n in range(np.size(Efor,1)):
         stommel.plot_truth(ax, Efor[:,n,:], yy)
         
+    #Add equilibrium based on unperturbed initial conditions. 
     model.ens_member=0
     stommel.plot_eq(ax, HMM.tseq, model, stommel.prob_change(Efor) * 100.)
     
diff --git a/dapper/mods/Stommel/ref_forcing_da.py b/dapper/mods/Stommel/ref_forcing_da.py
index 2fdafed0..a139660a 100644
--- a/dapper/mods/Stommel/ref_forcing_da.py
+++ b/dapper/mods/Stommel/ref_forcing_da.py
@@ -12,26 +12,31 @@
 import os
 
 def exp_ref_forcing_da(N=100, seed=1000):
-    # Timestepping. Timesteps of 1 day, running for 10 year.
+    # Timestepping. Timesteps of 1 day, running for 200 year.
     Tda = 20 * stommel.year #time period over which DA takes place. 
     kko = np.arange(1, int(Tda/stommel.year)+1)
     tseq = modelling.Chronology(stommel.year, kko=kko, 
                                 T=200*stommel.year, BurnIn=0)  # 1 observation/year
     # Create default Stommel model
     model = stommel.StommelModel()
-    #Heat air fluxes 
+    #Switch on heat exchange with atmosphere. 
+    #Start with default stationary atm. temperature.
     functions = stommel.default_air_temp(N)
+    #Add white noise with std dev of 2C over both pole and equator basin separately. 
     noised = [stommel.add_noise(func, seed=seed+n*20+1, sig=np.array([2.,2.])) 
               for n,func in enumerate(functions)]
     functions = [stommel.merge_functions(Tda, noised[0], func2) 
                  for func2 in noised]
+    #Switch on the atm. heat fluxes. 
     model.fluxes.append(stommel.TempAirFlux(functions))
     #Salinity air fluxes 
     functions = stommel.default_air_salt(N)
+    #Add white with std dev. of .2 ppt. 
     noised = [stommel.add_noise(func, seed=seed+n*20+2, sig=np.array([.2,.2])) 
               for n,func in enumerate(functions)]
     functions = [stommel.merge_functions(Tda, noised[0], func2) 
                  for func2 in noised]
+    #Switch on salinity fluxes. 
     model.fluxes.append(stommel.SaltAirFlux(functions))
     # Initial conditions
     x0 = model.x0
@@ -70,6 +75,7 @@ def exp_ref_forcing_da(N=100, seed=1000):
     for n in range(np.size(Efor,1)):
         stommel.plot_truth(ax, Efor[:,n,:], yy)
         
+    #Add equilibrium based on unperturbed initial conditions. 
     model.ens_member=0
     stommel.plot_eq(ax, HMM.tseq, model, stommel.prob_change(Efor) * 100.)
     
diff --git a/dapper/mods/Stommel/ref_init.py b/dapper/mods/Stommel/ref_init.py
index 86f08934..1c686689 100644
--- a/dapper/mods/Stommel/ref_init.py
+++ b/dapper/mods/Stommel/ref_init.py
@@ -14,14 +14,17 @@
 
 # Number of ensemble members
 N = 100
-# Timestepping. Timesteps of 1 day, running for 10 year.
+# Timestepping. Timesteps of 1 day, running for 200 year.
 kko = np.array([])
 tseq = modelling.Chronology(stommel.year, kko=kko, 
                             T=200*stommel.year, BurnIn=0)  # 1 observation/year
 # Create default Stommel model
 model = stommel.StommelModel()
+#Switch on heat exchange with atmosphere. Assume stationary air temperatures. 
 model.fluxes.append(stommel.TempAirFlux(stommel.default_air_temp(N)))
+#Switch on salinity exchange with atmosphere. Assume stationary air salinity. 
 model.fluxes.append(stommel.SaltAirFlux(stommel.default_air_salt(N)))
+#Use default initial conditions.
 default_init = model.init_state
 # Initial conditions
 x0 = model.init_state.to_vector()
@@ -36,7 +39,7 @@
        'model': model.step,
        'noise': 0
        }
-# Observation
+#Default observations. 
 Obs = model.obs_ocean()
 # Create model.
 HMM = modelling.HiddenMarkovModel(Dyn, Obs, tseq, X0)
@@ -47,10 +50,10 @@
     xx,yy=HMM.simulate()
     stommel.plot_truth(ax, xx, yy)
     
+#Plot equilibria values based on unperturbed initial conditions.
 model.ens_member=0
 stommel.plot_eq(ax, tseq, model)
     
-
 #Save figure 
 fig.savefig(os.path.join(stommel.fig_dir,'ocean_noise.png'),format='png',dpi=500)
 
diff --git a/dapper/mods/Stommel/single_run.py b/dapper/mods/Stommel/single_run.py
index 78f730d8..5645e445 100644
--- a/dapper/mods/Stommel/single_run.py
+++ b/dapper/mods/Stommel/single_run.py
@@ -12,17 +12,20 @@
 
 # Number of ensemble members
 N = 0
-# Timestepping. Timesteps of 1 day, running for 10 year.
+# Timestepping. Timesteps of 1 day, running for 400 year.
 tseq = modelling.Chronology(stommel.year*10, kko=np.array([],dtype=int), 
                             T=400*stommel.year, BurnIn=0)  # 1 observation/year
 # Create default Stommel model
 model = stommel.StommelModel()
+#Switch on heat exchange with atmosphere. Assume stationary air temperatures. 
 model.fluxes.append(stommel.TempAirFlux(stommel.default_air_temp(N)))
+#Switch on salinity exchange with atmosphere. Assume stationary air salinity. 
 model.fluxes.append(stommel.SaltAirFlux(stommel.default_air_salt(N)))
-# Initial conditions
+# Adjust default initial conditions.
 x0 = model.init_state
 x0.temp += np.array([[-1.,1.]])
 x0.salt += np.array([[0.,.4]]) #Move initial state away from equilibrium
+#Set initial conditions. 
 x0 = x0.to_vector()
 X0 = modelling.GaussRV(C=np.zeros_like(x0), mu=x0)
 #Print model information.