$\frac{\partial T_{ml}} {\partial t} \rho_w c_w h_{ml}=Q_{srf}$
Compute the time evolution of the ocean mixed layer temperature $T_{ml}(t)$ for different
Estimate the typical time scale for stationarity and select appropriate time step $\Delta t$ for model integration.
We rearrange the Taylor series $$f(t+\Delta t)=f(t)+\frac{\Delta t}{1!} f'(t)+O(\Delta t^2)$$ to obtain the time derivative of the function f $$f'(t)=\frac{f(t+\Delta t)-f(t)}{\Delta t}$$ whereas the truncation error was neglected.
We identify the mixed layer temperature $T_{ml}$ as the function $f$. Thus it follows for the differential equation
$$\frac{\partial T_{ml}} {\partial t} \rho_w c_w h_{ml}=Q_{srf}$$$$\frac{ T_{ml}(t+\Delta t)-T_{ml}(t)}{\Delta t}\rho_w c_w h_{ml}=Q_{srf}$$$$\frac{ T_{ml}(t+\Delta t)-T_{ml}(t)}{\Delta t}\rho_w c_w h_{ml}=Q^\downarrow + Q_{LW}^{\uparrow}$$The following equation allows to calculate the temperature $T_{ml}(t+\Delta t)$ from the previous time step $T_{ml}(t)$:
$$T_{ml}(t+\Delta t)=\frac{\Delta t}{\rho_w c_w h_{ml}} (Q^\downarrow + Q_{LW}^{\uparrow})+T_{ml}(t)$$Thus, it is straightforward to do the integration within one time loop.
import seawater
S=34 # salinityy [PSU}
T=0.0 # Temperature [deg C]
P=0 # Pressure in [dbar]
rho0=seawater.dens(S,T,P) # Density of sea water
cp=seawater.cp(S,T,P) # Heat capacity of sea water
print 'Density [kg/m^3]', rho0
print 'Heat capacity [J/kg/K]',cp
Density [kg/m^3] 1027.29893846 Heat capacity [J/kg/K] 3992.61688157
import scipy.constants
def outgoing(T,epsilon=0.97):
""""Calculate outgoing radiation from Stefan-Boltzmann law
Input parameters:
T temperature in K
epsilon infrared emissvity, default 0.97
return radiation in W/m^2
"""
return epsilon*scipy.constants.sigma*T**4
# Important variables (to play around with)
h_ml=10 # [m] mixed layer depth
T_ml_0=20 # [deg C] initial temperature
dt=60*60*24*10 # [s] time step, 10 days
# Not so important variables (treated as constants)
import seawater
rho0=seawater.dens(34.0,T_ml_0,0) # Density of sea water, salinity 34 PSU, pressur 0 dbar
cp=seawater.cp(34.0,T_ml_0,0) # Heat capacity of sea water
epsilon=0.97 # Infrared emissivity
# Constants
import scipy.constants
sig=scipy.constants.sigma #Stefan-Boltzmann constant 5.670373e-08 W m^-2 K^-4