source: trunk/LMDZ.MARS/doc/principe.tex @ 3026

Last change on this file since 3026 was 1954, checked in by emillour, 7 years ago

Mars GCM:

  • Make a "doc" subdirectory to store the documentation source files with the code.

EM

File size: 13.6 KB
RevLine 
[1954]1\chapter{Main features of the model}
2
3\label{sc:apercu}
4
5\section{Basic principles}
6The General Circulation Model (GCM) calculates the temporal evolution of
7the different {\bf variables} (listed below)
8that control or describe the Martian meteorology and climate
9at different points of a {\bf 3D ``grid" } (see below) that covers
10the entire  atmosphere.
11
12From an initial state, the model calculates the evolution of these variables,
13timestep by timestep:
14\begin{itemize}
15\item At instant $t$, we know variable $X_t$ (temperature for example)
16at one point in the atmosphere.
17
18\item We calculate the evolution (the {\bf tendencies})
19$(\frac{\partial X}{\partial t})_1$ ,
20$(\frac{\partial X}{\partial t})_2$ , etc.
21arising from each physical phenomenon,
22calculated by a {\bf parameterization} of each of these phenomenon
23(for example, heating due to absorption of solar radiation).
24
25\item At the next time step $t + \delta t$, we can calculate $X_{t+ \delta t}$
26from  $X_t$ and $(\frac{\partial X}{\partial t})$.
27This is the {\bf``integration''} of the variables in time.
28(For example, $X_{t+ \delta t}=X_t +
29 \delta t(\frac{\partial X}{\partial t})_1 +
30 \delta t(\frac{\partial X}{\partial t})_2 + ...$)
31
32\end{itemize} 
33
34{\bf The main task of the model is to calculate these tendencies}
35$(\frac{\partial X}{\partial t})$
36{\bf arising from the different parameterized phenomenon.}
37
38\section{Dynamical-Physical separation}
39
40In practice, the 3D model operates in two parts:\\
41- one {\bf dynamical part} containing the numerical solution of
42the general equations for atmospheric circulation.
43This part (including the programming) is common to the Earth and Martian model,
44and in general for all atmospheres of the terrestrial type.\\
45- a second {\bf physical part} that is specific to the planet in question and
46which calculates the forced circulation and the climate details at each point.
47
48The calculations for the dynamical part are made on a 3D grid
49with horizontal exchanges between the grid boxes,
50whereas the physical part can be seen as a juxtaposition of atmosphere
51``columns'' that do not interact with each other (see diagram~\ref{fg:fidyn}).
52
53The dynamical and physical parts deal with variables of different natures,
54and operate on grids that are differently constructed.
55The temporal integration of the variables is based
56on different numerical schemes
57(simple, such as the one above for the physical part,
58and more complicated, the ``Matsuno-Leapfrog'' scheme for the dynamical part).
59The timesteps are also different.
60The physical timestep  is {\tt iphysiq} times longer than the dynamical
61timestep, as the solution of the dynamic equations requires a shorter timestep
62than the forced calculation for the physical part.
63
64In practice, the main program that handles the whole model (\verb+gcm.F+)
65is located in the dynamical part.
66When the temporal evolution is being calculated,
67at each timestep the program calls the following:
68\begin{enumerate}
69\item Call to the subroutine that handles the total tendency calculation
70$(\frac{\partial X}{\partial t})$ arising from the dynamical part
71(\verb+caldyn.F+)
72\item Integration of these dynamical tendencies to calculate the evolution
73of the variables at the following timesteps (subroutine \verb+integrd.F+)
74\item Every {\tt iphysiq} dynamical timestep, a call to the interface
75subroutine (\verb+calfis.F+) with the physical model (\verb+physiq.F+),
76that calculates the evolution of some of the purely physical variables
77(e.g: surface temperature {\tt tsurf}) and returns the tendencies 
78$(\frac{\partial X}{\partial t})$
79arising from the physical part.
80\item Integration of the physical variables (subroutine \verb+addfi.F+)
81\item Similarly, calculation and integration of tendencies due to
82the horizontal dissipation and the ``sponge layer'' is done
83every {\tt idissip} dynamical time step.
84\end{enumerate} 
85{\em {Remark:} The physical part can be run separately for a 1-D calculation
86for a single column using program \verb+testphys1d.F+.} 
87 
88
89\section{Grid boxes:} Examples of typical grid values are 64x48x25,
9064x48x32 or 32x24x25 in longitudexlatitudexaltitude.
91For Mars (radius$\sim$3400~km), a 64x48 horizontal grid corresponds
92to grid boxes of the order of 330x220 kilometers near the equator.
93
94
95\subsection{Horizontal grids}
96
97\begin{figure}
98\begin{center}
99\centerline{\framebox{\epsfig{figure=Fig/didi.pdf,width=12.cm,clip=true}}}
100\caption{Physical/dynamical interface}
101\label{fg:fidyn}
102\end{center}
103\end{figure}
104
105\begin{figure}
106\centerline{\framebox[1.4\textwidth][c]{\includegraphics[width=1.2\textwidth]
107{Fig/grid.pdf}}}
108\caption{Dynamical and physical grids for a 6 $\times$ 7 horizontal resolution.
109In the dynamics (but not in the physics) winds u and v are on specific
110staggered grids. Other dynamical variables are on the dynamical ``scalar'' grid.
111The physics uses the same ``scalar'' grid for all the variables,
112except that nodes are indexed in a single vector containing
113NGRID=2+(JM-1)$\times$IM points when counting from the north pole.\\
114N.B.: In the Fortran program, the following variables are used:
115 {\tt iim=IM , iip1=IM+1, jjm=JM , jjp1=JM+1}.}
116\label{fg:grid}
117\end{figure}
118
119
120
121
122
123
124Dynamics and physics use different grids.
125Figure~\ref{fg:grid} shows the correspondance and indexing
126of the physical and dynamical grids
127as well as the different locations of variables on these grids.
128To identify the coordinates of a variable (at one grid point up,
129down, right or left)
130we use coordinates {\tt rlonu}, {\tt rlatu},
131 {\tt rlonv}, {\tt rlatv} (longitudes and
132latitudes, in {\bf radians}).
133
134On the dynamical grid, values at {\tt i=1} are the same as at {\tt i=IM+1}
135as the latter node is a redundant point (due to the periodicity in longitude,
136these two nodes are actualy located at the same place).
137Similarly, the extreme {\tt j=1} and {\tt j=JM+1} nodes on the
138dynamical grid (respectively corresponding to North and South poles)
139are duplicated IM+1 times.\\
140In contrast, the physical grid does not contain redundant points
141(only one value for each pole and no extra point along longitudes),
142as shown in figure~\ref{fg:grid}.
143In practice, computations relative to the physics are made
144for a series of {\tt ngrid} atmospheric
145columns, where {\tt NGRID}={\tt IM}x({\tt JM}-1)+2.
146
147\subsection{Vertical grids}
148
149\begin{figure}[htbp]
150                                                                               
151{\includegraphics[scale=0.85,clip=true]{Fig/coord.pdf}}
152\caption[hybrides]
153 {Sketch illustrating the difference between hybrid and non-hybrid coordinates}
154 \label{fg:hybrid}
155\end{figure}
156
157
158The GCM was initially programmed using sigma coordinates $\sigma = p/ps$
159(atmospheric pressure over surface pressure ratio)
160which had the advantage of using a constant domain
161($\sigma=1$ at the surface and $\sigma=0$ at the top of the atmosphere)
162whatever the underlying topography.
163However, it is obvious that these coordinates significantly disturb
164the stratospheric dynamical representation as the topography is propagated
165to the top of the model by the coordinate system.
166This problem can elegantly be solved by using a hybrid sigma-P (sima-pressure)
167hybrid coordinate which is equivalent to using $\sigma$ coordinates
168near the surface and gradualy shifting to purely pressure $p$ 
169coordinates with increasing altitude.
170Figure~\ref{fg:hybrid} illustrates the importance of using these
171hybrid coordinates compared to simple $\sigma$ coordinates.
172The distribution of the vertical layers is irregular,
173to enable greater precision at ground level.
174In general we use 25 levels to describe the atmosphere to a height of 80 km,
17532 levels for simulations up to 120~km, or 50 levels to rise
176up to thermosphere.
177The first layer describes the first few meters above the ground,
178whereas the upper layers span several kilometers.
179Figure~\ref{fg:sigma} describes the vertical grid representation
180and associated variables.
181
182\begin{figure}
183\begin{verbatim}
184DYNAMICS                                               PHYSICS
185--------                                               -------
186[coordinates ap(),bp()]                                [pressures]
187
188ap(llm+1)=0,bp(llm+1)=0  ****************************  plev(nlayer+1)=0
189aps(llm),bps(llm)        .. llm .......... nlayer ...  play(nlayer)
190ap(llm),bp(llm)          ****************************  plev(nlayer)
191aps(llm-1),bps(llm-1)    .. llm-1 ........ nlayer-1 .  play(nlayer-1)
192ap(llm-1),bp(llm-1)      ****************************  plev(nlayer-1)
193                             :                :
194                             :                :
195                             :                :
196aps(2),bps(2)            ... 2  ............. 2  ....  play(2)
197ap(2),bp(2)              ****************************  plev(2)
198aps(1),bps(1)            ... 1  ............. 1  ....  play(1)
199ap(1)=1,bp(1)=0          **********surface***********  plev(1)=Ps
200\end{verbatim}
201\caption{Vertical grid description of the {\tt llm (or nlayer)}
202atmospheric layers
203in the programming code ({\tt llm} is the variable used in the dynamical part,
204and {\tt nlayer} is used in the physical part).
205Variables {\tt ap, bp} and {\tt aps, bps} indicate the hybrid levels at
206the interlayer levels and at middle of the layers respectively.
207Pressure at the interlayer is  $Plev(l)=ap(l)+bp(l) \times Ps$ and pressure
208in the middle of the layer is defined by $Play(l)=aps(l)+bps(l) \times Ps$,
209(where $Ps$ is surface pressure).
210Sigma coordinates are merely a specific case of hybrid
211coordinates such that $aps=0$ and $bps=P/Ps$.
212Note that for the hybrid coordinates, $bps=0$ above $\sim50$~km, leading to
213purely pressure levels.
214The user can choose whether to run the model using hybrid coordinates or not
215by setting variable {\tt hybrid} in run.def to True or False.}
216\label{fg:sigma}
217\end{figure}
218
219\section{Variables used in the model}
220
221\subsection{Dynamical variables}
222The dynamical state variables are the atmospheric temperature,
223surface pressure, winds and tracer concentrations.
224In practice, the formulation selected to solve the equations in
225the dynamics %(see chapter~\ref{sc:dynamic})
226is optimised using the following less ``natural'' variables:
227
228\begin{description}
229\item - {\bf potential temperature} $\theta$ ({\tt teta} in the code),
230linked to temperature {\bf T} by $\theta = T{(P/Pref)}^{-\kappa}$
231with $\kappa = R/C_p$
232(note that $\kappa$ is called {\tt kappa} in the dynamical code, and {\tt rcp}
233in the physical code). We take $Pref=610$ Pa on  Mars.
234\item - {\bf surface pressure} ({\tt ps} in the code).
235\item - {\bf mass} the atmosphere mass in each grid box ({\tt masse}
236in the code).
237\item - {\bf the covariant meridional and zonal winds}
238{\tt ucov} and {\tt vcov}.
239These variables are linked to the "natural" winds by
240\verb+ucov = cu * u+ and \verb+vcov = cv * v+, where
241\verb+cu+ and \verb+cv+ are constants that only depend on the latitude.
242%(see Appendix A).
243\item - {\bf mixing ratio of tracers} in the atmosphere, typically
244expressed in kg/kg (array {\tt q} in the code).
245\end{description} 
246
247{\tt ucov} and {\tt vcov}, ``vectorial'' variables, are stored on
248``scalari'' grids u and v respectively, in the dynamics
249(see section \ref{fg:grid}).
250{\tt teta}, {\tt q}, {\tt ps}, {\tt masse}, ``scalar variables'',
251are stored on the ``scalar'' grid of the dynamics.
252
253
254
255\subsection{Physical variables}
256
257In the physics, the state variables of the dynamics are transmitted
258via an interface that interpolates the winds on the scalar grid
259(that corresponds to the physical grid) and transforms the dynamical variables
260into more ``natural'' variables.
261Thus we have winds {\bf u} and  {\bf v} (m.s$^{-1}$),
262temperature {\bf T} (K), pressure at the middle of the layers
263{\bf play} (Pa) and at interlayers {\bf plev} (Pa), tracers
264{\bf q}, etc. (kg/kg) on the same grid.
265
266Furthermore, the physics also handle the evolution of the purely
267physical state variables:
268\begin{description}
269\item - {\bf co2ice} CO$_2$ ice on the surface (kg.m$^{-2}$)
270\item - {\bf tsurf} surface temperature (K),
271\item - {\bf tsoil} temperature at different layers under the surface (K),
272\item - {\bf emis} surface emissivity,
273\item - {\bf q2} wind variance, or more precisely the square root of
274the turbulent kinetic energy.
275\item - {\bf qsurf} tracer on the surface (kg.m$^{-2}$).
276\end{description}
277
278\subsection{Tracers}
279The model may include different types of tracers:
280\begin{description}
281\item - dust particles, which may have several modes
282\item - chemical species which depict the chemical composition of the
283atmosphere %\\ \\
284%$co2      =  nqchem\_min\\
285%co       =  nqchem\_min + 1\\
286%o        =  nqchem\_min + 2\\
287%o(1d)    =  nqchem\_min + 3\\
288%o2       =  nqchem\_min + 4\\
289%o3       =  nqchem\_min + 5\\
290%h        =  nqchem\_min + 6\\
291%h2       =  nqchem\_min + 7\\
292%oh       =  nqchem\_min + 8\\
293%ho2      =  nqchem\_min + 9\\
294%h2o2     =  nqchem\_min + 10\\
295%n2       =  nqchem\_min + 11\\
296%ar       =  nqchem\_min + 12\\ \\ $
297%If you choose not to handle dust particles (general case), then nqchem\_min=1
298\item - water vapor and water ice particles %\\
299%h2o corresponds to the nqmx tracer (nqmx is chosen when compiling
300%with the option {\tt -t nqmx})\\
301%If the option "ice" is chosen, water ice corresponds to the (nqmx-1)th tracer.
302\item - ...
303\end{description} 
304
305
306In the code, all tracers are stored in one three-dimensional array {\tt q},
307the third index of which corresponds to each individual tracer.
308In input and output files (``start.nc'', ``startfi.nc'',
309see Section~\ref{loc:contact1}) tracers are stored seperately using their
310individual names. Loading specific tracers requires that
311the approriate tracer names are set in the {\tt traceur.def} file
312(see Section~\ref{sc:traceur.def}), and specific computations
313for given tracers (e.g.: computing the water cycle, chemistry in the
314upper atmosphere, ...) is controled by setting coresponding options in the
315{\tt callphys.def} file (see Section~\ref{sc:callphys.def}).
Note: See TracBrowser for help on using the repository browser.