Dust Cycle in Mars PCM5

2025-08-22T14:27:41Z

Emillour: Created page with "This page details at different levels how to use and what is featured in the dust cycle in Mars PCM5. == Brief Overview == The reference work leading to the implementation in..."

PCM vertical coordinate

2025-08-21T15:07:06Z

Emillour:

== Overview ==
The PCM vertical coordinate, also called "model levels" is rather specific in the PCM and typically consists of "sigma" or hybrid "sigma-pressure" coordinates. What one needs to keep in mind is that model levels are not at fixed altitude, nor in most cases at fixed pressure. The pressure $$P$$ of a model level $$k$$ at time $$t$$ is:
$$
\begin{align}
P(k,t) & = ap(k) + bp(k) . Ps(t)
\end{align}
$$
Where $$ap(k)$$ and $$bp(k)$$ are respectively hybrid pressure and hybrid sigma coefficients. Note that $$ap()$$ and $$bp()$$ are time-independent and $$Ps(t)$$ is surface pressure (which typically varies with time).

=== sigma coordinates ===
If in the general equation above one sets $$ap(k)=0$$ for all $$k$$, then $$P(k,t)=bp(k).Ps(t)$$, which can be re-written as $$bp(k)=P(k,t)/Ps(t)$$, which is the expression of a "sigma" coordinate (often represented with the notation $$\sigma=P/Ps$$). In case of a hydrostatic system, pressure monotonically decreases with altitude and thus "sigma" is indeed a coordinate, which monotonically decreases from 1 at the surface (where $$P=P_s$$) to 0 at the top of the atmosphere (where $$P=0$$). "Sigma" coordinates are also often called "terrain-following" coordinates because for given model level $$k$$ the pressure is a constant factor of the surface pressure, which to first order depends on the topography.

=== pressure coordinates ===
If in the general equation above one sets $$bp(k)=0$$ for all $$k$$, then $$P(k,t)=ap(k)$$, which implies that a layer k is always at the same pressure. While this is quite fine at high enough altitudes, it is quite easy to see that in the near surface and presence of topography (or if fact even without topography when weather systems will impact on local surface pressure) it will be problematic to use a fixed pressure grid where some points will be undefined.

=== hybrid sigma-pressure coordinates ===
The two approaches above can be combined by an appropriate choice of the values of $$ap(k)$$ and $$bp(k)$$, i.e. enforcing that:
* near the surface (small values of $$k$$) $$ap(k) \simeq 0 $$ (i.e. small compared to $$ bp(k) . Ps(t)$$), then the vertical coordinate there is close to being a purely "sigma" coordinate
* high enough (i.e. high enough above the topography, high values of $$k$$) $$bp(k) \simeq 0$$ and $$ bp(k). Ps(t) << ap(k)$$, then the vertical coordinate there is close to being a purely "pressure" coordinate

== In practice ==
Information about the vertical coordinate is only partially stored in the PCM (lon-lat) dynamics start files which contain the number of atmospheric layer (llm), the ap() and bp() coefficients (along with the the $$preff$$ and $$pa$$ parameters discussed below, which are in the "controle" array). Note that the DYNAMICO dynamics start file do not contain this information.

In practice the vertical layer discretisation is generated at run time, using information from the input [[The z2sig.def Input File|z2sig.def file]] and a couple of related parameters : $$preff$$, reference surface pressure, (in Pa), and $$pa$$, a nameless parameter homogeneous to a pressure (also expressed in Pa) which roughly corresponds to a "reference transition pressure" from which hybrid-sigma coordinates become essentially pressure coordinates. Based on these parameters, and the specified atmospheric scale height and related target pseudo-altitudes of model levels from file [[The z2sig.def Input File|z2sig.def]] the model iteratively (requires solving a non-linear problem with no analytic solution) computes and builds the sought model levels.

An additional input parameter that will impact on the generation of model levels is the $$hybrid$$ (logical) parameter set in [[The run.def Input File|run.def]], e.g.:
<pre>
# use hybrid vertical coordinate (else will use sigma levels)
hybrid=.true.
</pre>
If $$hybrid$$ is set to $$.false.$$ the generated model levels will be purely sigma levels, i.e. the $$ap()$$ coefficients will be zero.

You can check out the code where this is done, in the ''disvert_noterre'' routine: https://trac.lmd.jussieu.fr/Planeto/browser/trunk/LMDZ.COMMON/libf/dyn3d_common/disvert_noterre.F and/or its "mirror" version in the dynamico-LMDZ interface ''disvert_icosa_lmdz'' https://trac.lmd.jussieu.fr/Planeto/browser/trunk/ICOSA_LMDZ/src/disvert_icosa_lmdz.f90

[[Category:Generic-Model]]
[[Category:Mars-Model]]
[[Category:Venus-Model]]
[[Category:Titan-Model]]
[[Category:Pluto-Model]]

PCM vertical coordinate

2025-08-21T14:58:34Z

Emillour:

== Overview ==
The PCM vertical coordinate, also called "model levels" is rather specific in the PCM and typically consists of "sigma" or hybrid "sigma-pressure" coordinates. What one needs to keep in mind is that model levels are not at fixed altitude, nor in most cases at fixed pressure. The pressure $$P$$ of a model level $$k$$ at time $$t$$ is:
$$
\begin{align}
P(k,t) & = ap(k) + bp(k) . Ps(t)
\end{align}
$$
Where $$ap(k)$$ and $$bp(k)$$ are respectively hybrid pressure and hybrid sigma coefficients (note that $$ap()$$ and $$bp()$$ are time-independent) and $$Ps(t)$$ is surface pressure (and typically varies with time).

=== sigma coordinates ===
If in the general equation above one sets $$ap(k)=0$$ for all $$k$$, then $$P(k,t)=bp(k).Ps(t)$$, which can be re-written as $$bp(k)=P(k,t)/Ps(t)$$, which is the expression of a "sigma" coordinate (often represented with the notation $$\sigma=P/Ps$$). In case of a hydrostatic system, pressure monotonically decreases with altitude and thus "sigma" is indeed a coordinate, which monotonically decreases from 1 at the surface (where $$P=P_s$$) to 0 at the top of the atmosphere (where $$P=0$$). "Sigma" coordinates are also often called "terrain-following" coordinates because for given model level $$k$$ the pressure is a constant factor of the surface pressure, which to first order depends on the topography.

=== pressure coordinates ===
If in the general equation above one sets $$bp(k)=0$$ for all $$k$$, then $$P(k,t)=ap(k)$$, which implies that a layer k is always at the same pressure. While this is quite fine at high enough altitudes, it is quite easy to see that in the near surface and presence of topography (or if fact even without topography when weather systems will impact on local surface pressure) it will be problematic to use a fixed pressure grid where some points will be undefined.

=== hybrid sigma-pressure coordinates ===
The two approaches above can be combined by an appropriate choice of the values of $$ap(k)$$ and $$bp(k)$$, i.e. enforcing that:
* near the surface (small values of $$k$$) $$ap(k) \simeq 0 $$ (i.e. small compared to $$ bp(k) . Ps(t)$$), then the vertical coordinate there is close to being a purely "sigma" coordinate
* high enough (i.e. high enough above the topography, high values of $$k$$) $$bp(k) \simeq 0$$ and $$ bp(k). Ps(t) << ap(k)$$, then the vertical coordinate there is close to being a purely "pressure" coordinate

== In practice ==
Information about the vertical coordinate is only partially stored in the PCM (lon-lat) dynamics start files which contain the number of atmospheric layer (llm), the ap() and bp() coefficients (along with the the $$preff$$ and $$pa$$ parameters discussed below, which are in the "controle" array). Note that the DYNAMICO dynamics start file do not contain this information.

In practice the vertical layer discretisation is generated at run time, using information from the input [[The z2sig.def Input File|z2sig.def file]] and a couple of related parameters : $$preff$$, reference surface pressure, (in Pa), and $$pa$$, a nameless parameter homogeneous to a pressure (also expressed in Pa) which roughly corresponds to a "reference transition pressure" from which hybrid-sigma coordinates become essentially pressure coordinates. Based on these parameters, and the specified atmospheric scale height and related target pseudo-altitudes of model levels from file [[The z2sig.def Input File|z2sig.def]] the model iteratively (requires solving a non-linear problem with no analytic solution) computes and builds the sought model levels.

An additional input parameter that will impact on the generation of model levels is the $$hybrid$$ (logical) parameter set in [[The run.def Input File|run.def]], e.g.:
<pre>
# use hybrid vertical coordinate (else will use sigma levels)
hybrid=.true.
</pre>
If $$hybrid$$ is set to $$.false.$$ the generated model levels will be purely sigma levels, i.e. the $$ap()$$ coefficients will be zero.

You can check out the code where this is done, in the ''disvert_noterre'' routine: https://trac.lmd.jussieu.fr/Planeto/browser/trunk/LMDZ.COMMON/libf/dyn3d_common/disvert_noterre.F and/or its "mirror" version in the dynamico-LMDZ interface ''disvert_icosa_lmdz'' https://trac.lmd.jussieu.fr/Planeto/browser/trunk/ICOSA_LMDZ/src/disvert_icosa_lmdz.f90

[[Category:Generic-Model]]
[[Category:Mars-Model]]
[[Category:Venus-Model]]
[[Category:Titan-Model]]
[[Category:Pluto-Model]]

Dissipation

2025-08-21T14:47:57Z

Emillour:

Note that there are some slides from the LMDZ Earth Climate Model tutorial sessions:
https://lmdz.lmd.jussieu.fr/utilisateurs/formation/winter_2024
in the "Dynamics: grid/temporal discretization/stability/dissipation" course which can be helpful to learn about dissipation in the PCM,
why it is there, how it is handled, and what are the underlying parameters.

== Description==
In the LMDZ grid point model, nonlinear interactions between explicitly resolved scales
and subgrid-scale processes are parameterized by applying a scale-selective horizontal dissipation operator based on an $$n$$ time iterated Laplacian $$\Delta^n$$. For the grid point model, for
instance, this can be written:

\begin{align}
\label{def:Wns}
\frac{\partial q}{\partial t} = \frac{(-1)^n}{\tau_{diss}}(\delta x)^{2n}\Delta^nq
\end{align}

where $$q$$ is a field component on which dissipation is applied, $$\delta x$$ is the smallest horizontal distance represented in the model and $$\tau_{diss}$$ is the dissipation timescale for a structure of scale $$\delta x$$. These operators are necessary to ensure the grid point model numerical stability. In practice, the operator is separately applied to three components :
* the divergence of the flow,
* the vorticity of the flow,
* potential temperature.

We classically use n = 1 for the divergence of the flow (<code>nitergdiv=1</code> in [[The_run.def_Input_File|run.def]]) and n = 2 for flow vorticity and potential temperature (<code>nitergrot=2</code> and <code>niterh=2</code> in [[The_run.def_Input_File|run.def]]).

== How to change it in the model ==

In practice, the values of $$n$$ and $$\tau_{diss}$$ are prescribed in the [[The_run.def_Input_File|run.def]] with the keys:
*nitergdiv
*nitergrot
*niterh
for the values of $$n$$ on each field, and the associated time scales $$\tau$$ (in s):
*tetagdiv
*tetagrot
*tetatemp

In [[The_run.def_Input_File|run.def]], there is also a key ''idissip'' (depreciated) or ''dissip_period'' which is the rate (in dynamical steps) at which to apply the dissipation (and thus impacts on the overall time step of the dissipation).

In addition there are multiplicative factors on the dissipation that are applied going to the upper layers of the atmosphere (usually more dissipation is required in the upper atmosphere than the lower atmosphere since perturbations traveling upwards will grow in amplitude). These coefficients are:
* dissip_fac_mid : dissipation multiplicative factor in the middle atmosphere
* dissip_fac_up : dissipation multiplicative factor in the upper atmosphere

At which pressure/altitude the middle and upper atmosphere referred to above depends on another input flag: vert_prof_dissip
* if <code>vert_prof_dissip=0</code> (default setup for most planets) then the bottom of the transition zone between mid and upper atmosphere is set at the pressure (in Pa) provided by the input parameter dissip_pupstart (set in [[The_run.def_Input_File|run.def]]) and extends over dissip_deltaz km (also set in [[The_run.def_Input_File|run.def]]), assuming a transition zone scale height of dissip_hdelta km.
* if <code>vert_prof_dissip=1</code> (default setup if <code>planet_type=="mars"</code>!) then the transition between mid and upper atmosphere is set to start at a pseudo-altitude of 70 km, with a transition region of 30km (see inidissip.F for the details)

When the PCM runs it yields a digest of the dissipation multiplicative coefficients and timesteps (look for keyword "inidissip" in the output!) of the likes of:
<pre>
Dissipation :
Multiplication de la dissipation en altitude :
dissip_fac_mid = 2.0000000000000000
dissip_fac_up = 20.000000000000000
Transition mid /up: Pupstart,delta = 1000.0000000000000 Pa 10.000000000000000 (km)
inidissip: Time constants for lateral dissipation
inidissip: dissip_period= 5 dtdiss= 870.83333333333326 dtvr= 174.16666666666666

pseudoZ(km) zvert dt(tetagdiv) dt(tetagrot) dt(divgrad)
0.0 1.0000004E+00 8.7083367E-02 8.7083367E-02 8.7083367E-02
0.0 1.0000096E+00 8.7084169E-02 8.7084169E-02 8.7084169E-02
0.1 1.0000957E+00 8.7091668E-02 8.7091668E-02 8.7091668E-02
0.2 1.0005771E+00 8.7133585E-02 8.7133585E-02 8.7133585E-02
0.5 1.0024163E+00 8.7293751E-02 8.7293751E-02 8.7293751E-02
0.9 1.0079124E+00 8.7772369E-02 8.7772369E-02 8.7772369E-02
1.4 1.0217973E+00 8.8981516E-02 8.8981516E-02 8.8981516E-02
2.1 1.0526661E+00 9.1669671E-02 9.1669671E-02 9.1669671E-02
3.1 1.1136968E+00 9.6984432E-02 9.6984432E-02 9.6984432E-02
4.3 1.2193074E+00 1.0618135E-01 1.0618135E-01 1.0618135E-01
5.7 1.3732887E+00 1.1959055E-01 1.1959055E-01 1.1959055E-01
7.5 1.5542353E+00 1.3534799E-01 1.3534799E-01 1.3534799E-01
9.6 1.7216407E+00 1.4992621E-01 1.4992621E-01 1.4992621E-01
12.1 1.8454911E+00 1.6071151E-01 1.6071151E-01 1.6071151E-01
14.9 1.9221312E+00 1.6738559E-01 1.6738559E-01 1.6738559E-01
18.1 1.9631130E+00 1.7095442E-01 1.7095442E-01 1.7095442E-01
21.4 1.9826551E+00 1.7265621E-01 1.7265621E-01 1.7265621E-01
24.8 1.9916537E+00 1.7343984E-01 1.7343984E-01 1.7343984E-01
28.1 1.9959119E+00 1.7381066E-01 1.7381066E-01 1.7381066E-01
31.4 1.9979711E+00 1.7398998E-01 1.7398998E-01 1.7398998E-01
34.8 1.9989834E+00 1.7407814E-01 1.7407814E-01 1.7407814E-01
38.1 1.9994871E+00 1.7412200E-01 1.7412200E-01 1.7412200E-01
41.4 1.9997400E+00 1.7414402E-01 1.7414402E-01 1.7414402E-01
44.8 1.9998677E+00 1.7415514E-01 1.7415514E-01 1.7415514E-01
48.1 1.9999325E+00 1.7416079E-01 1.7416079E-01 1.7416079E-01
51.4 1.9999655E+00 1.7416366E-01 1.7416366E-01 1.7416366E-01
54.8 1.9999823E+00 1.7416513E-01 1.7416513E-01 1.7416513E-01
58.1 1.9999910E+00 1.7416588E-01 1.7416588E-01 1.7416588E-01
61.4 1.9999954E+00 1.7416626E-01 1.7416626E-01 1.7416626E-01
64.8 1.9999976E+00 1.7416646E-01 1.7416646E-01 1.7416646E-01
68.1 1.9999988E+00 1.7416656E-01 1.7416656E-01 1.7416656E-01
71.4 1.9999997E+00 1.7416664E-01 1.7416664E-01 1.7416664E-01
74.8 2.0000019E+00 1.7416683E-01 1.7416683E-01 1.7416683E-01
78.1 2.0000163E+00 1.7416809E-01 1.7416809E-01 1.7416809E-01
81.4 2.0001218E+00 1.7417728E-01 1.7417728E-01 1.7417728E-01
84.8 2.0009008E+00 1.7424511E-01 1.7424511E-01 1.7424511E-01
88.1 2.0066545E+00 1.7474616E-01 1.7474616E-01 1.7474616E-01
91.4 2.0490550E+00 1.7843854E-01 1.7843854E-01 1.7843854E-01
94.8 2.3562671E+00 2.0519160E-01 2.0519160E-01 2.0519160E-01
98.1 4.3369577E+00 3.7767674E-01 3.7767674E-01 3.7767674E-01
101.4 1.1438614E+01 9.9611261E-01 9.9611261E-01 9.9611261E-01
104.8 1.8031963E+01 1.5702835E+00 1.5702835E+00 1.5702835E+00
108.1 1.9705847E+01 1.7160508E+00 1.7160508E+00 1.7160508E+00
111.4 1.9959620E+01 1.7381503E+00 1.7381503E+00 1.7381503E+00
114.8 1.9994525E+01 1.7411898E+00 1.7411898E+00 1.7411898E+00
118.5 1.9999423E+01 1.7416164E+00 1.7416164E+00 1.7416164E+00
123.5 1.9999971E+01 1.7416641E+00 1.7416641E+00 1.7416641E+00
129.8 1.9999999E+01 1.7416666E+00 1.7416666E+00 1.7416666E+00
137.1 2.0000000E+01 1.7416667E+00 1.7416667E+00 1.7416667E+00
144.5 2.0000000E+01 1.7416667E+00 1.7416667E+00 1.7416667E+00
</pre>

== Good to know rules and rules of thumb ==
* If your simulation shows numerical instabilities, a good idea is to increase dissipation. This means decreasing parameters $$\tau$$ and/or increasing dissip_fac_up (and to a lesser extent dissip_fac_mid).
* Optimal values for the dissipation timescales depends on the resolution of the horizontal grid. The higher the resolution, the more dissipation we need.
* Because dissipation is applied using an Explicit Euler time marching scheme, it is liable to be unstable if the time step is too large. This is tested at run-time by the model which will generate an error message of the likes of:
<pre>
STOP : lateral dissipation is too intense and will
generate instabilities in the model !
You must increase tetah (or increase dissip_period
or increase day_step)
</pre>

[[Category:Generic-Model]]
[[Category:Mars-Model]]
[[Category:Venus-Model]]

Dissipation

2025-08-21T13:57:05Z

Emillour:

2025-04-30T09:00:56Z

Emillour:

The ''diagfi.def'' file is an optional file that is read by the GCM at run time. If it is present then it will dictate which fields will be included in the ''diagfi.nc'' output file.

A clarification: generating the ''diagfi.nc'' output file is only possible when running with the lon-lat (LMDZ.COMMON) dynamical core (or in 1D), with the Generic or Mars physics packages.

== ''diagfi.def '' format ==
Each line of this ASCII file should contain the name (case sensitive!) of a variable to output.

Note that the variable name to use is the one set in the code when calling the ''writediagfi'' or ''write_output'' routines. e.g.
<syntaxhighlight lang="fortran">
call writediagfi(ngrid,"temperature","temperature","K",3,zt)
</syntaxhighlight>
or
<syntaxhighlight lang="fortran">
call write_output("temperature","temperature","K",zt)
</syntaxhighlight>

implies that the variable ''temperature'' can be outputted in ''[[diagfi.nc]]''

== Example of ''diagfi.def'' file ==

<pre>
aire
altitude
ap
aps
ASR
ASRcs
beta
bp
bps
CLF
CLFt
controle
Declin
dt_ekman1
dt_ekman2
dt_diff1
dt_diff2
DYN
evap_surf_flux
fluxsurf_rad
GND
h2o_ice
h2o_ice_col
h2o_ice_surf
h2o_vap
h2o_vap_col
h2o_vap_surf
H2Oice_reffcol
ISR
latentFlux
latitude
longitude
Ls
Lss
mass_evap_col
OLR
OLRcs
p
pctsrf_sic
phisinit
ps
RA
rain
reevap
RH
rnat
sea_ice
sensibFlux
shad
snow
soildepth
tau_col
temp
Time
Tsat
tsea_ice
tslab1
tslab2
tsurf
u
v
w
</pre>

== Various related comments and remarks ==
* If there is no ''diagfi.def'' file found at run-time, then '''ALL''' variables will be written in ''diagfi.nc'', which will rapidly become huge. You have been warned. Note however that it is a convenient way to have access to the names of all available variables that may be outputted in ''diagfi.nc''.
* The rate at which the outputs are made in ''diagfi.nc'' is set by the ''diagfi_output_rate'' parameter which is usually set in the ''callphys.def'' file and expressed in terms of physical time steps for the Generic PCM. For the Mars PCM, a different flag sets the ''diagfi.nc'' output rate using ''outputs_per_sol''.

[[Category:Inputs]]
[[Category:WhatIs]]
[[Category:Generic-Model]]
[[Category:Generic-LMDZ]]
[[Category:Generic-1D]]
[[Category:Mars-Model]]
[[Category:Mars-LMDZ]]