\documentclass{article}
\usepackage{a4wide}
\usepackage{graphicx}
\textheight=23cm
\topmargin -1cm
%
%\input{title}
\begin{document}

\bibliographystyle{apalike}
\newcommand{\dt}[1]{\frac{\partial #1 }{\partial t} }
\newcommand{\av}[2]{{\overline{#1}}^{ #2 }}
\newcommand{\avt}[1]{{\overline{#1}}}
\newcommand{\avg}[1]{\left< #1 \right>}
% \newcommand{\emis}[2]{E \left( #1 , #2 \right) }
\newcommand{\der}[2]{\frac{\partial #1 }{\partial #2} }
\newcommand{\dep}[1]{\left( #1 \right) }
\newcommand{\depb}[1]{\left[ #1 \right] }
\newcommand{\resol}[3]{ {\em #1 longitudes $\times$  #2 latitudes $\times$ #3  vertical levels} }
\def\t{\theta}
\def\cd{C_D}

% cross-references
\def\Sc#1{Section~\ref{sc:#1}}
\def\an#1{Appendix~\ref{an:#1}}
\def\ch#1{Chapter~\ref{ch:#1}}
\def\fig#1{Fig.~\ref{fg:#1}}
\def\fg#1{Fig.~\ref{fg:#1}}
\def\figs#1{Figs.~\ref{fg:#1}}
\def\eq#1{Eq.~\ref{eq:#1}}
\def\eqs#1{Eqs.~\ref{eq:#1}}
\def\tb#1{Table~\ref{tb:#1}}

\def\gcm{general circulation model}
\def\adhoc{{\em ad hoc}}
\def\ie{{\em i.~e.}}
\def\apost{{\em a posteriori}}
\def\apriori{{\em a priori}}
\def\eg{{\em e.~g.}}
\def\etal{{\em et~al.}}

\def\pun{Paper~1}






\title{The 20-parameter general circulation model
\author{
Fr\'ed\'eric Hourdin\\
hourdin@lmd.ens.fr
}\date{\today}}


\bibliographystyle{apalike}

\maketitle





\section{The planetary atmospheric general circulation model\label{sc:gcm}}

Schematically, atmospheric general circulation models
consist of two parts. (1) The "dynamical part" is dedicated to
a temporal and spatial integration of the 
hydrodynamical equations for the explicitly resolved scales (generally
of the order from few tens to few hundreds km on the horizontal
for classical terrestrial applications).
Terrestrial models are generally based on the "primitive equations of
meteorology" which are valid for Earth-like planets (solid
body covered by a thin atmosphere assumed to remain vertically in
hydrostatic balance). This part of the code can thus be
used for various planets by just changing some parameters such as the
gravity or rotation rate or some boundary conditions such as the
specification of the surface height. (2) The "physical part" is dedicated
to the diabatic forcing of the atmospheric circulation and to the
modeling of sub-grid scale processes.
The corresponding parametrizations are generally much
less universal and must be significantly modified from one planet to the
other. For example, the complete intergation of atmospheric radiative transfer
(which includes absorption by complex molecular vibration-rotation bands
with complex temperature and pressure dependencies) is by far out of
reach of present computational means, since one
wants to integrate those equations several time per day and for
thousands various locations on the planet. Numerous
approximations are therefore used in general circulation models which are
generally strongly dependent on the particular planetary conditions.
In fact the three planetary versions of the Planetary Atmospheric
General Circulation Model developed at LMD essentially differ by the way
these radiative transfer calculations are performed.

\subsection{Dynamics}

The dynamics of the planetary model are taken from the terrestrial
general circulation model used at LMD for climate studies
Classical {\em primitive equations} of meteorology (Navier Stokes equations
for a perfect gas plus hydrostatic hypothesis and thin layer approximation,
see \eg\ Holton 1979~\nocite{Holt:79})
are integrated using finite differences in space and an
explicit centered (or leapfrog) temporal scheme. The original formulation
of this model is presented in details \nocite{Sado:84} by Sadourny and
Laval, 1984. We use in fact a
more recent formulation by Sadourny, coded in FORTRAN~77 by Le~Van,
in which latitudes and longitudes of
the grid can be chosen arbitrarily, which allows to zoom in a particular
region on the sphere.
For our experiments,  we chose a uniform area grid with points 
equidistant in longitude $\lambda$ and sine of latitude $\phi$ so that
all meshes have the same area.
Two horizontal resolutions were used, one with 32 longitudes and 24
latitudes and the other one with 64 longitudes and 48 latitudes.
The vertical discretization is based on a hybrid $\sigma-p$ coordinate ($\sigma$ beeing pressure
normalized by its surface value $\sigma=p/p_s$). The $\sigma-levels$ are
regularly distributed in $\ln{\sigma}$ except near the lower
and upper boundaries.

The discretization scheme conserves both potential enstrophy for barotropic
divergent flows \cite{Sado:75,Sado:75a} and total angular momentum for
stationary axi-symmetric flows.
The latter property was not included in the original terrestrial version but
was introduced in the Martian  atmospheric circulation model
\cite{Hour:93} in which, because of the non numerical conservation of
angular momentum, spurious prograde zonal winds were produced
in equatorial regions.
The new conservative formulation has 
since then been included in the terrestrial climate model
and was found to strongly improve the
simulation of the tropical trade-winds.
Accurate conservation of angular momentum is, of
course, crucial for the study of atmospheric superrotation
since the budget of this particular quantity plays a central role in
the phenomenon.
The numerical scheme and conservation properties are presented in details in
\an{dyn}.

Non linear interactions between explicitly resolved scales
and subgrid-scale processes are
parametrized applying an horizontal dissipation operator on potential
temperature and wind.
This dissipation is based on an iterated Laplacian:
\begin{equation}
\frac{\partial q}{\partial t} =
(-1)^{n_{diss}} \frac{1}{\tau_{diss}}
{(\delta x)}^{2n_{diss}} \triangle^{n_{diss}} q
\end{equation}
where $\delta x$ is the smallest horizontal distance represented in the
model and $\tau_{diss}$ the dissipation time of a structure of scale $\delta x$,
a technique commonly used in the modeling of two-dimensional turbulence
and in many general circulation models.

The iteration is introduced to be more selective in smaller
scales, not acting directly on the explicitly resolved flow (in the
present simulations, $n_{diss}=2$ and $\tau_{diss}=6\times 10^4$~s). Too
strong a dissipation could produce some efficient process to transport
angular momentum equatorward which would artificially reinforce superrotation.
This was responsible for some peculiar results in the past: a strong
superrotation was produced in some terrestrial simulations~\cite{Hunt:73}
which were due to a strong latitudinal filter on the zonal wind, acting
as a very efficient horizontal dissipation.

\subsection{Physical parametrizations}

The main problem for the determination of an appropriate set of physical
parametrizations was the development of
a realistic but simple radiative scheme.
Many fundamental numerical studies of the atmospheric general
circulation are based on the use of a so-called
Newtonian forcing in which the temperature field is simply relaxed
toward a prescribed "radiative equilibrium"
state~\cite{Held:80,Lind:88,Jame:92}.
This approach has the advantage of simplicity and is characterized by a
time constant which can easily be compared to the other time constants
of the circulation, but it presents important drawbacks:
it is physically incorrect and
moreover, the
equilibrium temperature field is prescribed and not related to
the relevant external parameters of the planetary environment. In such
experiments the
zonal wind field is also indirectly
relaxed to that in thermal wind balance with the
equilibrium temperature field whereas, in real atmospheres, the coupling
between temperature and wind is also important to determine the
equilibrium temperature field itself.

For the present model, a complete but very simple
radiative transfer code was developed, based on some kind of grey absorption
approximation. The parametrization depends upon only two independent and
external parameters: the total atmospheric transmission for the
visible and thermal parts of the spectrum ($\tau_{IR}$ and $\tau_{Vis}$
respectively).
Finally, the parametrization accounts for the effect of temperature
variations of the Planck function (which is fundamental in order to
obtain a radiative equilibrium) but not for the main complexity 
of atmospheric opacities.
We also introduced (as is usually done in terrestrial models)
parametrizations for the vertical turbulent diffusion,
convective adjustment (in order to prevent supper-adiabatic
vertical temperature profiles) and thermal conduction under the surface.
All those parametrizations are presented
in~\an{phys}.
In the present study, it was not attempted to introduce any kind
of hydrological cycle or
phase change as occur on most known planetary atmospheres with
significant impact on the atmospheric circulation.

\subsection{The planetary parameters}

In this simple planetary atmospheric general circulation model,
the description of the planet only depends upon 19 independent
parameters:
\begin{itemize}
\item 6 astronomical parameters: the radius $a$ and gravity $g$ of the
planet; its rotation rate $\Omega$; the length of year $yr$,
obliquity $\delta$ and solar constant $F_0$ (the orbit is assumed to be
circular).
\item 4 parameters for the surface properties: emissivity $\epsilon_s$,
albedo $A_s$, drag coefficient $C_D$ and thermal inertia $I$, assumed to
be constant over the whole surface of the planet.
\item 3 parameters to describe the atmosphere: mean surface pressure
$p_s$, molecular mass $\mu_{gas}$ (which enters in the
definition of the gas constant $R=\mu R^*$, where $R^*$ is the universal
constant for perfect gas) and the adiabatic
coefficient $\kappa=R/C_p$.
\item the 2 transmissions for the visible and thermal part of
the spectrum ($\tau_{IR}$ and $\tau_{Vis}$).
\item 4 coefficients for the sub-grid scale processes: a mixing length
$l_{mix}$
and a minimum turbulent kinetic energy $e_{mix}$ for the vertical 
mixing and the two coefficients $n_{diss}$
and $\tau_{diss}$ of the horizontal dissipation operator.
\end{itemize}
The choice of this particular set of parameters is of course rather arbitrary.
For
instance, other astronomical parameters must be introduced to study the
effects of strong excentricity as arise on Mars.
Some choices are also strongly model dependent such as those
parameters which are related to sub-grid scale processes.
Finally, for the purpose of general investigation of the possible atmospheric
circulation regimes, it will be useful to reduce those parameters 
through dimensional analysis, which has not been attempted here.

\appendix


\section{Physical parametrizations\label{an:phys}}

The physical parametrizations are presented in details below.

\subsubsection{Radiative transfer}

\def\ir{{_{IR}}}
\def\vis{{_{Vis}}}
\def\upw{^{\uparrow}}
\def\dnw{^{\downarrow}}
%\newcommand{\trans}[2]{ \epsilon \left( #1 , #2 \right) }
\newcommand{\trans}[2]{ \tau_{  #1 , #2 } }
\def\dz{\delta_z}

The time evolution of
the temperature ($T_{k+\frac{1}{2}}$)
at the middle of a layer (or $s$ level),
due to absorption and emission of infra-red radiations, is
computed as
\begin{equation}
\left( \dt{ T_{k+\frac{1}{2}} } \right)\ir =
\frac{g}{C_p}
\frac{ F_k^{\upw}-F_k^{\dnw}  - F_{k+1}^{\upw}+F_{k+1}^{\dnw} }
{ p_k -p_{k+1} }
\end{equation}
where the upward $F_k^{\upw}$ and downward $F_k^{\dnw}$ fluxes  are computed
using a spectrally averaged transmissivity
$\trans{}{}$ in the following form
\begin{equation}
F_k^{\upw}  =  \epsilon_{s} \sigma {T_s}^4  \trans{0}{k}
+ \sum_{l=0}^k { \sigma T_{l+\frac{1}{2}}^4
\left[ \trans{l+1}{k} -\trans{l}{k} \right]  }
\end{equation}
\begin{equation}
F_k^{\dnw} = \sum_{l=k}^K 
{ \sigma T_{l+\frac{1}{2}}^4
\left[ \trans{k}{l} -\trans{k}{l+1} \right]  }
\end{equation}
(where $\epsilon_s$ is the surface emissivity).
Various simple forms can be chosen for the transmissivity.
First tests, performed with a grey absorption approximation, led to strong
temperature discontinuities at the surface (between the ground
temperature and the air above the surface), as already found by other
authors at
the beginning of climate numerical modeling~\cite{Mana:61}.

Finally, the transmissivities were computed in the frame
of the strong absorption approximation, better adapted to the
tropospheric conditions where most of the efficient molecular bands are
strongly saturated. In this case, assuming that lines are randomly
distributed over the spectrum, and neglecting temperature variations
of the absorbing coefficients along the optical path, the transmissivity of the
whole band, between two layers $k$ and $l$, can be computed, assuming a
completely isotropic radiation, in the form (for
those who are not familiar with such models see Goody, 1964):\nocite{Good:64}
%
\def\up{\overline{up}}
\def\u{\overline{u}}
\begin{equation}
\trans{k}{l}=\exp{ \left( - k_\ir\up \right) }
\end{equation}
%
where $\up$ is a pressure weighted absorber amount
%
\begin{equation}
\up = \int { p \rho dz}=\left( p_k -p_l \right) 
\frac{\left( p_k +p_l \right) }{2}
\end{equation}
%
and $k_\ir$ a constant describing the efficiency of the atmospheric absorption.

This formula produces again some temperature discontinuity at the surface
but less than 10~K.
The transmissivity depends on one free parameter only,
for example the total optical transmission of the atmosphere
$\tau_\ir=\exp\dep{-k_\ir p_s^2/2}$ (fraction of the surface thermal radiation
directly emitted to space). The model accounts for the fundamental effect of
temperature variations linked to the Planck function, essential
for the thermal equilibrium and radiative forcing of the
circulation, but not for the much more complex pressure and temperature
dependencies linked to the molecular absorbing coefficients.\\

A classical grey absorption approximation is used for the absorption of
solar radiations.
For reason of simplicity, no attempt was
made to model some kind of scattering by molecules or aerosols.
In such conditions, the downward solar radiation consists of a uniq
collimated beam and the grey transmission takes the simple form
\begin{equation}
\trans{k}{l}=\exp\dep{-\frac{k_\vis}{\cos\theta}\left| p_k - p_l \right| }
\end{equation}
where $\theta$ is the local solar zenith angle.
As for the infra-red radiative transfer, the model depends on one
free parameter, for instance the
total atmospheric transmission $\tau_\vis$, computed for a mean incidence
$\cos\theta=1/2$ (fraction of the global incoming solar
radiation reaching the surface), which is the fraction of the solar flux
reaching the surface of the planet. It influences greatly the vertical stability
of the atmosphere, an important dynamical parameter.

The reflected
radiation, which depends on a surface albedo $A_s$, is assumed to be
totally isotropic in the upward direction and is computed for a
mean zenith angle $\theta_0$ with $1/\cos\theta_0=1.66$.

The radiative forcing of the atmosphere also depends on the incoming solar
flux at the top of the atmosphere $F_0$ and on those of the orbital
parameters which determine the spatial and temporal variations of insolation 
(length of day, obliquity and length of year).

\subsubsection{Vertical turbulent mixing}

The formulation of the vertical turbulent mixing is taken almost
directly from a very simple scheme used in the LMD climate model. The
effect of turbulent mixing on both potential temperature and momentum is
computed as
\begin{equation}
\der{q}{t}=-\frac{1}{\rho}\der{F_q}{z}
\end{equation}
in term of the vertical divergence of a turbulent flux
\begin{equation}
F_q=-\rho K_z \der{q}{z}
\end{equation}
The turbulent mixing
coefficient $K_z$ is computed as
%
\begin{equation}
K_z=l e^{1/2}
\end{equation}
%
in term of mixing length $l$ and a diagnostic estimate of the turbulent
kinetic energy
%
\begin{equation}
e=Max\left\{ l^2
\left[ \dep{ \der{\vec{V}}{z}}^2 - 2.5 \frac{g}{\t} 
\der{\t}{z}  \right] , e_{min} \right\}
\end{equation}
%
where the minimum value of the kinetic energy was set to
$e_{min}=1\times 10^{-6}m^2s^{-2}$.
In the case where $e>e_{min}$, $K_z$,expressed as a function
of the Richardson number
%
\begin{equation}
Ri=\frac{g\partial\t/\partial z}{\t\dep{\partial \vec{V}/\partial
z}^2},
\end{equation}
%
is equal to
%
\begin{equation}
K_z=l^2 \left\| \der{\vec{V}}{z}\right\| \sqrt{1-2.5 Ri}
\end{equation}

The surface flux ${F_q}_s$ is simply computed as
%
\begin{equation}
{F_q}_s=-{\cd}_0 \left[ V_0 + \|\vec{V_1}\| \right] \dep{q_1-q_s}
\end{equation}
where $q_1$ (resp. $q_s$) is the value of $q$ in the first atmospheric
layer (resp. at the surface),
%
$V_1=\sqrt{u_1^2+v_1^2}$ and $V_0$ is a
minimum wind value fixed to 1~m~s$^{-1}$.


\subsubsection{Convective adjustment}

The occurrence of vertical unstable temperature profiles
\begin{equation}
\der{\t}{z}<0
\end{equation}
is prevented using a simple dry convective adjustment.
If such an unstable profile is produced by the model,
an adiabatic profile is immediately restored
with a simple energy conserving scheme. If the resulting temperature
profile is unstable at its upper or lower limit this mechanism is
instantaneously extended in such a way that the final profile
is entirely stable.

This convective adjustment is in fact achieved in a real atmosphere
by parcel exchange through vertical convective motions.
These motions not only transport potential temperature but also
momentum.
The momentum exchange can only be estimated from a diagnostic of the
vertical mass flux which can not be obtained from such an unphysical
parametrization.
Nevertheless, an estimate of the degree of instability of the vertical
profile can be obtained from the
relative enthalpy exchange necessary to restore the diabatic profile
$\bar{\t}$ from the original profile $\t$
\begin{equation}
\alpha=\int{\left| \t-\bar{\t} \right| \rho dz}
/\int{\t \rho dz}
\end{equation}
When $\alpha<1$ (this condition is always satisfied in the
simulations) the momentum is entirely mixed on a fraction
$\alpha$ of the mesh. This arbitrary choice is qualitatively
acceptable in the sense that a greater instability will produce
a greater momentum mixing.

\subsubsection{Surface processes}

Temperature at the surface of a planet is governed by the balance
between incoming fluxes (direct solar insolation,
thermal radiation from the atmosphere and the surface itself and
turbulent heat fluxes) 
and  thermal conduction in the soil.
The time evolution of the temperature under the surface is given
by a classical conduction equation
%
\begin{equation}
\der{T}{t}=-\frac{1}{C}\der{F_c}{z}
\end{equation}
%
where the conductive flux $F_c$ is given by
%
\begin{equation}
F_c=-\lambda \der{T}{z}
\end{equation}
%
and where $\lambda$ and $C$ are the soil conductivity and volumetric
heat, respectively.
In the simple case of a vertically homogeneous soil
(which is assumed in the model)
it can be shown that, as far as the evaluation of surface
temperatures and fluxes is concerned, the model depends on only one free
parameter: the so-called thermal inertia $I=\sqrt{\lambda C}$.

The conduction in the soil is computed using a finite difference scheme
analogous to that developed by Jacobsen and Heise,
1982~\nocite{Jaco:82}.
Some complementary information on the numerical scheme is given
by Hourdin \etal, 1993. \nocite{Hour:93}

\def\unefig#1#2{ #1) & \includegraphics[width=13cm]{FIGS/#2.pdf}}
\begin{figure}
\centerline{
\begin{tabular}{cc}
\unefig{a}{tterre} \\
\unefig{b}{lcterre} \\
\unefig{c}{vterre} \\
\unefig{d}{uterre}
\end{tabular}
}
\caption{ \label{fg:terre}
%titre de la legende:
Results for a pseudo dry Earth for December-January-February
}

{\bf (a)} Temperature in $^oC$ (upper panel), {\bf (b)} stream function of the mean meridional
circulation (10$^9$~kg~s$^{-1}$), {\bf (c)} meridional wind in m~s$^{-1}$
(central panel) and {\bf (d) } zonal wind in m~s$^{-1}$ (lower panel) have been
average in time over three months and in longitudes.
\end{figure}

% \bibliography{../biblio/fred,/u1/gcm/biblio/lmd}

\bibliography{/home/hourdin/TEX/BIBLIO/fred}


\end{document}
