% set global definitions % \newcommand{\D}{\displaystyle} % \section{Surface Fluxes and Vertical Diffusion} \subsection{Surface Fluxes \label{surflux}} The bulk aerodynamic formulas are used to parameterize surface fluxes of zonal and meridional momentum (wind stress) $F_u$ and $F_v$, sensible heat $F_T$ and latent heat $L \, F_q$, where $F_q$ is the surface flux of moisture and $L$ is the latent heat of vaporisation $L_v$, or, depending on temperature, the latent heat of sublimation $L_s$: \begin{equation}\label{fluxes} \begin{array}{rcl} \D F_u & = & \D \rho \, C_m \, |\vec{v}| \, u \\ && \\ \D F_v & = & \D \rho \, C_m \, |\vec{v}| \, v \\ && \\ \D F_T & = & \D c_p \, \rho \, C_h \, |\vec{v}| \, (\gamma T - T_S ) \\ && \\ \D L \, F_q & = & \D L\, \rho \, C_h \, C_w \, |\vec{v}| \, (\delta q - q_S ) \end{array} \end{equation} All fluxes are positive in downward direction. $\rho$ denotes the density, $c_p$ is the specific heat for moist air at constant pressure ($c_p= c_{pd} \, [1+(c_{pv}/c_{pd}-1)\, q]$, where $c_{pd}$ and $c_{pv}$ are the specific heats at constant pressure for dry air and water vapor, respectively). $C_m$ is the drag coefficient, $C_h$ is the transfer coefficient for heat, $T_S$ is the surface temperature, $q_S$ is the surface specific humidity and $|\vec{v}|$ is the absolute value of the horizontal velocity at the lowermost level with a prescribed minimum (default= 1 m/s) to avoid numerical problems. The wetness factor $C_w$ accounts for different evaporation efficiencies due to surface characteristics (Section \ref{hydro}). $u$, $v$, $T$ and $q$ are the zonal and meridional wind components, the temperature and the specific humidity, respectively, of the lowermost model level. The factors $\gamma$ and ${\delta}$ are used to relate the model quantities to the respective near surface values. $\delta$ is set to 1 and $\gamma$ is set to give a potential temperature: \begin{equation}\label{gamma} \gamma = \left(\frac{p_S}{p}\right)^{\frac{R_d}{c_{pd}}} \end{equation} where p is the pressure of the lowermost model level, $p_S$ is the surface pressure and $R_d$ is the gas constant for dry air. While $\gamma$, $\rho$, $C_m$, $C_h$, $|\vec{v}|$, $T_S$ and $q_S$ apply to time level $t - \Delta t$, values for $u^{t+ \Delta t}$, $v^{t+ \Delta t}$, $T^{t+ \Delta t}$ and $q^{t+ \Delta t}$ are computed implicitly from the discretized tendency equations: \begin{equation} \begin{array}{rcccl} \D \frac{u^{t+\Delta t}-u^{t-\Delta t}}{2 \Delta t} & = & \D - \, \frac{1}{\rho \, \Delta z}\, F_u^{t+\Delta t} & = & \D - \, \frac {g \, \rho \, C_m \, |\vec{v}|}{p_S \, \Delta \sigma} \, u^{t + \Delta t} \\ &&&& \\ \D \frac{v^{t+\Delta t}-v^{t-\Delta t}}{2 \Delta t} & = & \D - \, \frac{1}{\rho \, \Delta z}\, F_v^{t+\Delta t} & = & \D - \, \frac {g \, \rho \, C_m \, |\vec{v}|}{p_S \, \Delta \sigma} \, v^{t + \Delta t} \\ &&&& \\ \D \frac{T^{t+\Delta t}-T^{t-\Delta t}}{2 \Delta t} & = & \D -\, \frac{1}{c_p \, \rho \, \Delta z} \, F_T^{t+\Delta t} & = & \D - \, \frac{g \, \rho \, C_h \, |\vec{v}|} {p_S \, \Delta \sigma} \, (\gamma T^{t + \Delta t} - T_S) \\ &&&& \\ \D \frac{q^{t+\Delta t}-q^{t-\Delta t}}{2 \Delta t} & = & \D - \, \frac{1}{\rho \, \Delta z} \, F_q^{t+\Delta t}& = & \D - \, \frac{g \, \rho \, C_h \, C_w \, |\vec{v}|} {p_S \, \Delta \sigma} \, (\delta q^{t + \Delta t} - q_S) \end{array} \end{equation} where $g$ is the gravitational acceleration and $\Delta \sigma = \Delta p/p_S $ is the thickness of the lowermost model layer. In addition to the tendencies, the surface fluxes of momentum, sensible and latent heat and the partial derivative of the sensible and the latent heat flux with respect to the surface temperature are computed: \begin{equation}\label{fluxes2} \begin{array}{rcl} \D F_u & = & \D \rho \, C_m \, |\vec{v}| \, u^{t+\Delta t} \\ && \\ \D F_v & = & \D \rho \, C_m \, |\vec{v}| \, v^{t+ \Delta t}\\ && \\ \D F_T & = & \D c_p \, \rho \, C_h \, |\vec{v}| \, (\gamma T^{t + \Delta t} - T_S ) \\ && \\ \D L \, F_q & = & \D L\, \rho \, C_h \, C_w \, |\vec{v}| \, (\delta q^{t+\Delta t} - q_S ) \\ && \\ \D \frac{\partial F_T}{\partial T_S} & = & \D - c_p \, \rho \, C_h \, |\vec{v}| \\ && \\ \D \frac{\partial (L \, F_q)}{\partial T_S} & = & \D - L\, \rho \, C_h \, C_w \, |\vec{v}| \, \frac{\partial q_S(T_S)}{\partial T_S} \end{array} \end{equation} The derivatives of the fluxes may be used, for examples, for an implicit calculation of the surface temperature (see Section \ref{surtemp}). \subsection*{ Drag and transfer coefficients} The calculation of the drag and the transfer coefficient $C_m$ and $C_h$ follows the method described in Roeckner et al.~(1992) for the ECHAM-3 model, which bases on the work of Louis (1979) and Louis et al.~(1982). A Richardson number dependence of $C_m$ and $C_h$ in accordance to the Monin-Obukhov similarity theory is given by \begin{equation} \begin{array}{rcl} C_m & = & \left( \frac{k}{\ln (z/z_0)}\right)^{2} \, f_m (Ri, z/z_0) \\ &&\\ C_h & = & \left( \frac{k}{\ln (z/z_0)}\right)^{2} \, f_h (Ri, z/z_0) \end{array} \end{equation} where $k$ is the von Karman constant ($k$ = 0.4) and $z_0$ is the roughness length, which depends on the surface characteristics (Section~\ref{landsurf} and Section~\ref{seasurf}). The Richardson number $Ri$ is defined as \begin{equation} Ri=\frac{g\, \Delta z \,(\gamma_E T - \gamma_E T_S)}{\gamma_E T \, |\vec{v}|^2} \end{equation} where $\gamma_E$ transfers temperatures to virtual potential temperatures to include the effect of moisture. \begin{equation}\label{gammaE} \gamma_E = \left(1- \left(\frac{R_v}{R_d}-1\right)\, q \right) \,\left(\frac{p_S}{p}\right)^{\frac{R_d}{c_{pd}}} \end{equation} where $q$ refers to the respective specific humidities and $R_v$ is the gas constant for water vapor. Different empirical formulas for stable ($Ri \ge 0$) and unstable ($Ri < 0$) situations are used. For the stable case, $f_m$ and $f_h$ are given by \begin{equation}\label{fmfh1} \begin{array}{rcl} \D f_m & = &\D \frac{1}{1+(2\,b\,Ri) /\sqrt{\,1+ d\, Ri}} \\ && \\ \D f_h & = &\D \frac{1}{1+(3\,b\,Ri) /\sqrt{\,1+ d\, Ri}} \end{array} \end{equation} while for the unstable case, $f_m$ and $f_h$ are \begin{equation}\label{fmfh2} \begin{array}{rcl} \D f_m & = & \D 1- \frac{2\,b\,Ri}{1+3\,b\,c\, [ \frac{k}{\ln(z/z_0+1)}]^2\sqrt{-Ri\, (z/z_0+1)}} \\ && \\ \D f_h & = &\D 1- \frac{3\,b\,Ri}{1+3\,b\,c\, [ \frac{k}{\ln(z/z_0+1)}]^2\sqrt{-Ri\, (z/z_0+1)}} \end{array} \end{equation} where $b$, $c$, and $d$ are prescribed constants and set to default values of $b$ = 5, $c$ = 5 and $d$ = 5. As in ECHAM-3 for unstable condition over oceans the empirical formula from Miller et al. (1992) is used to compute $C_h$ \begin{equation}\label{miller} C_h=C_{mn} \cdot (1-C_R^{\delta})^{1/\delta} \end{equation} with \begin{equation} C_R=\frac{0.0016\cdot (\Delta \Theta_v)^{1/3}}{C_{mn}\cdot |\vec{v}|} \end{equation} and \begin{equation} C_{mn}=\left(\frac{k}{\ln(z/z_0)}\right)^2 \end{equation} $\delta$ is set to 1.25. \subsection{Vertical Diffusion \label{vdiff}} Vertical diffusion representing the non resolved turbulent exchange is applied to the horizontal wind components $u$ and $v$, the potential temperature $\theta$ ($= T (p_S/p)^{R_d/c_{pd}}$) and the specific humidity $q$. The tendencies due to the turbulent transports are given by \begin{equation} \begin{array}{rcccl} \D \frac{\partial u}{\partial t} & = & \D \frac {1}{\rho}\frac{\partial J_u}{\partial z} & = & \D \frac {1}{\rho}\frac{\partial }{\partial z} ( \rho\, K_m \, \frac{\partial u}{\partial z}) \\ &&&& \\ \D \frac{\partial v}{\partial t} & = & \D \frac {1}{\rho}\frac{\partial J_v}{\partial z} & = & \D \frac {1}{\rho}\frac{\partial }{\partial z} ( \rho \, K_m \, \frac{\partial v}{\partial z} )\\ &&&& \\ \D \frac{\partial T}{\partial t} & = &\D \frac {1}{\rho}\frac{\partial J_T}{\partial z} & = & \D \frac {1}{\rho}\frac{\partial }{\partial z} ( \rho \, K_h \, (\frac{p}{p_S})^{R_d/c_{pd}}\,\frac{\partial \theta}{\partial z}) \\ &&&& \\ \D \frac{\partial q}{\partial t} & = & \D \frac {1}{\rho}\frac{\partial J_q}{\partial z} & = & \D \frac {1}{\rho}\frac{\partial }{\partial z}( \rho\, K_h \, \frac{\partial q}{\partial z} ) \end{array} \end{equation} where p is the pressure, $p_S$ is the surface pressure, $R_d$ is the gas constant for dry air and $c_{pd}$ is the specific heat for dry air at constant pressure. Here, the turbulent fluxes (positive downward) of zonal and meridional momentum $J_u$ and $J_v$, heat $c_{pd} \, J_T$ and moisture $J_q$ are parameterized by a linear diffusion along the vertical gradient with the exchange coefficients $K_m$ and $K_h$ for momentum and heat, respectively. $K_m$ and $K_h$ depend on the actual state (see below). As the effect of the surface fluxes are computed separately (Section \ref{surflux}), no flux boundary conditions for the vertical diffusion scheme are assumed at the top and the bottom of the atmosphere but the vertical diffusion is computed starting with initial values for $u$, $v$, $q$ and $T$ which include the tendencies due to the surface fluxes. As for the surface fluxes, the equations are formulated implicitely with exchange coefficients applying to the old time level. This leads to sets of linear equations for $u^{t+\Delta t}$, $v^{t+\Delta t}$, $T^{t+\Delta t}$ and $q^{t+\Delta t}$, which are solved by a back substitution method. \subsection*{Exchange coefficients} The calculation of the exchange coefficient $K_m$ and $K_h$ follows the mixing length approach as an extension of the similarity theory used to define the drag and transfere coefficients (Section \ref{surflux} and Roeckner et al.~1992): \begin{equation} \begin{array}{rcl} \D K_m & = & \D l_m^2\, \left| \frac{\partial\vec{v}}{\partial z} \right| \, f_m(Ri) \\ &&\\ \D K_h & = & \D l_h^2\, \left| \frac{\partial\vec{v}}{\partial z} \right| \, f_h(Ri) \end{array} \end{equation} where the functional dependencies of $f_m$ and $f_h$ on $Ri$ are the same as for $C_m$ and $C_h$ (Eq.~\ref{fmfh1} and Eq.~\ref{fmfh2}), except that the term \begin{equation} \left[\frac{k}{\ln(z/z_0+1)}\right]^2\sqrt{(z/z_0+1)} \end{equation} is replaced by \begin{equation} \frac{l^2}{(\Delta z)^{3/2} \, z^{1/2}}\left[ \left( \frac{z+\Delta z}{z}\right)^{1/3} -1 \right]^{3/2} \end{equation} The Richardson number $Ri$ is defined as \begin{equation} Ri=\frac{g}{\gamma T} \frac{\partial (\gamma_E T)}{\partial z} \left| \frac{\partial \vec{v}}{\partial z} \right|^{-2} \end{equation} with $\gamma$ from Eq.~\ref{gamma} and $\gamma_E$ from Eq.~\ref{gammaE}. According to Blackadar (1962), the mixing lengths $l_m$ and $l_h$ are given by \begin{equation} \begin{array}{rcl} \D \frac{1}{l_m} & = & \D \frac{1}{k\, z} +\frac{1}{\lambda_m} \\ &&\\ \D \frac{1}{l_h} & = & \D \frac{1}{k\, z} +\frac{1}{\lambda_h} \end{array} \end{equation} with $\lambda_h = \lambda_m\sqrt{(3 d)/2}$. The parameters $\lambda_m$ and $d$ are set to default values of $\lambda_m = 160~m$ and $d= 5$. \newpage \section{Horizontal Diffusion} The horizontal diffusion parameterization based on the ideas of Laursen and Eliasen (1989), which, in the ECHAM-3 model (Roeckner et al.~1992), improves the results compared with a ${\nabla}^k$ horizontal diffusion. The diffusion is done in spectral space. The contribution to the tendency of a spectral prognostic variable $X_n$ is \begin{equation} \frac{\partial X_n}{\partial t} = -k_X L_n X_n \end{equation} where $n$ defines the total wave number. $L_n$ is a scale selective function of the total wave number and is chosen such that large scales are not damped while the damping gets stronger with increasing $n$: \begin{equation} L_n = \left\{ \begin{array}{lcl} (n-n_{\star})^{\alpha} & \mbox{for} & n > n_{\star} \\ &&\\ 0 & \mbox{for} & n \le n_{\star} \end{array} \right. \end{equation} where $n_{\star}$ is a cut-off wave number. For T21 resolution the parameters $n_{\star}$ and $\alpha$ are set to default values of $n_{\star}$~=~15 and $\alpha$~=~2 similar to the ECHAM-3 model (Roeckner et al.~1992). The diffusion coefficient $k_X$ defines the timescale of the damping and depends on the variable. In the model, $k_X$ is computed from prescribed damping time scales $\tau_X$ for the smallest waves. Default values of $\tau_D$~=~0.2~days for divergence, $\tau_{\xi}$~=~1.1~days for vorticity and $\tau_T$~=~15.6~days for temperature and $\tau_q$~=~0.1~days for humidity are chosen, which are comparable with the respective values in the T21 ECHAM-3 model exept for humidity where here a considerable smaller value is used. In contrast to ECHAM-3 no level or velocity dependent additional damping is applied. For T42 resolution the respective defaults are: $n_{\star}$~=~16, $\alpha$~=~4, $\tau_D$~=~0.06~days, $\tau_{\xi}$~=~0.3~days, $\tau_T$~=~0.76~days and $\tau_q$~=~0.1~days. \newpage \section{Radiation} \subsection{Short Wave Radiation} The short wave radiation scheme bases on the ideas of Lacis and Hansen (1974) for the cloud free atmosphere. For the cloudy part, either constant albedos and transmissivities for high- middle- and low-level clouds may be prescribed or parameterizations following Stephens (1978) and Stephens et al.~(1984) may be used. The downward radiation flux density $F^{\downarrow SW}$ is assumed to be the product of the extrateristical solar flux density $E_0$ with different transmission factors for various processes: \begin{equation} F^{\downarrow SW}= \mu_0 \, E_0 \cdot {\cal T}_R \cdot {\cal T}_O \cdot {\cal T}_W \cdot {\cal T}_D \cdot {\cal T}_C \cdot {\cal R}_S \end{equation} Here, $\mu_0$ refers to the cosine of the solar zenith angle and the factor ${\cal R}_S$ incorporates different surface albedo values. The Indices of the transmissivities ${\cal T}$ denote Rayleigh scattering ($R$), ozone absorption ($O$), water vapor absorption ($W$) and absorption and scattering by aerosols (dust; $D$) and cloud droplets ($C$), respectively. $E_0$ and $\mu_0$ are computed following Berger (1978a, 1978b). The algorithm used is valid to 1,000,000 years past or hence. The numeric to compute $E_0$ and $\mu_0$ is adopted from the CCM3 climate model (Kiehl et al.~1996, coding by E.~Kluzek 1997). The calculation accounts for earths orbital parameters and the earths distance to the sun, both depending on the year and the time of the year. In default mode the model runs with daily averaged insolation but a diurnal cycle can be switched on. Following, for example, Stephens (1984) the solar spectral range is divided into two regions: (1) A visible and ultraviolet part for wavelengths $\lambda < 0.75$ $\mu$m with pure cloud scattering, ozone absorption and Rayleigh scattering, and without water vapor absorption. (2) A near infrared part for wavelengths $\lambda > 0.75$ $\mu$m with cloud scattering and absorption and with water vapor absorption. Absorption and scattering by aerosols is neglected in the present scheme. Dividing the total solar energy $E_0$ into the two spectral regions results in the fractions ${E_1}$~=~0.517 and $E_2$~=~0.483 for spectral ranges 1 and 2, respectively. \subsection*{Clear sky} For the clear sky part of the atmospheric column parameterizations following Lacis and Hansen (1974) are used for Rayleigh scattering, ozone absorption and water vapor absorption. {\bf Visible and ultraviolet spectral range ($\lambda < 0.75$ $\mu$m)} In the visible and ultraviolet range, Rayleigh scattering and ozone absorption are considered for the clear sky part. Rayleigh scattering is confined to the lowermost atmospheric layer. The transmissivity for this layer is given by \begin{equation} {\cal T}_{R1}=1 - \frac{0.219}{1+0.816\mu_0} \end{equation} for the direct beam, and \begin{equation} {\cal T}_{R1}=1 - 0.144 \end{equation} for the scattered part. Ozone absorption is considered for the Chappuis band in the visible ${\cal A}^{vis}$ and for the ultraviolet range ${\cal A}^{uv}$. The total transmissivity due to ozone is given by \begin{equation} {\cal T}_{O1} = 1 - {\cal A}^{vis}_O - {\cal A}^{uv}_O \end{equation} with \begin{equation} {\cal A}^{vis}_O = \frac{ 0.02118x}{1+0.042x+0.000323x^2} \end{equation} and \begin{equation} {\cal A}^{uv}_O=\frac{1.082x}{(1+138.6x)^{0.805}}+\frac{ 0.0658x}{1+(103.6x)^3} \end{equation} where the ozone amount traversed by the direct solar beam, $x$, is \begin{equation} x=M \; u_{O_3} \end{equation} with $u_{O_3}$ being the ozone amount [cm] in the vertical column above the considered layer, and $M$ is the magnification factor after Rodgers (1967) \begin{equation} M= \frac{35}{(1224 {\mu_0}^2 +1)^{\frac{1}{2}}} \end{equation} The ozone path traversed by diffuse radiation from below is \begin{equation} x^{*}=M \; u_{O_3}+\overline{M} \; (u_t -u_{O_3}) \end{equation} where $u_t$ is the total ozone amount above the main reflecting layer and $\overline{M}$=1.9 is the effective magnification factor for diffusive upward radiation. {\bf Near infrared ($\lambda > 0.75$ $\mu$m)} In the near infrared solar region absorption by water vapor is considered only. The transmissivity is given by \begin{equation} {\cal T}_{W2}=1-\frac{2.9 y}{(1+141.5y)^{0.635} + 5.925y} \end{equation} where $y$ is the effective water vapor amount [cm] including an approximate correction for the pressure and temperature dependence of the absorption and the magnification factor $M$. For the direct solar beam, $y$ is given by \begin{equation} y=\frac{M}{g} \int\limits^p_0 0.1 \; q \left(\frac{p}{p_0}\right)\left(\frac{T_0}{T}\right)^ {\frac{1}{2}} dp \end{equation} while for the reflected radiation reaching the layer from below, $y$ is \begin{equation} y=\frac{M}{g} \int\limits^{p_S}_0 0.1 \; q \left(\frac{p}{p_0}\right)\left(\frac{T_0}{T}\right)^ {\frac{1}{2}} dp + \frac{\beta_d}{g} \int\limits^{p_S}_{p} 0.1 \; q \left(\frac{p}{p_0}\right)\left(\frac{T_0}{T}\right)^ {\frac{1}{2}} dp \end{equation} with the acceleration of gravity $g$, the surface pressure $p_S$, a reference pressure $p_0$~=~1000~hPa, a reference temperature $T_0$~=~273~K, the specific humidity $q$ [kg/kg] and the magnification factor for diffuse radiation $\beta_d$~=~1.66. \subsection*{Clouds} Two possibilities for the parameterization of the effect of clouds on the short wave radiative fluxes are implemented: (1) prescribed cloud properties and (2) a parameterization following Stephens (1978) and Stephens et al. (1984), which is the default setup. {\bf Prescribed cloud properties} Radiative properties of clouds are prescribed depending on the cloud level. Albedos ${\cal R}_{C1}$ for cloud scattering in the visible spectral range ($\lambda < 0.75$ $\mu$m), and albedos ${\cal R}_{C2}$ for cloud scattering and absorptivities ${\cal A}_{C2}$ for cloud absorption in the near infrared part ($\lambda > 0.75$ $\mu$m) are defined for high, middle and low level clouds. The default values are listed in Table \ref{tabcl1}. {\protect \begin{table}[h] \begin{center} \begin{tabular}{|c|c|c|c|}\hline Cloud & Visible range &\multicolumn{2}{c|}{Near infrared} \\ Level & ${\cal R}_{C1}$ & ${\cal R}_{C2}$ & ${\cal A}_{C1}$ \\ \hline &&& \\ High & 0.15 & 0.15 & 0.05 \\ Middle & 0.30 & 0.30 & 0.10 \\ Low & 0.60 & 0.60 & 0.20 \\ \hline \end{tabular} \end{center} \caption{\label{tabcl1} Prescribed cloud albedos ${\cal R}_{C}$ and absorptivities ${\cal A}_{C}$} for spectral range 1 and 2 \end{table} } {\bf Default: Parameterization according to Stephens (1978) and Stephens et al. (1984)} Following Stephens (1978) and Stephens et al. (1984) cloud parameters are derived from the cloud liquid water path $W_L$ [g/m$^2$] and the cosine of the solar zenith angel $\mu_0$. In the visible and ultraviolet range cloud scattering is present only while in the near infrared both, cloud scattering and absorption, are parameterized. {\bf Visible and ultraviolet spectral range ($\lambda < 0.75$ $\mu$m)} For the cloud transmissivity ${\cal T}_{C1}$ Stephens parameterization for a non absorbing medium is applied: \begin{equation} {\cal T}_{C1}=1- \frac{\beta_1\tau_{N1}/\mu_0}{1+\beta_1\tau_{N1} /\mu_0} = \frac{1}{ 1+\beta_1 \tau_{N1}/\mu_0} \end{equation} $\beta_1$ is the backscatter coefficient, which is available in tabular form. In order to avoid interpolation of tabular values the following interpolation formula is used \begin{equation} \beta_1 = f_{b1} \; \sqrt{\mu_0} \end{equation} where the factor $f_{b1}$ comprises a tuning opportunity for the cloud albedo and is set to a default value of 0.0641 for T21L10 (0.02 T21L5 and 0.085 T42L10). $\tau_{N1}$ is an effective optical depth for which Stephens (1979) provided the interpolation formula \begin{equation} \tau_{N1}= 1.8336 \; (\log{W_L})^{3.963} \end{equation} which is approximated by \begin{equation} \tau_{N1}= 2\; (\log{W_L})^{3.9} \end{equation} to be used also for the near infrared range (see below). {\bf Near infrared ($\lambda > 0.75$ $\mu$m)} The transmissivity due to scattering and absorption of a cloud layer in the near infrared spectral range is \begin{equation} {\cal T}_{C2}=\frac{4u}{R} \end{equation} where u is given by \begin{equation} u^2=\frac{(1- \tilde{\omega}_0+2\; \beta_2 \; \tilde{\omega}_0)}{(1- \tilde{\omega}_0)} \end{equation} and R by \begin{equation} R=(u+1)^2 \exp{(\tau_{eff})} -(u-1)^2 \exp{(-\tau_{eff})} \end{equation} with \begin{equation} \tau_{eff}=\frac{\tau_{N2}}{\mu_0}\sqrt{(1- \tilde{\omega}_0)(1-\tilde{\omega}_0 + 2 \; \beta_2 \; \tilde{\omega}_0)} \end{equation} where the original formulation for the optical depth $\tau_{N2}$ by Stephens (1978) \begin{equation} \tau_{N2}=2.2346 \; (\log{W_L})^{3.8034} \end{equation} is, as for the visible range, approximated by \begin{equation} \tau_{N2}= 2 \; (\log{W_L})^{3.9} \end{equation} Approximations for the table values of the back scattering coefficient $\beta_2$ and the single scattering albedo $\tilde{\omega}_0$ are \begin{equation} \beta_2=\frac{f_{b2}\; \sqrt{\mu_0}}{\ln{(3+0.1\; \tau_{N2})}} \end{equation} and \begin{equation} \tilde{\omega}_0=1- f_{o2}\;\mu_0^2\;\ln{(1000/\tau_{N2})} \end{equation} where $f_{b2}$ and $f_{o2}$ provide a tuning of the cloud properties and are set to default values of $f_{b2}$=0.045 and $f_{o2}$=0.0045 for T21L10 (0.004 T21L5, 0.0048 T42L10). The scattered flux is computed from the cloud albedo ${\cal R}_{C2}$ which is given by \begin{equation} {\cal R}_{C2}=[\exp{(\tau_{eff})}-\exp{(-\tau_{eff})}] \; \frac{u^2- 1}{R} \end{equation} \subsection*{Vertical integration} For the vertical integration, the adding method is used (e.g. Lacis and Hansen 1974, Stephens 1984). The adding method calculates the reflection ${\cal R}_{ab}$ and transmission ${\cal T}_{ab}$ functions for a composite layer formed by combining two layers one (layer $a$) on top of the other (layer $b$). For the downward beam ${\cal R}_{ab}$ and ${\cal T}_{ab}$ are given by \begin{eqnarray}\label{LH31} {\cal R}_{ab} & = &{\cal R}_{a}+{\cal T}_{a}{\cal R}_b{\cal T}^{*}_a/(1- {\cal R}^*_a{\cal R}_b) \nonumber \\ {\cal T}_{ab} & = &{\cal T}_a{\cal T}_b/(1-{\cal R}^*_a{\cal R}_b) \end{eqnarray} where the denominator accounts for multiple reflections between the two layers. For illumination form below ${\cal R}^*_{ab}$ and ${\cal T}^*_{ab}$ are given by \begin{eqnarray}\label{LH32} {\cal R}^*_{ab} & = &{\cal R}^*_b+{\cal T}^*_b{\cal R}^*_a{\cal T}_b/(1- {\cal R}^*_a{\cal R}_b) \nonumber \\ {\cal T}^*_{ab} & = &{\cal T}^*_a{\cal T}_b/(1-{\cal R}^*_a{\cal R}_b) \end{eqnarray} The following four steps are carried out to obtain the radiative upward and downward fluxes at the boundary between two layers from which the total flux and the absorption (heating rates) are calculated: 1) ${\cal R}_l$ and ${\cal T}_l$, $l=1, L$ are computed for each layer and both spectral regions according to the parameterizations. 2) The layers are added, going down, to obtain ${\cal R}_{1,l}$ and ${\cal T}_{1,l}$ for $L=2,L+1$ and ${\cal R}^*_{1,l}$ and ${\cal T}^*_{1,l}$ for $L=2,L$. 3) Layers are added one at the time, going up, to obtain ${\cal R}_{L+1-l,L+1}$, $l=1,L-1$ starting with the ground layer, ${\cal R}_{L+1} = {\cal R}_S$ which is the surface albedo and ${\cal T}_{L+1}$=0. 4) The upward $F^{\uparrow SW}_l$ and downward $F^{\downarrow SW}_l$ short wave radiative fluxes at the interface of layer ($1,l$) and layer (l+1,L+1) are determined from \begin{eqnarray} F^{\uparrow SW}_l & = &{\cal T}_{1,l}\;{\cal R}_{l+1,L+1}/(1- {\cal R}^*_{1,l}\;{\cal R}_{l+1,L+1}) \nonumber \\ F^{\downarrow SW}_l & = &{\cal T}_{1,l}/(1- {\cal R}^*_{1,l}\;{\cal R}_{l+1,L+1}) \end{eqnarray} The net downward flux at level $l$, $F_l^{\updownarrow SW}$, is given by \begin{equation} F_l^{\updownarrow SW}=F_l^{\downarrow SW}-F_l^{\uparrow SW} \end{equation} Finally, the temperature tendency for the layer between $l$ and $l+1$ is computed: \begin{equation} \frac{\Delta T_{l+\frac{1}{2}}}{2\Delta t} = - \frac{g}{c_p \, p_S}\frac{F_{l+1}^{\updownarrow SW}-F_{l}^{\updownarrow SW}}{\Delta \sigma} \end{equation} \newpage \subsection{Long Wave Radiation} {\bf Clear sky} For the clear sky long wave radiation, the broad band emissivity method is employed (see, for example, Manabe and M\"oller 1961, Rodgers 1967, Sasamori 1968, Katayama 1972, Boer et al. 1984). Using the broad band transmissivities ${\cal T}_{(z,z^{\prime})}$ between level $z$ and level $z^{\prime}$, the upward and downward fluxes at level $z$, $F^{\uparrow LW}(z)$ and $F^{\downarrow LW}(z)$, are \begin{equation} \begin{array}{rcl} \D F^{\uparrow LW}(z) & = &\D {\cal A}_S \, B(T_S) {\cal T}_{(z,0)} + \int\limits_0^z B(T^{\prime}) \frac{\partial {\cal T}_{(z,z^{\prime})}}{\partial z^{\prime}} d z^{\prime} \\ & & \\ \D F^{\downarrow LW}(z) & = & \D \int\limits_{\infty}^z B(T^{\prime}) \frac{\partial {\cal T}_{(z,z^{\prime})}}{\partial z^{\prime}} d z^{\prime} \end{array} \end{equation} where $B(T)$ denotes the black body flux ($ B(T) = \sigma_{SB} T^4$) and ${\cal A}_S$ is the surface emissivity. The effect of water vapor, carbon dioxide and ozone is included in the calculations of the transmissivities ${\cal T}$ (with ${\cal T} = 1 - {\cal A}$, where ${\cal A}$ is the absoroptivity/emissivity). The transmissivities for water vapor ${\cal T}_{H_2O}$, carbon dioxide ${\cal T}_{CO_2}$ and ozone ${\cal T}_{O_3}$ are taken from Sasamori (1968): \begin{equation}\label{taus} \begin{array}{lcl} \D {\cal T}_{H_2O}& = & \D 1-0.846\;(u_{H_2O}+3.59 \cdot 10^{-5})^{0.243} -6.90\cdot 10^{-2} \\ && \\ \multicolumn{3}{l}{\mbox{for } u_{H_2O} < 0.01\mbox{ g, and }}\\ && \\ \D {\cal T}_{H_2O}& = & \D 1-0.240\log{(u_{H_2O}+0.010)}+0.622 \\ && \\ \multicolumn{3}{l}{\mbox{else.}}\\ &&\\ &&\\ \D {\cal T}_{CO_2}& = & \D 1-0.0825\;u_{CO_2}^{0.456}\\ && \\ \multicolumn{3}{l}{\mbox{for } u_{CO_2} \le 0.5 \mbox{ cm, and }}\\ && \\ \D {\cal T}_{CO_2}& = & \D 1-0.0461\log{(u_{CO_2})}+0.074 \\ && \\ \multicolumn{3}{l}{\mbox{else.}}\\ &&\\ &&\\ \D {\cal T}_{O_3} & = & \D 1-0.0122\log{(u_{O_3}+6.5 \cdot 10^{-4})}+0.0385 \end{array} \end{equation} where $u_{H_2O}$, $u_{CO_2}$ and $u_{O_3}$ are the effective amounts of water vapor, carbon dioxide and ozone, respectively, which are obtained from: \begin{equation} u(p, p^{\prime}) = \frac{f}{g} \int\limits_{p}^{p^{\prime}} q_X \left(\frac{p^{\prime\prime}}{p_0}\right) dp^{\prime\prime} \end{equation} where $q_X$ denotes the mixing ratios [kg/kg] of water vapor, carbon dioxide and ozone, respectively, $g$ is the gravitational acceleration, $p$ is pressure and $p_0$~=~1000 hPa is the reference pressure. The factor $f$ is used to transfer the units to g/cm$^2$ for $u_{H_2O}$ and cm-STP for $u_{CO_2}$ and cm-STP for $u_{O_3}$, which are used in Eq.~\ref{taus}. To account for the overlap between the water vapor and the carbon dioxide bands near 15~$\mu$m, the CO$_2$ absorption is corrected by a H$_2$O transmission at 15~$\mu$m, ${\cal T}_{H_2O}^{15\mu m}$, with ${\cal T}_{H_2O}^{15\mu m}$ given by \begin{equation} {\cal T}_{H_2O}^{15\mu m} = 1.33-0.832 \, (u_{H_2O} + 0.0286)^{0.26} \end{equation} Water vapour continuum absorption is parameterized by \begin{equation} {\cal T}_{H_2O}^{cont} = 1. - \exp(-k_{cont} u_{H_2O}) \end{equation} with a constant $k_{cont}$ (default =0.03 for T21L10, 0.035 T21L5,T42L10) {\bf Clouds} Clouds can be either treated as gray bodies with a prescribed cloud flux emissivity (grayness) or the cloud flux emissivity is obtained from the cloud liquid water contend. If the cloud flux emissivity (grayness) ${\cal A}^{cl}$ is externally prescribed, the value is attributed to each cloud layer. Otherwise, which is the default, ${\cal A}^{cl}$ is calculated from the cloud liquid water (e.g. Stephens 1984) \begin{equation} {\cal A}^{cl}=1.-\exp{(-\beta_d \; k^{cl} \; W_L)} \end{equation} where $\beta_d$~=~1.66 is the diffusivity factor, $k^{cl}$ is the mass absorption coefficent (with is set to a default value of 0.1~m$^2$/g (Slingo and Slingo 1991)) and $W_L$ is the cloud liquid water path. For a single layer between $z$ and $z^{\prime}$ with fractional cloud cover $cc$, the total transmissivity ${\cal T}^*_{(z, z^{\prime})}$ is given by \begin{equation} {\cal T}^*_{(z, z^{\prime})} = {\cal T}_{(z, z^{\prime})} \, (1 - cc \, {\cal A}^{cl}) \end{equation} where ${\cal T}_{(z, z^{\prime})}$ is the clear sky transmissivity. When there is more than one cloud layer with fractional cover, random overlapping of the clouds is assumed and ${\cal T}^*_{(z, z^{\prime})}$ becomes \begin{equation} {\cal T}^*_{(z, z^{\prime})} ={\cal T}_{(z, z^{\prime})} \, \prod_j (1 - cc_j \, {\cal A}^{cl}_{j}) \end{equation} where the subscript $j$ denotes the cloud layers. \subsection*{Vertical discretization} To compute the temperature tendency for a model layer resulting form the divergence of the radiative fluxes, the vertical discretization scheme of Chou et al. (2002) is used. The upward and downward fluxes, $F_l^{\uparrow LW}$ and $F_l^{\downarrow LW}$, at level $l$, which is the interface between two model layers, are computed from \begin{equation}\label{rad1} \begin{array}{rcllcl} \D F_l^{\uparrow LW}& = &\D \sum\limits_{l^\prime = l}^{L} B_{l^\prime+\frac{1}{2}} [{\cal T}^*_{(l , l^\prime)}- {\cal T}^*_{(l^\prime+1,l)}] & & \D \;\; l=1, \cdots , L\\ &&&&\\ && \D +{\cal T}^*_{(l ,L+1)} \; F_{L+1}^{\uparrow LW} & &\\ & & & &\\ &&&&\\ \D F_l^{\downarrow LW}& = &\D \sum\limits_{l^\prime=1}^{l- 1} B_{l^\prime+\frac{1}{2}} [{\cal T}^*_{(l^\prime+1,l)}-{\cal T}^*_{(l^\prime,l)}] & & \D \;\; l=2, \cdots , L+1 \end{array} \end{equation} where ${\cal T}^*_{(l,l^{\prime})}$ denotes the transmissivity of the layer from level $l$ to level $ l^{\prime}$ (see above) and $B_{l+\frac{1}{2}}$ is the black body flux for level $l+\frac{1}{2}$. The downward flux at the top of the atmosphere, $F_0^{\downarrow LW}$, and the upward flux at the surface ,$F_{L+1}^{\uparrow LW}$, are given by \begin{equation} \begin{array}{rcl} \D F_0^{\downarrow LW} & = & \D 0 \\ &&\\ \D F_{L+1}^{\uparrow LW} & = & \D {\cal A}_S \; B(T_S) + (1-{\cal A}_S) \; F_{L+1}^{\downarrow LW} \end{array} \end{equation} where ${\cal A}_S$ denotes the surface emissivity and $T_S$ is the surface temperature. Note, that for a more convenient discription of the scheme, $l+\frac{1}{2}$ denotes a so called full level, where the temperatures are defined. This may be in contrast to the convention in most of the other sections where a full level is indicated by $l$. Eqs.~\ref{rad1} can be rearranged to give \begin{equation}\label{vertical1} \begin{array}{rclcl} \D F_l^{\uparrow LW}& = &\D B_{l+\frac{1}{2}} + \sum\limits_{l^\prime=l+1}^{L+1} {\cal T}^*_{(l^\prime , l)} \, [B_{l^\prime+\frac{1}{2}} - B_{l^\prime- \frac{1}{2}}] & & \D \;\; l=1, \cdots , L\\ &&&&\\ & & + {\cal T}^*_{(l,L+1)} \, (1-{\cal A}_S) \; F_{L+1}^{\downarrow LW} & & \\ & & & & \\ &&&&\\ \D F_{l^\prime}^{\downarrow LW}& = & \D B_{l^\prime- \frac{1}{2}} - \sum\limits_{l=1}^{l^\prime-1}{\cal T}^*_{(l^\prime , l)} \, [B_{l+\frac{1}{2}} - B_{l-\frac{1}{2}}] & & \D \;\; l^\prime=2, \cdots , L+1 \end{array} \end{equation} with the boundary conditions \begin{equation} \begin{array}{rcl} \D B_{L+\frac{3}{2}}& = & \D {\cal A}_S \, B(T_S) \\ &&\\ \D B_{\frac{1}{2}} & = & \D 0 \end{array} \end{equation} The net downward flux at level $l$, $F_l^{\updownarrow LW}$, is given by \begin{equation} F_l^{\updownarrow LW}=F_l^{\downarrow LW}-F_l^{\uparrow LW} \end{equation} Finally, the temperature tendency for the layer between $l$ and $l+1$ is computed: \begin{equation} \frac{\Delta T_{l+\frac{1}{2}}}{2\Delta t} = - \frac{g}{c_p \, p_S}\frac{F_{l+1}^{\updownarrow LW}-F_{l}^{\updownarrow LW}}{\Delta \sigma} \end{equation} {\bf Emission of a layer} As pointed out by Chou et al.~(2002), the difference between the upward and downward emission of a layer will be large, if the layer is rather opaque and the temperature range across the layer is large. This, in particular, holds for coarse vertical resolution as in the default version of the model. Therefore, the upward and the downward emission of a layer is computed separately following the ideas of Chou et al.~(2002): The contribution of the upward flux at level $p$ from the adjecant layer below can be written as \begin{equation} \label{FUPLW} \Delta F^{\uparrow LW}(p) = -\int\limits^{p+\Delta p}_{p} B(p^{\prime}) \; \frac{\partial {\cal T}_{(p,p^{\prime})}}{\partial p^{\prime}} \; dp^{\prime} = B^u \; (1-{\cal T}_{ (p+\Delta p, p)}) \end{equation} where $\Delta p$ is the thickness of the adjacent layer, $B^u$ is the effective Planck flux for the adjacent layer, and ${\cal T}_{(p+\Delta p, p)}$ is the flux transmittance between $p$ and $p +\Delta p$. Assuming that the Planck function varies linearly with pressure and the transmittance decreases exponentially with pressure away from $p$ it follows \begin{equation} B(p^{\prime})= B(p) + \frac{(B(p)-B(p+\Delta p)) (p^{\prime} - p)}{\Delta p} \end{equation} and \begin{equation} {\cal T}_{(p, p^{\prime})} = \exp{(-c\; (p^{\prime}-p))} \end{equation} with $c$ ia a constant. From Eq.~\ref{FUPLW} the effective Planck flux for the adjacent layer $B^u$ is \begin{equation} B^u=\frac{B(p)-B(p+\Delta p)\;{\cal T}_{(p+\Delta p, p)}}{1-{\cal T}_{(p+\Delta p, p)}} +\frac{B(p)- B(p+\Delta p)}{\ln({\cal T}_{(p+\Delta p, p)})} \end{equation} Similarly, for the downward flux at the lower boundary of the layer, the effective Planck function of the layer $B^d$ is \begin{equation} B^d=\frac{B(p+\Delta p)-B(p)\;{\cal T}_{(p+\Delta p, p)}}{1-{\cal T}_{(p+\Delta p, p)}} +\frac{B(p+\Delta p)-B(p)}{\ln({\cal T}_{(p+\Delta p, p)})} \end{equation} Replacing the respective Planck functions in Eqs.~\ref{vertical1} by $B^u$ and $B^d$ results in \begin{equation}\label{vertical2} \begin{array}{rclcl} \D F_l^{\uparrow LW}& = &\D B^u_{l+\frac{1}{2}} + \sum\limits_{l^\prime=l+1}^{L+1} {\cal T}^*_{(l^\prime , l)} \, [B^u_{l^\prime+\frac{1}{2}} - B^u_{l^\prime- \frac{1}{2}}] & & \D \;\; l=1, \cdots , L\\ &&&&\\ & & \D + {\cal T}^*_{(l,L+1)} \, (1-{\cal A}_S) \; F_{L+1}^{\downarrow LW} & & \\ & & & & \\ &&&&\\ \D F_{l^\prime}^{\downarrow LW}& = & \D B^d_{l^\prime- \frac{1}{2}} - \sum\limits_{l=1}^{l^\prime-1} {\cal T}^*_{(l^\prime , l)} \, [B^d_{l+\frac{1}{2}} - B^d_{l-\frac{1}{2}}] & & \D \;\; l^\prime=2, \cdots , L+1 \end{array} \end{equation} where \begin{equation} \begin{array}{lcl} \D B^d_{l^{\prime}-\frac{1}{2}} & = & \D \frac{B_{l^{\prime}}-B_{l^{\prime}-1} \; {\cal T}_{(l^{\prime},l^{\prime}-1)}}{1-{\cal T}_{(l^{\prime},l^{\prime}-1)}} + \frac{B_{l^{\prime}} - B_{l^{\prime}-1}}{\ln({\cal T}_{(l^{\prime},l^{\prime}-1)})} \\ && \\ \D B^u_{l^{\prime}-\frac{1}{2}} & = & \D (B_{l^{\prime}} + B_{l^{\prime}-1} ) - B^d_{l^{\prime}-\frac{1}{2}} \end{array} \end{equation} For the calculation of the effective Plank function, the mean transmissivity for a layer partially filled with clouds is given by \begin{equation} {\cal T}_{(l^{\prime},l^{\prime}-1)} = f_{{\cal T}} \; {\cal T}^{cs}_{(l^{\prime},l^{\prime}-1)} \; (1 - cc_{(l^{\prime},l^{\prime}-1)}{\cal A}^{cl}_{(l^{\prime},l^{\prime}-1)}) \end{equation} with the cloud emissivity ${\cal A}^{cl}$ and the clear sky transmissivity ${\cal T}^{cs}$ being defined above, and the factor $f_{{\cal T}}$ provides a tuning opportunity. When a model layer spans a region where the temperature lapse rate changes signs, the linearity of $B$ with respect to $p$ can not longer be assumed and $B^d$ and $B^u$ are simply computed from \begin{equation} B^u_{l+\frac{1}{2}}=B^d_{l-\frac{1}{2}}= 0.5 \; B_{l+\frac{1}{2}} + 0.25 \; (B_{l} + B_{l^{\prime}}) \end{equation} \subsection{Ozone} Ozone concentration is prescribed. Either a three dimensional ozone distribution can be externally provided or an idealized annual cycle of ozone concentration can be used. The idealized distribution bases on the analytic ozone distribution of Green (1964): \begin{equation} u_{O_3}(h)=\frac{a+a \; \exp{(-b/c)}}{1+\exp((h-b)/c)} \end{equation} where $u_{O_3}(h)$ is the ozone amount [cm-STP] in a vertical column above the altitude $h$, $a$ is the total ozone amount in a vertical column above the ground, $b$ the altitude at which the ozone concentration has its maximum. While for $a$~=~0.4~cm, $b$~=~20~km and $c$~=~5~km this distribution fits close to the mid-latitude winter ozone distribution, an annual cycle and a latitudinal dependence is introduced by varying $a$ with time $t$ and latitude $\phi$: \begin{equation} a(t,\phi)=a0+a1\cdot|\sin(\phi)|+ac\cdot \sin(\phi) \cdot \cos(2\pi(d-doff)/ndy) \end{equation} where $d$ is the actual day of the year, $doff$ an offset and $ndy$ the number of days per year. The defaults for the involved parameters are: $a0=0.25$, $a1=0.11$ and $ac=0.08$. \subsection{Additional Newtonian cooling} For the standard setup with a vertical resolution of five equally spaced sigma-levels, the model produces a strong bias in the stratospheric (uppermost level) temperatures. This may be attributed to the insufficient representation of the stratosphere and its radiative and dynamical processes. The bias also effects the tropospheric circulation leading, for example, to a misplacement of the dominant pressure centers. To enable the simulation of a more realistic tropospheric climate, a Newtonian cooling can be applied to the uppermost level. Using this method, the model temperature $T$ is relaxed towards a externally given distribution of the temperature $T_{NC}$ which results in additional temperature tendencies $\dot{T}$ for the uppermost model level of \begin{equation} \dot{T}=\frac{T_{NC}-T}{\tau_{NC}} \end{equation} where $\tau_{NC}$ is the time scale of the relaxation, which has a default value of ten days. \newpage \section{Moist Processes and Dry Convection} \subsection{Correction of Negative Humidity} Local negative values of specific humidity are an artifact of spectral models. In the model, a simple procedure corrects these negative values by conserving the global amount of water. The correction of negative moisture is performed at the beginning of the grid-point parameterization scheme. A negative value of specific humidity is reset to zero. Accumulation of all corrections defines a correction factor. A hierarchical scheme of three steps is used. First, the correction is done within an atmospheric column only. If there are atmospheric columns without sufficient moisture, a second correction step is done using all grid points of the respective latitude. Finally, if there is still negative humidity remaining, a global correction is performed. \subsection{Saturation Specific Humidity} For parameterizations of moist processes like cumulus convection and large scale condensation the computation of the saturation specific humidity $q_{sat}(T)$ and its derivative with respect to temperature $dq_{sat}(T)/dT$ is needed at several places. In the model, the Tetens formula (Lowe 1977) is used to calculate the saturation pressure $e_{sat} (T)$ and its derivative with respect to temperature $de_{sat}(T)/dT$: \begin{equation} \begin{array}{rcl} \D e_{sat}(T) & = & \D a_1 \exp{\left(a_2 \, \frac{T-T_0}{T-a_3}\right)} \\ && \\ \D \frac{de_{sat}(T)}{dT} & = & \D \frac{a_2 \, (T_0 - a_3)}{(T-a_3)^2} \, e_{sat}(T) \end{array} \end{equation} with the constants $a_1$~=~610.78, $a_2$~=~17.2693882, $a_3$~=~35.86 and $T_0$~=~273.16. The saturation specific humidity $q_{sat}(T)$ and its derivative $dq_{sat}(T)/dT$ are given by \begin{equation}\label{qsat} \begin{array}{rcl} \D q_{sat}(T) & = & \D \frac{\epsilon \, e_{sat}(T)}{p-(1-\epsilon )\, e_{sat} (T)} \\ &&\\ \D \frac{dq_{sat}(T)}{dT} & = & \D \frac{p \, q_{sat}(T)}{p-(1-\epsilon)\, e_{sat} (T)} \frac{de_{sat}(T)}{dT}\\ \end{array} \end{equation} where $p$ is the pressure and $\epsilon$ is the ration of the gas constants for dry air $R_d$ and water vapor $R_v$ ($\epsilon = R_d / R_v$). \subsection{Cumulus Convection} The cumulus convection is parameterized by a Kuo-type convection scheme (Kuo 1965, 1974) with some modifications to the original Kuo-scheme. The Kuo-scheme considers the effect of cumulus convection on the large scale flow applying the following assumptions. Cumulus clouds are forced by mean low level convergence in regions of conditionally unstable stratification. The production of cloud air is proportional to the net amount of moisture convergence into one grid box column plus the moisture supply by surface evaporation. In a modification to the original scheme, the implemented scheme also considers clouds which originate at upper levels where moisture convergence is observed. This type of cloud may occur in mid-latitude frontal regions. Therefore, only the moisture contribution which takes place in the layer between the lifting level and the top of the cloud is used instead of the whole column. Thus, the total moisture supply $I$ in a period $2 \Delta t$ is given by \begin{equation}\label{cli} I= \frac{2 \Delta t \, p_S}{g} \int\limits_{\sigma_{Top}}^{\sigma_{Lift}} A_q \, d \sigma \end{equation} where $A_q$ is the moisture convergence plus the surface evaporation if the lifting level $\sigma_{Lift}$ is the lowermost model level. $\sigma_{Top}$ is the cloud top level, $p_S$ is the surface pressure and $g$ is the gravitational acceleration. Lifting level, cloud base and cloud top are determined as follows. Starting form the lowermost level, the first level with positive moisture supply $A_q$ is considered as a lifting level. If the lowermost level $L$ is considered to be a lifting level and the surface layer is dry adiabatic unstable ($\theta_S > \theta_L$ where $\theta$ denotes the potential temperature), the convection starts from the surface. Air from the lifting level ($l+1$) is lifted dry adiabatically up to the next level ($l$) by keeping its specific humidity. A cloud base is assumed to coincide with level $l+\frac{1}{2}$ if the air is saturated at $l$. Above the cloud base the air is lifted moist adiabatically. Distribution of temperature $T_{cl}$ and of moisture $q_{cl}$ in the cloud is found by first lifting the air dry adiabatically \begin{equation}\label{clad} \begin{array}{rcl} \D (T_{cl})_l^{Ad} & = & \D (T_{cl})_{l+1} \left(\frac{\sigma_l}{\sigma_{l+1}}\right)^{\frac{R_ d}{c_{pd}}} \\ &&\\ \D (q_{cl})_l^{Ad} & = & \D (q_{cl})_{l+1} \end{array} \end{equation} and then by correcting temperature and moisture values due to the condensation of water vapor \begin{equation}\label{ccc1} \begin{array}{rcl} \D (T_{cl})_l & = & \D (T_{cl})_l^{Ad} + \frac{L}{c_p} \, \frac{(q_{cl})_l^{Ad} - q_{sat} [(T_{cl})_l^{Ad}]}{1+\frac{L}{c_p}\, \frac{dq_{sat}[(T_{cl})_l^{Ad}]}{dT}} \\ &&\\ \D (q_{cl})_l & = & \D (q_{cl})_l^{Ad} -\frac{(q_{cl})_l^{Ad} - q_{sat}[(T_ {cl})_l^{Ad}]}{1+\frac{L}{c_p}\, \frac{dq_{sat}[(T_{cl})_l^{Ad}]}{dT}} \end{array} \end{equation} where the suturation specific humidity $q_{sat}$ and its derivative with respect to temperature $dq_{sat}/dT$ are computed from Eqs.~\ref{qsat}. $L$ is either the latent heat of vapourisation $L_v$ or the latent heat of sublimation $L_s$ depending on the temperature. $c_p$ is the specific heat for moist air at constant pressure ($c_p= c_{pd} \, [1+(c_{pv}/c_{pd}-1)\, q]$ where $c_{pd}$ and $c_{pv}$ are the specific heats at constant pressure for dry air and water vapor, respectively) and $R_d$ in Eq.~\ref{clad} is the gas constant for dry air. For reasons of accuracy the calculation (\ref{ccc1}) is repeated once where $(T_{cl})^{Ad}$ and $(q_{cl})^{Ad}$ are now replaced by the results of the first iteration. Cumulus clouds are assumed to exist only if the environmental air with temperature $T_e$ and moisture $q_e$ is unstable stratified with regard to the rising cloud parcel: \begin{equation} (T_{cl})_l > (T_e)_l \end{equation} The top of the cloud $\sigma_{Top}$ is then defined as \begin{equation} \sigma_{Top}=\sigma_{l+\frac{1}{2}} \; \mbox{if } \left\{ \begin{array}{lcll} (T_{cl})_l & \le & (T_{e})_l & \mbox{and} \\ &&& \\ (T_{cl})_{l+1} & > & (T_{e})_{l+1} & \end{array} \right. \end{equation} Cumulus clouds do exist only if the net moisture accession $I$ as given by Eq.~\ref{cli} is positive. Once this final check has been done, the heating and moistening of the environmental air and the convective rain are computed. In the model either the original scheme proposed by Kuo (1968) or the modified scheme with the parameter $\beta$ (Kuo 1974) can be chosen, where $\beta$ determines the partitioning of heating and moistening of the environmental air. In the scheme without $\beta$ the surplus $P$ of total energy of the cloud against the environmental air is given by \begin{equation} P=\frac{p_s}{g} \int\limits_{\sigma_{Top}}^{\sigma_{Base}} (c_p\, (T_{cl} -T_{e}) + L\, (q_{sat}(T_e)-q_{e})) d\sigma \end{equation} The clouds produced dissolve instantaneously by artificial mixing with the environmental air, whereby the environment is heated and moistened by \begin{equation}\label{handm} \begin{array}{rcl} \D (\Delta T)^{cl} & = & \D a \, (T_{cl} -T _e) \\ &&\\ \D (\Delta q)^{cl} & = & \D a \, (q_{sat}(T_e) -q _e) \end{array} \end{equation} where $a$ is the fractional cloud area being produced by the moisture supply: \begin{equation} a=L\, \frac{I}{P} \end{equation} In the scheme with $\beta$ the fraction 1-$\beta$ of the moisture is condensed, while the remaining fraction $\beta$ is stored in the atmosphere. The parameter $\beta$ depends on the mean relative humidity and, in the present scheme, is given by \begin{equation} \beta = \left( 1 - \frac{1}{\sigma_{Base}-\sigma_{Top}} \int\limits_{\sigma_{Top}}^{\sigma_{Base}} \frac{q_e}{q_{sat}(T_e)} d\sigma \right)^3 \end{equation} Instead of Eq.~\ref{handm}, the temperature and moisture tendencies are now \begin{equation}\label{handm2} \begin{array}{rcl} \D (\Delta T)^{cl} & = & \D a_T \, (T_{cl} -T _e) \\ &&\\ \D (\Delta q)^{cl} & = & \D a_q \, (q_{sat}(T_e) -q _e) \end{array} \end{equation} where $a_T$ and $a_q$ are given by \begin{equation} \begin{array}{rcl} \D a_T & = & \D \frac{(1-\beta )\, L\, I }{c_p\, \frac{p_S}{g}\, \int\limits_{\sigma_{Top}}^{\sigma_{Base}} (T_{cl} - T_e)\, d\sigma} \\ &&\\ \D a_q & = & \D \frac{\beta \, I }{\frac{p_S}{g}\, \int\limits_{\sigma_{Top}}^{\sigma_{Base}} (q_{sat}(T_e) - q_e) \, d\sigma} \end{array} \end{equation} The final tendencies for moisture $\partial q / \partial t$ and temperature $\partial T / \partial t$ which enter the diabatic leap frog time step are given by \begin{equation} \begin{array}{rcl} \D \frac{\partial q}{\partial t} & = & \D \frac{(\Delta q)^{cl}}{2 \Delta t} - {\delta}^{cl} A_q \\ && \\ \D \frac{\partial T}{\partial t} & = & \D \frac{(\Delta T)^{cl}}{2 \Delta t} \end{array} \end{equation} where ${\delta}^{cl}$ is specified by \begin{equation} {\delta}^{cl} = \left\{ \begin{array}{ll} 1 & \mbox{if } \; \sigma_{Top} \le \sigma \le \sigma_{Lift} \\ & \\ 0 & \mbox{otherwise} \end{array} \right. \end{equation} and $2\Delta t$ is the leap frog time step of the model. The convective precipitation rate $P_{c}$ [m/s] of each cloud layer is \begin{equation} P_{c} = \frac{c_p\, \Delta p}{L\, g \, \rho_{H_2O}} \frac{(\Delta T)^{cl}}{2\Delta t} \end{equation} where $\Delta p$ is the pressure thickness of the layer and $\rho_{H_2O}$ is the density of water. $(\Delta T)^{cl}$ is computed from Eq.~\ref{handm} or Eq.~\ref{handm2}, respectively. \subsection{Shallow Convection} In addition to deep convection a shallow convection scheme is included. Following Tiedtke (1983) shallow convection is parameterized by means of a vertical diffusion of moisture and potential temperature (and, optional, momentum). It is only applied, when the penetrative convection is not operating due to the lack of moisture or (optional) if the unstable layer is below a given threshold height(default is 700hPa). The numerical scheme is similar to that of the normal vertical diffusion (see section \ref{vdiff} but with a constant diffusion coefficient $K$ which is set to default of 10 m$^2$/s within the cloud layer and \begin{equation} 10 \cdot \frac{rh_k -0.8}{1-0.8}(rh_{k+1}-rh_k) \end{equation} at cloud top (here $rh_k$ and $rh_{k+1}$ denote the relative humidity at level above the cloud and the uppermost cloud level, respectively).$K=0.$ elsewhere. The diffusion is limited to the lower part of the atmosphere up to a given pressure (set to a default of 700hPa). For the five level version, the shallow convection is switched off. \subsection{Large Scale Precipitation} Large scale condensation occurs if the air is supersaturated ($q > q_{sat}(T)$). Condensed water falls out instantaneously as precipitation. No storage of water in clouds is considered. An iterative procedure is used to compute final values ($T^*$, $q^*$) starting from the supersaturated state ($T$, $q$): \begin{equation}\label{lsp1} \begin{array}{rcl} \D T^* & = & \D T + \frac{L}{c_p} \, \frac{q - q_{sat} (T)}{1+\frac{L}{c_p}\, \frac{dq_{sat}(T)}{dT}} \\ &&\\ \D q^* & = & \D q -\frac{q - q_{sat}(T)}{1+\frac{L}{c_p}\, \frac{dq_{sat}(T)}{dT}} \end{array} \end{equation} where the suturation specific humidity $q_{sat}$ and its derivative with respect to temperature $dq_{sat}/dT$ are computed from Eqs.~\ref{qsat}. $L$ is either the latent heat of vapourisation or the latent heat of sublimation depending on the temperature. $c_p$ is the specific heat for moist air at constant pressure ($c_p= c_{pd} \, [1+(c_{pv}/c_{pd}-1)\, q]$ where $c_{pd}$ and $c_{pv}$ are the specific heats at constant pressure for dry air and water vapor, respectively). This calculation is repeated once using ($T^*$, $q^*$) as the new initial state. Finally, The temperature and moisture tendencies and the precipitation rate $P_{l}$ [m/s] are computed: \begin{equation} \begin{array}{rcl} \D \frac{\partial T}{\partial t} & = & \D \frac{T^*-T}{2\Delta t} \\ &&\\ \D \frac{\partial q}{\partial t} & = & \D \frac{q^*-q}{2\Delta t} \\ && \\ \D P_{l} & = & \D \frac{p_S \, \Delta \sigma}{g \, {\rho}_{H_2O}} \frac{ (q-q^*)}{2\Delta t} \end{array} \end{equation} where $p_S$ is the surface pressure, $\rho_{H_2O}$ is the density of water, $\Delta \sigma$ is the layer thickness and $2\Delta t$ is the leap frog time step of the model. \subsection{Cloud Formation} Cloud cover and cloud liquid water content are diagnostic quantities. The fractional cloud cover of a grid box, $cc$, is parameterized following the ideas of Slingo and Slingo (1991) using the relative humidity for the stratiform cloud amount $cc_s$ and the convective precipitation rate $P_{c}$ [mm/d] for the convective cloud amount $cc_c$. The latter is given by \begin{equation} cc_c= 0.245 + 0.125 \ln{(P_c)} \end{equation} where $0.05 \le cc_c \le 0.8$. Before computing the amount of stratiform clouds, the relative humidity $rh$ is multiplied by $(1-cc_c)$ to account for the fraction of the grid box covered by convective clouds. If $cc_c \ge 0.3$ and the cloud top is higher than $\sigma = 0.4$ ($\sigma=p/p_S)$, anvil cirrus is present and the cloud amount is \begin{equation} cc_s=2\; (cc_c-0.3) \end{equation} High-, middle- and low-level stratiform cloud amounts are computed from \begin{equation} cc_s=f_{\omega} \left(\frac{rh-rh_c}{1-rh_c}\right)^2 \end{equation} where $rh_c$ is a level depending critical relative humidity. Optionally, a restriction of low-level stratiform cloud amount due to subsidence can by introduced by the factor $f_{\omega}$ where $f_{\omega}$ is depends on the vertical velocity $\omega$. In the default version, $f_{\omega}$~=~1. Cloud liquid water content $q_{H_2O}$ [kg/kg] is computed according to Kiehl et al. (1996): \begin{equation} q_{H_2O} = \frac{q^0_{H_2O}}{\rho} \exp{(-z/h_l)} \end{equation} where the reference value $q^0_{H_2O}$ is $0.21\cdot 10^{-3}$~kg/m$^3$, $\rho$ is the air density, $z$ is the height and the local cloud water scale height $h_l$~[m] is given by vertically integrated water vapor (precipitable water) \begin{equation} h_l= 700 \ln{\left(1 + \frac{1}{g} \int\limits^{p_s}_0 q dp \; \right)} \end{equation} \subsection{Evaporation of Precipitation and Snow Fall} Possible phase changes of convective or large scale precipitation within the atmosphere are considered by melting or freezing of the precipitation depending on the respective level temperature (using 273.16K as a threshold), and by evaporation parameterized in terms of the saturation deficit according to \begin{equation} E_0=-\frac{1}{\Delta t}\cdot \frac{\gamma\cdot(q_0-q_s)}{1+\frac{L\cdot dq_s/dT}{C_{pd}(1+(\delta-1)q_v)}} \end{equation} $\gamma$ is set to a default of 0.01 for T21L10 (0.006 T21L5, 0.007 T42L10). \subsection{Dry Convective Adjustment} Dry convective adjustment is performed for layers which are dry adiabatically unstable, e.g. $\partial \theta / \partial p > 0$ where $\theta$ denotes the potential temperature. The adjustment is done so that the total sensible heat of the respective column is conserved. Wherever dry convection occurs, it is assumed that the moisture is completely mixed by the convective process as well. The adjustment is done iteratively. The atmospheric column is scanned for unstable regions. A new neutral stable state for the unstable region is computed which consists of a potential temperature $\theta_N$ and specific humidity $q_N$: \begin{equation} \begin{array}{lcl} \D \theta_N & = & \D \frac{\sum\limits_{l=l_1}^{l_2} T_l \; \Delta\sigma_l} {\sum\limits_{l=l_1}^{l_2} \sigma_l^{\kappa} \; \Delta\sigma_l} \\ &&\\ \D q_N & = & \D \frac{\sum\limits_{l=l_1}^{l_2} q_l \; \Delta\sigma_l} {\sum\limits_{l=l_1}^{l_2} \Delta\sigma_l} \\ \end{array} \end{equation} where $l_1$ and $l_2$ define the unstable region, $\sigma = (p/p_S)$ is the vertical coordinate, $T$ and $q$ are temperature and specific humidity, respectively, and $\kappa$ is $R_d$/$c_{pd}$ where $R_d$ and $c_{pd}$ are the gas constant and the specific heat for dry air, respectively. The procedure is repeated starting from the new potential temperatures und moistures until all unstable regions are removed. The temperature and moisture tendencies which enter the diabatic time steps are then computed from the final $\theta_N$ and $q_N$ \begin{equation} \begin{array}{lcl} \D \frac{T_l^{t+\Delta t}-T_l^{t-\Delta t}}{2\Delta t} & = & \D \frac{\theta_N \; \sigma_l^{\kappa} -T_l^{t-\Delta t}}{2\Delta t} \\ &&\\ \D \frac{q_l^{t+\Delta t}-q_l^{t-\Delta t}}{2\Delta t} & = & \D \frac{q_N - q_l^{t-\Delta t}}{2\Delta t} \end{array} \end{equation} \newpage \section{Land Surface and Soil \label{landmod}} The parameterizations for the land surface and the soil include the calculation of temperatures for the surface and the soil, a soil hydrology and a river transport scheme. In addition, surface properties like the albedo, the roughness length or the evaporation efficiency are provided. As, at the moment, coupling to an extra glacier module is not available, glaciers are treated like other land points, but with surface and soil properties appropriate for ice. Optionally, A simple biome model can be used (simba). \subsection{Temperatures}\label{surtemp} The surface temperature $T_S$ is computed from the linearized energy balance of the uppermost $z_{top}$ meters of the ground: \begin{equation} \label{land.1} c_{top} \; z_{top} \; \frac{\Delta T_S}{\Delta t} = F_S - G + \Delta T_S \; \frac{\partial (Q_a -F_g)}{\partial T_S} - F_m \end{equation} $z_{top}$ is a prescribed parameter and set to a default value of $z_{top}$~=~0.20~m. $Q_a$ denotes the total heat flux from the atmosphere, which consists of the sensible heat flux, the latent heat flux, the net short wave radiation and the net long wave radiation. $Q_g$ is the flux into the deep soil. $Q_a$ and $Q_g$ are defined positive downwards. $Q_m$ is the snow melt heat flux and $c_{top}$ is the volumetric heat capacity. Depending on the snow pack, $z_{top}$ can partly or totally consist of snow or soil solids: $z_{top}=z_{snow}+z_{soil}$. Thus, the heat capacity $c_{top}$ is a combination of snow and soil heat capacities: \begin{equation} c_{top} =\frac{ c_{snow} \; c_{soil} \; z_{top}}{c_{snow} \; z_{soil} + c_{soil} \; z_{snow}} \end{equation} The default value of $c_{snow}$ is 0.6897~$\cdot$~10$^{6}$~J/(kg K) using a snow density of 330~kg/m$^3$. $c_{soil}$ is set to a default value of 2.07~$\cdot$~10$^{6}$~J/(kg K) for glaciers and to a value of 2.4~$\cdot$~10$^{6}$~J/(kgK) otherweise. Below $z_{top}$ the soil column is discretized into $N$ layers with thickness $\Delta z_i$, where layer $1$ is the uppermost of the soil layers. The default values for the model are $N$~=~5 and $\Delta z$~=~(0.4~m, 0.8~m, 1.6~m, 3.2~m, 6.4~m). The heat flux into layer 1, $Q_g$, is given by \begin{equation} Q_g=\frac{2 k_{1}}{\Delta z_{1}} (T_S - T_{1}) \end{equation} where $k_{1}$ and $T_{1}$ are the thermal conductivity and the temperature. If the snow depth is greater than $z_{top}$, the thermal properties of snow are blended with the first soil layer to create a snow/soil layer with thickness $z_{snow}-z_{top}+\Delta z_1$. The thermal conductivity $k_1$ and heat capacity $c_1$ of a snow/soil layer are \begin{equation} \begin{array}{rcl} \D k_1 & = & \D \frac{k_{snow}\; k_{soil}\; (\Delta z_1+z_{snow}-z_{top})}{k_{snow}\; \Delta z_1 + k_{soil}\; (z_{snow}-z_{top})} \\ &&\\ \D c_1 & = & \D \frac{c_{snow}\; c_{soil}\; (\Delta z_1+z_{snow}-z_{top})}{c_{snow}\; \Delta z_1 + c_{soil}\; (z_{snow}-z_{top})} \end{array} \end{equation} with default values of $k_{snow}$~=~0.31~W/(m K), $k_{soil}$~=~2.03~W/(m K) for glaciers and $k_{soil}$~=~7~W/(m K) otherweise. After the surface temperature $T_S$ has been calculated from Eq.~\ref{land.1}, snow melts when $T_S$ is greater than the freezing temperature $T_{melt}$. In this case, $T_S$ is set to $T_{melt}$ and a new atmospheric heat flux $Q_a(T_{melt})$ is calculated from $Q_a$ and $\partial Q_a/\partial T_S$. If the energy inbalance is positive ($Q_a(T_{melt}) > c_{top} \; z_{top}\; (T_{melt} - T_S^{t})/\Delta t$; where $T_S^{t}$ is the surface temperature at the previous time step), the snow melt heat flux $Q_m$ is \begin{equation}\label{melt} Q_m=\max(Q_a(T_{melt}) - \frac{c_{top} \; z_{top}}{\Delta t} \; (T_{melt} - T_S^{t}), \; \frac{W_{snow} \; L_f}{\Delta t}) \end{equation} where $W_{snow}$ is the mass of the snow water of the total snow pack and $L_f$ is the latent heat of fusion. Any excess of energy is used to warm the soil. With the heat flux $F_z$ at depth $z$ of the soil \begin{equation} F_z = -k \; \frac{\partial T}{\partial z} \end{equation} one dimensional energy conservation requires \begin{equation} c\; \frac{\partial T} {\partial t} = - \frac{\partial F_z}{\partial z}= \frac{\partial}{\partial z} \left [ k \; \frac{\partial T}{\partial z} \right] \end{equation} where $c$ is the volumetric soil heat capacity, $T$ is the soil temperature, and $k$ is the thermal conductivity. In the model, thermal properties (temperature, thermal conductivity, volumetric heat capacity) are defined at the center of each layer. Assuming the heat flux from $i$ to the interface $i$ and $i+1$ equals the heat flux from the interface to $i+1$, the heat flux $F_i$ from layer $i$ to layer $i+1$ (positive downwards) is given by \begin{equation} F_i= - \frac{2\; k_i\; k_{i+1} (T_i - T_{i+1})}{k_{i+1} \; \Delta z_i + k_i\; \Delta z_{i+1}} \end{equation} The energy balance for layer $i$ is \begin{equation} \frac{c_i\; \Delta z_i}{\Delta t} \; (T_i^{t+\Delta t} - T_i ^{t}) = F_i - F_{i-1} \end{equation} The boundary conditions are zero flux at the bottom of the soil column and heat flux $F_g$ at the top. This equation is solved implicitly using fluxes $F_i$ evaluated at $t+\Delta t$ \begin{equation} \begin{array}{lclcl} \D \frac{c_i \Delta z_i}{\Delta t} (T_i^{t+\Delta t} - T_i ^{t}) & = &\D \frac{k_i k_{i+1} (T_{i+1}^{t+\Delta t} - T_i^{t+\Delta t})}{k_{i+1} \Delta z_i + k_i \Delta z_{i+1}} + G & for & \D i = 1 \\ & & & & \\ \D \frac{c_i \Delta z_i}{\Delta t} (T_i^{t+\Delta t} - T_i ^{t}) & = &\D \frac{k_i k_{i+1} (T_{i+1}^{t+\Delta t} - T_i^{t+\Delta t})}{k_{i+1} \Delta z_i + k_i \Delta z_{i+1}} + \frac{k_i k_{i-1} (T_{i- 1}^{t+\Delta t} - T_i^{t+\Delta t})}{k_{i-1} \Delta z_i + k_i \Delta z_{i-1}} & for & \D 1 < i < N \\ & & & & \\ \D \frac{c_i \Delta z_i}{\Delta t} (T_i^{t+\Delta t} - T_i ^{t}) & = &\D \frac{k_i k_{i-1} (T_{i- 1}^{t+\Delta t} - T_i^{t+\Delta t})}{k_{i-1} \Delta z_i + k_i \Delta z_{i-1}} & for & \D i = N \end{array} \end{equation} resulting in a linear system for the $T_i^{t+\Delta t}$. \subsection{Soil Hydrology}\label{hydro} The parameterization of soil hydrology comprises the budgets for snow amount and the soil water amount. The water equivalent of the snow layer $z_{snow}^{H_2O}$ is computed over land and glacier areas from \begin{equation} \frac{\partial z_{snow}^{H_2O}}{\partial t} = F_q + P_{snow}-M_{snow} \end{equation} where $F_q$ is the evaporation rate over snow computed from Eq.~\ref{fluxes2}, $P_{snow}$ is the snow fall and $M_{snow}$ is the snow melt rate (all fluxes are positive downward and in m/s). $M_{snow}$ is related to the snow melt heat flux $Q_m$ (Eq.~\ref{melt}) by \begin{equation} M_{snow}=\frac{Q_m}{\rho_{H_2O}\, L_f} \end{equation} where $L_f$ is the latent heat of fusion. The soil water reservoir $W_{soil}$ [m] is represented by a single-layer bucket model (Manabe 1969). Soil water is increased by precipitation $P$ and snow melt $M_{snow}$ and is depleted by the surface evaporation $F_q$: \begin{equation} \frac{\partial W_{soil}}{\partial t} = P + M + F_q \end{equation} where all fluxes are defined positive downwards and in m/s. Soil water is limited by a field capacity $W_{max}$ which geographical distribution can be prescribed via an external input or is set to a default value of 0.5~m everywhere. If the soil water exceeds $W_{max}$ the excessive water builds the runoff $R$ and is provided to the river transport scheme (Section~\ref{runoff}). The ratio of the soil water and the field capacity defines the wetness factor $C_w$ which is used in Eq.~\ref{fluxes2} to compute the surface evaporation: \begin{equation}\label{cwgl} C_w=\frac{W_{soil}}{f_{Cw} \; W_{max}} \end{equation} where the factor $f_{Cw}$ (with a default value of 0.25) takes into account that maximum evaporation will take place even if the bucket is not completely filled. For land points covered by glaciers, $C_w$ is set to a constant value of 1. \subsection{River Transport}\label{runoff} The local runoff is transported to the ocean by a river transport scheme with linear advection (Sausen et al.~1994). For each grid box (both, land and ocean costal points) the river water amount $W_{river}$ [m$^3$] is computed from \begin{equation} \frac{\partial W_{river}}{\partial t}= ADV + area \; (R - S) \end{equation} where $R$ is the local runoff (Section~\ref{hydro}), $S$ is the input into the ocean, $ADV$ is the advection of river water and $area$ is the area of the respective grid box. The input into the ocean $S$ is given by \begin{equation} S=\left\{ \begin{array}{ll} 0 & \mbox{for land points} \\ &\\ ADV & \mbox{for ocean points} \end{array} \right. \end{equation} This ensures that $S$ is non-zero only for ocean costal points. The advection from grid box $(i,j)$ into grid box $(i',j')$, $ADV_{(i,j)\rightarrow (i',j')}$, is formulated using an upstream scheme: \begin{equation} \begin{array}{rcl} \D ADV_{(i,j)\rightarrow (i+1,j)} & = & \D \left\{ \begin{array}{ll} u_{i,j} W_{i,j}, & \mbox{if } \; u_{i,j} \ge 0 \\ & \\ u_{i,j} W_{i+1,j}, & \mbox{if } \; u_{i,j} < 0 \end{array} \right. \\ && \\ \D ADV_{(i,j)\rightarrow (i,j+1)} & = & \D \left\{ \begin{array}{ll} -v_{i,j} W_{i,j}, & \mbox{if } \; v_{i,j} \le 0 \\ & \\ -v_{i,j} W_{i,j+1}, & \mbox{if } \; v_{i,j} > 0 \end{array} \right. \\ \end{array} \end{equation} where $i$ and $j$ are the zonal and meridional indices of the grid box, which are counted from the west to the east and from the north to the south, respectively. The zonal and meridional advection rates $u_{i,j}$ and $v_{i,j}$ are defined at the interface of two grid boxes and depend on the slope of the orography: \begin{equation} \begin{array}{rcl} \D u_{i,j} & = & \D \frac{c}{\Delta x}\left[\frac{h_{i,j}-h_{i+1,j}}{\Delta x}\right]^{\alpha} \\ && \\ \D v_{i,j} & = & \D \frac{c}{\Delta y}\left[\frac{h_{i,j+1}-h{i,j}}{\Delta y}\right]^{\alpha} \end{array} \end{equation} where $\Delta x$ and $\Delta y$ are the distances between the grid points in the longitudinal and the meridional direction. $h$ is the height of the orography, which is modified in order to omit local minima at land grid points. The empirical constants $c$ and $\alpha$ are set to the values given by Sausen et al.~(1994) for T21 resolution ($c = $~4.2~m/s and $\alpha =$~0.18). \subsection{Other Land Surface Parameter}\label{landsurf} Some additional quantities characterizing the land surface of each grid box need to be specified for use in the model. The land-sea mask and the orography are read from an external file. Optionally, this file may also include other climatological surface parameter: the global distribution of the surface roughness length $z_0$, a background albedo ${\cal R}_S^{clim}$, a glacier mask for permanent ice sheets, the bucked size for the soil water $W_{max}$ (see section above) and a climatological annual cycle of the soil wetness $C^{clim}_w$ (which may be used instead of the computed $C_w$ from Eq.~\ref{cwgl}. If there is no input for the particular field in the file, the parameter is set to be horizontal homogeneous with a specific value. The following defaults are used: $z_0$~=~ 2~m, ${\cal R}_S^{clim}$~=~0.2, no glaciers, $W_{max}$~=~0.5 and $C^{clim}_w$~=~0.25. For snow covered areas, the background albedo is modified to give the actual albedo ${\cal R}_S$ which is used in the radiation scheme. For points, which are not covered by glaciers, ${\cal R}_S$ is given by \begin{equation} {\cal R}_S={\cal R}_S^{clim} + ({\cal R}_S^{snow}-{\cal R}_S^{clim}) \; \frac{z_{snow}}{z_{snow} + 0.01} \end{equation} where $z_{snow}$ is the snow depth, and the albedo of the snow, ${\cal R}_S^{snow}$, depends on the surface temperature $T_S$ \begin{equation}\label{rsnow} {\cal R}_S^{snow}={\cal R}_{max}^{snow} + ({\cal R}_{min}^{snow} - {\cal R}_{max}^{snow}) \; \frac{T_S - 263.16}{10} \end{equation} with ${\cal R}_{min}^{snow} \le {\cal R}_S^{snow} \le {\cal R}_{max}^{snow}$ and default values ${\cal R}_{min}^{snow}$~=~0.4 and ${\cal R}_{max}^{snow}$~=~0.8. For glaciers, ${\cal R}_S$ is given by ${\cal R}_S^{snow}$ from Eq.~\ref{rsnow} but with a default ${cal R}_{min}^{snow}$~=~0.6. The surface specific humidity $q_S$ is given by the saturation specific humidity at $T_S$: \begin{equation} q_S =q_{sat}(T_S) \end{equation} where $q_{sat}(T_S)$ is computed from Eq.~\ref{qsat}. \section{Sea Surface}\label{seasurf} Sea surface temperatures $T_{sea}$, sea ice distributions $c_{ice}$ and surface temperatures over sea ice $T_i$ are provided by the ocean and sea ice modules (Section HEIKO). From these quantities, the following additional parameter are computed which enter the atmospheric parameterizations. The prescribed surface albedo ${\cal R_S}$ for open water is set to a default value of 0.069. For sea ice ${\cal R}_S$ is given as a function of the ice surface temperature $T_{i}$: \begin{equation} {\cal R}_S=\min{({\cal R}_S^{max}, \, 0.5 + 0.025 \, (273. - T_{i}))} \end{equation} where the prescribed maximum sea ice background albedo ${\cal R}_S^{max}$ is set to a default value of 0.7. The surface specific humidity $q_S$ is given by the saturation specific humidity at the surface temperature $T_S$ which is either $T_{sea}$ or $T_{i}$: \begin{equation} q_S =q_{sat}(T_S) \end{equation} where $q_{sat}(T_S)$ is computed from Eq.~\ref{qsat}. The wetness factor $C_w$ which enters the calculation of the surface evaporation (Eq.~\ref{fluxes2}) is set to 1. The roughness length $z_0$ over sea ice is set to a constant value of $z_0$~=~0.001~m. Over open water, $z_0$ is computed from the Charnock (1955) formula: \begin{equation} z_0 = C_{char} \frac{u_{*}^{2}}{g} \end{equation} with a minimum value of $1.5 \cdot 10^{-5}$~m. $C_{char}$ denotes the Charnock constant and is set to 0.018. $g$ is the gravitational acceleration. The friction velocity $u_{*}$ is calculated from the surface wind stress at the previous time level: \begin{equation} u_{*}=\sqrt{\frac{|F_u, F_v|}{\rho} } \end{equation} where $|F_u, F_v| $ is the absolute value of the surface wind stress computed from Eq.~\ref{fluxes2} and $\rho$ is the density. \newpage \section{References} Berger, A., 1978a: A simple algorithm to compute long-term variations of daily insolation. Institute of Astronomy and Geophysic, Universite Catholique de Louvain, {\bf Contribution 18}, Louvain-la-Neuve, Belgium. Berger, A., 1978b: Long-term variations of daily insolation and quaternary climatic change. {\it J.~Atmos.~Sci.}, {\bf 35}, 2362- 2367. Blackadar, A.~K., 1962: The vertical distributionof wind and turbulent exchange in a neutral atmosphere. {\it J.~Geophys.~Res.}, {\bf 67}, 3095-3102. Boer, G.~J., N.~A.~McFarlane, R.~Laprise, J.~Henderson and J.-P.~Blanchet, 1984: The Canadian Climate Centre spectral atmospheric general circulation model. {\it Atmosphere-Ocean}, {\bf 22}, 397-429. Charnock, M., 1955: Wind stress on a water surface. {\it Q.~J.~R.~Meteorol.~Soc.}, {\bf 81}, 639-640. Chou, M.-D., M.~J.~Suarez, X.-Z.~Liang and M.~M.~H.~Yan, 2001 (Revised 2002): A thermal infrared radiation parameterization for atmospheric studies. {\it Technical Report Series on Global Modelling and Data Assimilation, M.~J.~Suarez Ed.}, {\bf NASA/TM-2001-104606, Vol. 19}, 55pp. Green, A.~E.~S., 1964: Attenuation by ozone and the earth's albedo in the middle ultraviolet. {\it AAppl.~Opt.}, {\bf 3}, 203-208. Katayama, A., 1972: A simplified scheme for computing radiative transfer in the troposphere. Department of Meteorology, {\it Tech.~Report}, {\bf No. 6}, University of California, Los Angeles, CA, 77 pp. Kiehl, J.~T., J.~J.~Hack, G.~B.~Bonan, B.~A.~Boville, B.~P.~Briegleb, D.~L.~Williamson and P.~J.~Rasch, 1996: Description of the NCAR Community Climate Model (CCM3). {\it NCAR Technical Note}, {\bf NCAR/TN-420+STR}, 152pp. Kuo, H.~L., 1965: On formation and intensification of tropical cyclones through latent heat release by cumulus convection. {\it J.~Atmos.~Sci.}, {\bf 22}, 40-63. Kuo, H.~L., 1974: Further studies of the parameterization of the influence of cumulus convection on large-scale flow. {\it J.~Atmos.~Sci.}, {\bf 31}, 1232-1240. Lacis, A.~A., and J.~E.~Hansen, 1974: A parameterization for the absorption of solar radiation in the Earth's atmosphere. {\it J.~Atmos.~Sci.}, {\bf 31}, 118-133. Laurson, L.~and E.~Eliasen, 1989: On the effects of the damping mechanisms in an atmospheric general circulation model. {\it Tellus}, {\bf 41A}, 385-400. Louis, J.~F., 1979: A parametric model of vertical eddy fluxes in the atmosphere. {\it Boundary Layer Meteorology}, {\bf 17}, 187-202. Louis, J.~F., M.~Tiedtke and M.~Geleyn, 1982: A short history of the PBL parameterisation at ECMWF. Proceedings, ECMWF workshop on planetary boundary layer parameterization, Reading, 25-27 Nov. 81, 59-80. Lowe, P.~R., 1977: An approximating polynomial for the computation of saturation vapour pressure. {\it J.~Appl.~Met.}, {\bf 16}, 100-103. Manabe, S.~and F.~M\"oller, 1961: on the radiative equilibrium and heat balance of the atmosphere. {\it Mon~Wea.~Rev.}, {\bf 89}, 503-532. Manabe, S., 1969: Climate and ocean circulation. I. The atmospheric circulation and the hydrology of the earth's surface. {\it Mon.~Wea.~Rev.}, {\bf 97}, 739-774. Miller, M.~J., A.~Beljaars and T.~N.~Palmer, 1992: The sensitivity of the ECMWF model to the parameterization of evaporation from the tropical oceans. {\it J.~Climate}, {\bf 5}, 418-434. Rodgers, C.~D., 1967: The use of emissivity in the atmospheric radiation calculation. {\it Q.~J.~R.~Meteorol.~Soc.}, {\bf 93}, 43-54. Roeckner, E., K.~Arpe, L.~Bengtsson, S.~Brinkop, L.~D\"umenil, M.~Esch, E.~Kirk, F.~Lunkeit, M.~Ponater, B.~Rockel, R.~Sausen, U.~Schlese, S.~Schubert and M.~Windelband, 1992: Simulation of the present-day climate with the ECHAM-3 model: Impact of model physics and resolution. Max-Planck Institut f\"ur Meteorologie, {\bf Report No.~93}, 171pp. Sasamori, T., 1968: The radiative cooling calculation for application to general circulation experiments. {\it J.~Appl.~Met.}, {\bf 7}, 721-729. Sausen, R., S.~Schubert and L.~D\"umenil, 1994: A model of river runoff for use in coupled atmosphere-ocean models. {\it Journal of Hydrology}, {\bf 155}, 337-352. Slingo, A., and J.~M.~Slingo, 1991: Response of the National Center for Atmospheric Research community climate model to improvements in the representation of clouds. {\it J.~Geoph.~Res.}, {\bf 96}, 341-357. Stephens, G.~L., 1978: Radiation profiles in extended water clouds. II: Parameterization schemes. {\it J.~Atmos.~Sci.}, {\bf 35}, 2123-2132. Stephens, G.~L., 1984: The parameterization of radiation for numerical weather prediction and climate models. {\it Mon.~Wea.~Rev.}, {\bf 112}, 826-867. Stephens, G.~L., S.~Ackermann and E.~A.~Smith, 1984: A shortwave parameterization revised to improve cloud absorption. {\it J.~Atmos.~Sci.}, {\bf 41}, 687-690. Tiedtke, M., 1983: The sensitivity of the time-mean large-scale flow to cumulus convection in the ECMWF model. Proceedings of the ECMWF Workshop on Convection in {\it Large-Scale Models, 28 November-1 December 1983, European Centre for Medium-Range Weather Forecasts, Reading, England}, 297-316.