% 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.