Planets - User contributions [en]

Dynamico-giant

2024-03-06T17:30:01Z

LTeinturier: /* I) Appel à la physique */

[transferred from github, now active here]

= 1. Introduction =

"Welcome to the dynamico-giant wiki"

Dynamico-giant est un dépôt git qui permet d'installer le modèle Dynamico avec les fichiers de configs propres aux planètes géantes gazeuses et géantes glacées. Ci-dessous, vous trouverez presque (car nous ne sommes pas parfaits) toutes les informations sur l'installation, la compilation du modèle et les paramètres que vous pouvez modifier pour faire tourner le modèle.

"Allons-y!"

= 2. Installation =

== A) Fichiers d'architecture ==

=== - Modules ===

Before installation, set environment (do it once) in your .bash_profile.

Modules used on Occigen:

<pre>ulimit -s unlimited
# modules
module purge
module load intel/17.0
module load intelmpi/2017.0.098
module load hdf5/1.8.17
module load netcdf/4.4.0_fortran-4.4.2
module load ncview
module load nco
module load qt
module load python</pre>
Modules used on Ciclad:

<pre>ulimit -s unlimited
# modules
module purge
module load gnu/4.9.3
module load intel/15.0.6.233
module load openmpi/1.6.5-ifort
module load hdf5/1.8.18-parallel-ifort
module load netcdf4/4.4.1.1-parallel-ifort</pre>
Beware, if you change version of model (newer or older version), it's possible that you have to change modules...

=== - Installation sur un nouveau cluster ===

En pratique LMDZ.COMMON, ICOSA_LMDZ et IOIPSL peuvent utiliser exactement le même fichier arch.fcm ; mais celui pour ICOSAGCM est légèrement différent (les %FPP_DEF diffèrent, peut-être aussi le %FPP).

make_icosa_lmdz doit être lancé avec -arch_path ../ARCH puisque les arch.env et arch.path communs se trouvent dans ../ARCH (l'option -arch_path ne concerne d'ailleurs que les fichiers arch.env et arch.path; le fichier arch.fcm recherché sera toujours celui dans le "arch" de chacun des modèles. Donc il faut bien mettre pour chacun des quatre modèles (LMDZ.COMMON, ICOSA_LMDZ,IOIPSL et ICOSAGCM) le arch.fcm dans le sous-dossier arch/ du modèle correspondant.

L'installation d'IOIPSL doit se faire à partir du bon script bash. Si vous utilisez un server autre qu'occigen, il vous faut changer le nom du script à utiliser dans install_ioipsl.sh. Par exemple, il vous faudra utiliser install_ioipsl_ciclad-ifort.bash pour CICLAD.

XIOS est écrit en C++, et pas en Fortran. Le fichier arch.fcm correspondant est donc nécessairement différent de celui des quatre autres modèles. Pour créér ce arch.fcm, prendre exemple sur les fichiers .fcm déjà présent dans XIOS/arch/, avec une architecture similaire à celle du nouveau cluster. En particulier, il faut utiliser des versions de gcc/gfortran > 6+. Il faut absolument avoir une bibliothèque HDF5 compilée en parallèle, ainsi que netcdf-C et netcdf-fortran (et peut être aussi netcdf-Cxx) compilé en parallèle. Sans ça, il sera peut-être possible de compiler le modèle, mais pas de lancer une simulation en utilisant XIOS.

== B. Download ==

Download structure

Take place on a repository and git clone the model.

On occigen:

<pre>cd $SCRATCHDIR
git clone https://github.com/aymeric-spiga/dynamico-giant.git [optional different name]</pre>
== C. Install ==

Install code

<pre>cd dynamico-giant
./install.sh
./install_ioipsl.sh</pre>
A login and a password are necessary to install IOIPSL (only for the first time, after that we can save them). Please contact TGCC computing center to get the login/password.

After that, please add "$PWD"/FCM_V1.2/bin/ to PATH environment variable.

Ant that's all, we can now change parameters of files or compile the code to do a test (part 6) or eat a beautiful tartiflette (can be also do during compilation).

= 3. General parameters =

== A) Mesh grid & resolution ==

On run_icosa.def:

<pre>nbp --> number of subdivision on a main triangle: integer (default=40)
nbp = sqrt((nbr_lat x nbr_lon)/10)
nbp 20 40 80 160
T-edge length (km) 500 250 120 60
Example: nbp(128x96)=35 -> 40
nbp(256x192)=70 -> 80
nbp(360x720)=160 -> 160
nsplit_i, nsplit_j --> sub splitting of main rhombus: integer
Example: for nbp=80, nsplit_i=4,nsplit_j=6
nbp/nsplit_{i,j} = 20 > 10 & 13 > 10 --> GOOD</pre>
=== - Remapping ===

Pour le remapping, il faut aller dans context_lmdz_physics.xml changer les paramètres ni_glo et nj_glo par exemple pour remapper sur du lat/lon à 360 pts en latitude et 720 pts en longitude (0.5°).

<code><domain id="dom_regular" ni_glo="720" nj_glo="360" type="rectilinear"></code>

== B) Time ==

* ndays is number of days to simulate.
* day_step is number of dynamical time step per day
* Dynamics called every day_length(in s) / day_step per day
* 1 ts (physical timestep) = (day_length(in s)/day_step) x itau_physics
* Physics called (day_step / itau_physics) per day
* Physics called day_length(in s) / ts per day
* Radiative called every Physics x iradia physical timestep
* Radiative called every iradia/Physics days

== C) Sponge layer ==

<pre>iflag_sponge=0 for no sponge (default)
iflag_sponge=1 for sponge over 4 topmost layers
iflag_sponge=2 for sponge from top to ~1% of top layer pressure
mode_sponge=1 for u,v --> 0
mode_sponge=2 for u,v --> zonal mean
mode_sponge=3 for u,v,h --> zonal mean
tau_sponge --> damping frequency at last layer
--> e-5 medium / e-4 strong yet reasonable / e-3 very strong</pre>
Spiga et al (2020): Shaw and Shepherd (2007) showed that the inclusion of sponge-layer parameterizations that do not conserve angular momentum (which is the case for Rayleigh drag), or allow for momentum to escape to space, implies a sensitivity of the dynamical results (especially zonal wind speed) to the choice for model top or drag characteristic timescale, because of spurious downward influence when momentum conservation is violated.

Il ne faut donc pas ajouter de sponge layer pour les planètes géantes (iflag_sponge = 0)

== D) Dissipation ==

Spiga et al (2020): A subgrid-scale dissipation term is included in our Saturn DYNAMICO GCM to prevent the accumulation of energy at scales close to the grid resolution, caused by the GCM not resolving the turbulent scales at which this energy is dissipated. This hyperviscosity term is written in our Saturn DYNAMICO model as an iterated Laplacian term on a given variable. The three variables denoted are vorticity, divergence, and potential temperature, chosen to set horizontal dissipation on respectively the rotational component of the flow (e.g. Rossby waves), the divergent component of the flow (e.g. gravity waves), and the diabatic perturbations (e.g. coming from the physical packages).

<pre>tau_graddiv --> dissipation timescale of smallest wvl: u,v (gradiv) : real (default=5000)
tau_gradrot --> dissipation timescale of smallest wvl: u,v (nxgradrot) : real (default=5000)
tau_divgrad --> dissipation timescale of smallest wvl: h (divgrad) : real (default=5000)
nitergdiv --> number of iterations for gradiv operator : integer (default=1)
nitergrot --> number of iterations for nxgradrot operator : integer (default=1)
niterdivgrad --> number of iterations for divgrad operator : integer (default=1)</pre>
tau_graddiv,tau_gradrot,tau_divgrad correspondent au temps de dissipation (plus c'est petit, plus la dissipation est efficace). nitergdiv,nitergrot,niterdivgrad correspondent à l'ordre du laplacien (plus c'est grand, plus on dissipe sélectivement les petites échelles).

Attention, Trop dissiper -> instabilité numérique

Pas assez dissiper -> instabilité dynamique trop forte

L'ordre du laplacien dans les fichiers .def est suffisant. S'il faut changer la dissipation, il vaut mieux changer les temps de dissipation que l'ordre du laplacien (risque d'instabilité numérique).

=== - Facteurs de dissipation ===

A partir de ces valeurs de temps de dissipation, il est également possible d'ajouter 2 facteurs de dissipation (fac_mid et fac_up) qui permettent d'augmenter la dissipation à partir de certains niveaux. Le fac_mid permet de diminuer le temps de dissipation d'un facteur sur toute l'atmosphère jusqu'au bottom. Il y a une zone transition de quelques niveaux entre le bottom (le temp de dissipation au bottom est sans facteur) et quelques niveaux au-dessus où la valeur du temps de dissipation est avec ce facteur jusqu'au top. On ne peut donc pas choisir l'altitude et la zone de transition avec ce facteur. Quant au fac_up, il permet de diminuer le temps de dissipation d'un facteur à partir d'une altitude et d'une épaisseur de zone de transition qu'on peut choisir. On peut utiliser la combinaison des 2 ou l'un des 2. Pour ne pas prendre en compte ces facteurs de dissipation, il faut qu'il soit mis à 1.

Il existe 2 modes possibles pour le fac_up: le mode martien (en fonction de l'altitude) et le mode vénusien (en fonction de la pression). Dans le 1er cas, il faut renseigner l'altitude où démarre la zone de transition et l'épaisseur de cette zone. Dans le 2ème cas, il faut renseigner la pression où démarre la zone de transition et l'échelle de hauteur de cette zone de transition.

Il est vivement conseillé de tracer votre profil de dissipation à partir des équations (dans vert_prof_dissip_icosa_lmdz.f90) afin d'être certain de leur allure.

== E) Rayleigh Friction ==

Spiga et al 2020: This drag plays the role devoted to surface friction on terrestrial planets, which allows to close the angular momentum budget through downward control (Haynes and McIntyre, 1987; Haynes et al., 1991). This could also be regarded as a zeroth-order parameterization for Magneto-HydroDynamic (MHD) drag as a result of Lorenz forces acting on jet streams putatively extending to the depths of Saturn's interior (Liu et al., 2008; Galanti et al., 2019), much deeper than the shallow GCM's model bottom.

Comme Liu and Schneider (2010), cette couche de frottement ne s'exerce pas aux régions équatoriales car à l’extérieur du cylindre tangeant du rayon équatorial, les cylindres convectifs ne coupent pas cette couche et ne subissent pas les frottements liés aux effets de la MHD.

Le paramètre rayleigh_limlat correspond à la latitude maximale de part et d'autre de l'équateur où cette friction n'est pas observée.

Pour le temps (rayleigh_friction_tau), il est similaire à celui utilisé par Liu and Schneider (2010). Elle est fixée à 100 jours terrestres. Mais attention, on fonction de la valeur, la dynamique peut évoluer dans votre simulation.

== F) Conservation ==

Dans run_icosa.def:

<pre>check_conservation = detailed
itau_check_conserv = 320</pre>
Le itau_check_conserv est comme itau_physics pour que la contribution de la physique à AAM ne soit pas nulle. de plus il est coûteux d'appeler les diagnostics trop souvent.

S. Lebonnois, C. Covey, A. Grossman, H. Parish, G. Schubert, R. Walterscheid, P. Lauritzen, and C. Jablonowski. Angular momentum budget in General Circulation Models of superrotating atmospheres: A critical diagnostic. Journal of Geophysical Research (Planets), 117:E12004, 2012.

P. H. Lauritzen, J. T. Bacmeister, T. Dubos, S. Lebonnois, and M. A. Taylor. Held-Suarez simulations with the Community Atmosphere Model Spectral Element (CAM-SE) dynamical core: A global axial angular momentum analysis using Eulerian and floating Lagrangian vertical coordinates. Journal of Advances in Modeling Earth Systems, 6:129-140, 2014.

== G) Traceur ==

=== - Branche master version du 04/04/2022 dans le dossier Jupiter (traceur) ===

La version actuelle des fichiers de réglage dans le dossier Jupiter est réglée pour utiliser des traceurs. Si vous souhaiter utiliser le modèle sans traceur, il vous faut modifier les fichiers suivants :

run_icosa.def

<pre> nqtot = 0</pre>
traceur.def

<pre> 0</pre>
context_dynamico.xml (ligne 125)

<pre> <field id="q_start" name="q" grid_ref="grid_q_start" prec="8"/> </pre>
== H) Bottom du modèle ==

* refaire tourner le 1D
** changer la pression <code>psurf</code> dans rcm1d.def
** changer <code>ichoice=1</code> et changer <code>tref</code> (ex: 10b: 330K)
* refaire tourner le 3D (avec makestart) en changeant <code>preff</code> dans <code>[nom_de_la_planète]_const.def</code>

= 4. Physical parameterizations =

== A) Cycle diurne ==

<code>diurnal = .false.</code> Car le temps radiatif est beaucoup plus long sur les planètes géantes.

== B) Anneaux ==

[Pour saturne uniquement] Il est possible d'ajouter l'ombre des anneaux dans le callphys.def:

<code>rings_shadow = .true.</code>

== C) Collision-induced absorption data ==

L'absorption induite par collision est importante dans le transfet radiatif sur les planètes géantes, notamment celle provenant de H2. Cependant, le rapport ortho-para du H2 est différent selon les planètes et peut modifier fortement le chauffage sur ces planètes. Pour le cas des géantes gazeuses, le rapport ortho-para est normal (rapport 3:1). Mais pour le cas des géantes glacées, le rapport ortho-para est à l'équilibre. Par défaut, le rapport ortho-para est normal (H2orthopara_mixture = normal). Pour les géantes glacées, il faut mettre H2orthopara_mixture = equilibrium .

Si vous voulez modifier le fichier utilisé dans un cas de CIA, il faut aller dans le code (dynamico-giant/code/LMDZ.GENERIC/libf/phystd/interpolate????.F90) et changer le nom du fichier (assurez-vous que ce fichier soit dans votre répertoire DATAGENERIC/continuum_data).

Dans la version actuelle, il n'existe pas d'option qui permet de ne pas utiliser certaines contributions de CIA. Par exemple, si votre atmosphère est composé de H2, He et CH4 et que vous ne voulez pas des contributions venant du CH4, il vous faut commenter dans le modèle ces contributions.

== D) K-correlated data ==

Les fichiers k-corrélées doivent se situer dans votre répertoire DATAGENERIC/corrk_data .

== E) Cpp mode ==

cpp_mugaz_mode=0 pour que la valeur de cpp et de mugaz proviennent de la dynamique. cpp_mugaz_mode=1 pour forcer la valeur de cpp et de mugaz dans le callphys.def (à utiliser dans le makestart uniquement) cpp_mugaz_mode=2 pour calculer automatiquement à partir de données de références à 300 K et du gases.def (à éviter)

== F) Generic n-layer aerosols (replaces the former 2-layer and NH3 cloud) ==

Ce mode permet de créer des couches d'aérosols/nuages ayant une opacité (aeronlay_tauref) à une longueur d'onde donnée (aeronlay_lamref), un rayon de particule fixe (aeronlay_size) situés à une pression fixe et en fonction des propriétés optiques des particules (optprop_aeronlay_vis et optprop_aeronlay_ir). Pour le cas de l'altitude, on peut choisir soit (aeronlay_choice = 1) une pression max (aeronlay_ptop) et une pression min (aeronlay_pbot), soit (aeronlay_choice = 2) une pression min (aeronlay_pbot) et une échelle de hauteur (aeronlay_sclhght). On peut également choisir la variance effective pour les rayons de particules.

Le nombre de couche de nuage (nlayero) est égal au nombre de scatterers dans la compilation (option -s ).

Les fichiers de propriétés optiques doivent être dans votre DATAGENERIC.

Un exemple pour 3 couches:

<pre>aeronlay = .true.
nlayaero = 3
aeronlay_tauref = 2.5 0.04 0.1
aeronlay_lamref = 0.8e-6 0.8e-6 0.16e-6
aeronlay_choice = 2 2 1
aeronlay_pbot = 1.5e5 1.6e5 20.
aeronlay_ptop = 1.1e5 1. 1.
aeronlay_sclhght = 0.1 2.0 1
aeronlay_size = 0.5e-6 0.05e-6 0.5e-6
aeronlay_nueff = 0.3 0.3 0.3
optprop_aeronlay_vis = optprop_aerosol2_vis.dat optprop_aerosol3_vis.dat optprop_carbon4_vis.dat
optprop_aeronlay_ir = optprop_aerosol2_ir.dat optprop_aerosol3_ir.dat optprop_carbon4_ir.dat</pre>
== G) Panaches thermiques ==

# Comme pour toutes nouvelles simulations, il faut commencer par un run 1D de plusieurs décennies permettant d'obtenir un profil de température (<code>temp_profile.txt</code>) et les coefficients ap et bp (<code>apbp.txt</code>) en équilibre radiatif-convectif pour la planète que l'on étudie. Ce run 1D s'effectue dans les dossiers jupiter1d, saturn1d, neptune1d et/ou uranus1d, avec ces options pour le callphys.def :

callrad = true calladj = true tous les autres mots clés à false (y compris <code>calltherm</code>).

<ol start="2" style="list-style-type: decimal;">
<li>Faire un "makestart" run : permet d'obtenir les fichiers restart à partir du profil de température initial. Il s'agit du run dont l'état initial est le profil de température créé par run 1D (<code>temp_profile.txt</code> est appliquer à chaque point de grille horizontale du modèle). L'état initial est ainsi une planète isotherme horizontalement mais qui varie verticalement. Pour ce run, aucun traceur ne va être utilisé dans le modèle. Néanmoins, il faut renseigner au modèle le nombre de traceurs que nous souhaitons utiliser avec le schéma des panaches thermiques afin qu'il puisse créer la dimension nq et le champ q dans les fichiers <code>restart_icosa.nc</code> et <code>restartfi.nc</code>. dans '''callphys.def''' :
<pre> traceur = true</pre>
dans '''run_icosa.def''' : <code>nqtot = 2</code> (par exemple, h2o_vap et h2o_ice) dans '''traceur.def''' :
<pre> 2
h2o_vap
h2o_ice</pre></li>
<li>Ajout des quantités pour chaque traceur. Ici, nous ajoutons les profils pour chacun des traceurs dans les fichiers <code>restart_icosa.nc</code> et <code>restartfi.nc</code> "à la main" en utilisant le programme python <code>/processing_codes/tracer_settings.py</code> pour obtenir des fichiers restart avec la bonne abondance d'eau dans le cas présent.</li>
<li>Commencer la simulation : Il ne reste plus qu'à lancer la simulation 3D avec les nouveaux fichiers restart et les réglages suivant : dans '''callphys.def''' :
<pre> traceur = true
calltherm = true
# thermal plume model options:
divmpl = true
r_aspect_thermals = 2.0
tau_thermals = 0.0
betalpha = 0.9
afact = 0.7
fact_epsilon = 2.e-4
alpha_max = 0.7
fomass_max = 0.5
pres_limit = 2.e5
water = true
watercond = true
waterain = true
evap_prec = true</pre></li></ol>

dans '''run_icosa.def''' : <code>nqtot = 2</code> dans '''traceur.def''' :

<pre> 2
h2o_vap
h2o_ice</pre>
Enfin, ajouter la déclaration et l'écriture des variables relatives à l'utilisation du schéma des thermiques (<code>h2o_vap</code>, <code>h2o_ice</code>, <code>w_plm</code>) dans les fichiers XML de la physique.

== H) Rotation ==

Dans [name_of_planet]_const.def, changer le taux de rotation.

Cela n'a jamais été testé avec 0, mais cela pourrait créer des problèmes (ex: beta). ''vérifier que omega dans saturn_const.def n'intervient pas dans la physique LMDZ.GENERIC/libf/phystd/''

== I) Appel à la physique ==

Il faut changer à ''deux endroits'' en réglant la même valeur: - dans run_icosa.def, changer itau_physics ; - dans run.def, changer iphysiq. Ces paramètres sont exprimés en pas de temps dynamique. Pour appeler la physique à chaque pas de temps dynamique, régler ces paramètres à 1.

'''Attention: ''' Lorsque l'on change dans le ''run_icosa.def'' la valeur de itau_physics, faire bien attention à la valeur du champ itau_adv (qui controle la frequence d'advection des traceurs par la dynamique, compté en pas de temps dynamique). Il faut que itau_physics soit un multiple de itau_adv. Cela veut dire, que l'on advecte notre champ de traceurs plus souvent (ou autant de fois) que l'on appelle la physique, qui peut transformer ces traceurs.

= 5. Ecriture des fichiers =

== A) Nombre d'écriture dans chaque Xhistins.nc ==

Prenons l'exemple de l'écriture d'un fichier sur une simu Uranus pour 1000j

Nombre d'écriture = day_step X ndays/(itau_physics X output_freq)

* ndays = 1000
* day_step = 200
* itau_physics = 50
* output_freq = 200

On a 20 sorties tous les 1000 jours.

Globalement, il est conseillé d'avoir 800 à 1000 sorties minimum par année planétaire (voir beaucoup plus en fonction de votre étude).

== B) Ajouter une variable ==

Pour ajouter une variable dans le fichier de sortie Xhistins.nc, il faut l'ajouter dans le field_group correspondant et s'assurer que cette variable est présente en tant que variable de sortie dans dynamico-giant/code/LMDZ.GENERIC/libf/phystd/physiq_mod.F90. Si ce n'est pas le cas, il faut faire un call writediagfi de la variable puis l'ajouter dans un call send_xios_field.

== C) Forcer XIOS à écrire tous les .... ==

Il faut utiliser l'attribut (ici par ex pour forcer l'écriture tous les ts (==time step) ; d'autres délais doivent être possible): <code>sync_freq="1ts"</code> dans le <code><file id=... .... ></code> concerné.

== D) XIOS server ou client ==

# XIOS client <code>use_server=False</code> Broadwell 24 processeurs sur 28
# XIOS server <code>use_server=true</code> Broadwell 24 processurs sur 28 + 4 processeurs pour XIOS

Le cas 1 est 3 min plus lent que le cas 2 -- sur 2h20...(!) parce qu'on fait peu de sorties.

= 6. Running the model =

== A) Compilation ==

Pour compiler le modèle, il suffit d'aller dans le dossier avec vos fichiers .def et .xml et d'utiliser le script compile_occigen.sh si on est sur occigen ou compile_ciclad.sh si on est sur ciclad

<code>./compile_occigen.sh</code>

Les options qui peuvent être modifiées/ajoutées sont - -s : le nombre de scatterers - -parallel: mpi ou mpi_omp - -arch: le nom du ficher d'architecture - -arch_path: le chemin vers ce fichier - -job: 8 (conseillé) - -full: compile le code entièrement - -debug: pour trouver un éventuel bug. A ENLEVER obligatoirement s'il n'y a plus de bug car multiplie par 5 votre temp de calcul.

En 1D, d'autres options sont nécessaires: - -t: nombre de traceurs - -d: nombre de niveaux verticaux - -parallel: none (en 1D uniquement)

=== - Versions fonctionnelles de dynamico-giant ===

'''Jupiter''' - XIOS revision 1583 
- ICOSAGCM revision 756 
- IOIPSL revision 339 
- FCM_V1.2 revision 12 
- Physics revision 2228 (LMDZ.COMMON, LMDZ.GENERIC, ICOSA_LMDZ, ARCH)

but check Physics version 2142 and 2180

'''Saturne (référence Spiga 2020, à vérifier)''' - DYNAMICO --> Revision: 756 - PHYSICS --> Revision: 2005 - XIOS --> Revision: 1583 - IOIPSL --> Revision: 310?

'''Saturne (simulation de référence sur 61 niveaux, Bardet et al. Icarus 2021)''' - XIOS revision 1583 
- ICOSAGCM revision 765 
- IOIPSL revision 310 
- FCM_V1.2 revision 12 
- Physics revision 2005 (LMDZ.COMMON, LMDZ.GENERIC, ICOSA_LMDZ, ARCH)

'''Saturne (simulation avec la GWD paramétrisation sur 61 niveaux, chapitre 6 PhD Bardet)''' - XIOS revision 1583 
- ICOSAGCM revision 765 
- IOIPSL revision 310 
- FCM_V1.2 revision 12 
- Physics revision 2213 (LMDZ.COMMON, LMDZ.GENERIC, ICOSA_LMDZ, ARCH)

'''Saturne (simulation sur 96 niveaux, Bardet et al. Nature Astronomy 2022)''' - XIOS revision 1583 
- ICOSAGCM revision 765 
- IOIPSL revision 431 
- FCM_V1.2 revision 12 
- Physics revision 2305 (LMDZ.COMMON, LMDZ.GENERIC, ICOSA_LMDZ, ARCH)

'''Saturne (simulation avec la GWD paramétrisation sur 96 niveaux, chapitre 6 PhD Bardet)''' - XIOS revision 1583 
- ICOSAGCM revision 765 
- IOIPSL revision 431 
- FCM_V1.2 revision 12 
- Physics revision 2403 (LMDZ.COMMON, LMDZ.GENERIC, ICOSA_LMDZ, ARCH)

'''Uranus & Neptune (old version)''' - XIOS revision 1944 - ICOSAGCM revision 765 - IOIPSL revision 431 - FCM_V1.2 revision 12 
- Physics revision 2413 (LMDZ.COMMON, LMDZ.GENERIC, ICOSA_LMDZ, ARCH)

'''Uranus & Neptune (new version)''' - XIOS revision 2203 - ICOSAGCM revision (20/08/2021) - IOIPSL revision 450 - FCM_V1.2 revision 12 
- Physics revision 2555 (LMDZ.COMMON, LMDZ.GENERIC, ICOSA_LMDZ, ARCH)

'''Jupiter, Saturne, Uranus & Neptune [OCCIGEN VERSION]''' - XIOS revision 2319 - ICOSAGCM revision (90f7138a60ebd3644fbbc42bc9dfa22923386385) - IOIPSL revision 453 - FCM_V1.2 revision 12 
- Physics revision 2655 (LMDZ.COMMON, LMDZ.GENERIC, ICOSA_LMDZ, ARCH)

'''Jupiter, Saturne, Uranus & Neptune [IRENE VERSION]''' - XIOS revision 2399 - ICOSAGCM revision (4fbd393a9051fd9c1a5b662683f6ad8d0dc2867c) - IOIPSL revision 6234 - FCM_V1.2 revision 12 
- Physics revision 2842 (LMDZ.COMMON, LMDZ.GENERIC, ICOSA_LMDZ, ARCH)

== B) Run 1D ==

Pour une simulation 1D, 1 seul CPU suffit.

To run on Occigen:

<code>sbatch job_mpi</code>

To run on Irene:

<code>ccc_msub job_mpi</code>

To see evolution of your job on Occigen:

<code>squeue -u username</code>

To see evolution of your job on Irene:

<code>ccc_mpp -u username</code>

To kill your job on Occigen:

<code>scancel id_of_your_job</code>

To kill your job on Irene:

<code>ccc_mdel id_of_your_job</code>

Une fois que votre simulation est fini, lancer le script python printapbp.py qui vous permettra d'obtenir les fichiers apbp.txt et temp_profile.txt nécessaires pour les simulations en 3D.

== C) Run 3D ==

Pour lancer une première simulation, il faut tout d'abord le faire une seule fois dans le dossier makestart de l'un des dossiers jupiter, saturne, uranus ou neptune. Ce dernier permettra de créer les premiers fichiers restart_icosa.nc et restartfi.nc. Ensuite, il suffit revenir dans le dossier parent (jupiter, saturne, uranus, neptune) et de lancer ses simulations.

=== - MPI ===

Avant de lancer une simulation, il est nécessaire de déterminer - le nombre de procs - le nombre de noeuds

Il faut avant tout savoir qu'il y a 10 x nsplit_i x nsplit_j domaines. Ces domaines sont divisés en plusieurs cellules. Il y a (nbp)²/(nsplit_i x nsplit_j) cellules par domaine. Ce nombre est égal au nombre de MPI process.

Il est important de ne pas faire des tuiles "trop petites" et donc de garder nsplit_i et nsplit_j tels que nbp/nsplit_{i,j} >=10-15 (plutôt 15). Car il y a des calculs "redondants" faits sur les bords du domaine. Admettons que ton domaine soit 10 x 10 (100 cellules par domaine). Le nombre de cellules au bord du domaine est égal à 2 x nsplit_i + 2 x nsplit_j - 4. Dans ce cas, le nombre de cellules au bord est de 36. Cela signifie que le bord du domaine correspond à 36% du domaine (36 points sur 100). Si le domaine est 15 x 15, alors le bord du domaine n'est plus que de 24.8 %. On aimerait avoir des domaines les plus grands possibles mais il faut aussi voir que chaque cellule correspond à une colonne à résoudre... Et c'est là ou il faut expérimenter un peu pour trouver l'optimum entre le nombre total de proc à employer et le gain effectif (temps total "facturé" au vu du nombre de procs monopolisés pour une simu donnée). Ehouarn

Il faut au total que le nombre de procs soit égal à 10 x nsplit_i x nsplit_j qu'il faut répartir entre les noeuds.

Pour les noeuds, il faut tenir compte qu'un noeud sur Occigen, c'est 24 (ou 28) procs. Disons 24. Quand on demande N procs, le système te donne M noeuds, soit M X 24 procs (les noeuds ne sont pas partagés avec d'autres applications). Et bien sûr on te facturera ces M noeuds, même si tu n'utilises pas tous les procs. Donc il faut s'efforcer de tomber juste et demander un multiple de 24 procs. Même raisonement si tu demandes à utiliser des noeuds de 28 procs.

Par exemple, nsplit_i = 4 et nsplit_j = 6. 4 x 6 est égal à 24. On utilisera donc 10 noeuds pour avoir 240 nprocs.

Prenons un autre exemple. nsplit_i = 5 et nsplit_j = 6. 5x6 = 30 Or il y a 24 procs par noeuds. Il faudra utiliser 13 noeuds. Et ce ne sont pas 300 procs qui seront facturés au total mais 312.

Mais comme vu plus haut, en fonction de la configuration de tes domaines, il est possible que le 2nd exemple coûte moins cher en heure CPU que le 1er exemple car il sera plus rapide.

=== - OpenMP / MPI ===

Dans le cas mixte OpenMP/MPI, il est nécessaire de renseigner le nombre de tâche OpenMP et de MPI threads. Dans le script pour le job, on a :

cpus-per-task et OMP_NUM_THREADS correspondent au nombre de tâche OpenMP ntask-per-node correspond au nombre de tâche par noeud

Pour déterminer la configuration à utiliser en OpenMP/MPI, il faut respecter ces règles: - ntask-per-node x nodes correspond au nombre total de MPI threads qu'il devrait y avoir (c'est-à-dire dans le cas MPI seul). - cpus-per-task x ntask-per-node doit être égale ou inférieur au nombre de procs par noeuds (24 ou 28- sur Occigen) - cpus-per-task devrait être ~égal au nombre de niveaux verticaux / 10 - OMP_NUM_THREADS doit être égal au cpus-per-task - omp_level_size (dans le fichier run_icosa.def) devrait être ~égal au nombre de cpus-per-task. (dans le cas MPI seul, il faut le laisser égal à 1)

Au total, le nombre total de nprocs est égal à cpus-per-task x ntask-per-node x nodes

En reprenant le premier exemple MPI (on suppose qu'on a 40 niveaux verticaux), on aura - cpus-per-task = 4 - ntask-per-norde = 6 On a bien cpus-per-task x ntask-per-norde <= 24. Le nombre de noeuds doit être égal à 40 car 6x40 = 240 nprocs MPI Le nombre de procs qu'on utilisera au total sera de 960.

=== - Comment déterminer le nombre d'heure CPU (consommation en heure)? ===

Heures CPU = durée_de_la_simu X nombre de procs total facturé

Si on veut connaître le nombre d'heure CPU qu'on consommera à partir d'un test, voici le calcul

Heures CPU = durée_de_la_simu_test X nombre de procs total facturé X nombre de jour qu'on veut simulé / (nombre de jour simulées pendant la simu test)

== D) Optimisation ==

Pour savoir quelle est la meilleure configuration en MPI seul, en OpenMP/MPI et entre les 2, des tests ont été effectués. Ces tests ont été réalisés sur Uranus à nbp80 sur les noeuds Broadwell (28 procs par noeud). 1000 jours uraniens ont été simulés à chaque fois.

'''MPI:''' [[File:fig1.png]]

'''OpenMP/MPI:''' [[File:fig2.png]]

La comparaison du MPI seul et du mixte OpenMP/MPI permet de tirer la conclusion que le mixte OpenMP/MPI est beaucoup plus rapide (facteur 2.5 à 3) mais consomme plus (facteur 1.1 à 1.7).

[[File:fig3.png]]

== E) Visualizing the output files ==

=== - Comment comparer deux fichiers? file1.nc file2.nc ===

ncdiff file1.nc file2.nc output.nc

=== - Comment coller plusieurs fichiers file1.nc file2.nc file3.nc ===

ncrcat file1.nc file2.nc file3.nc output.nc

Avec juste la vitesse u et v: ncrcat -v u -v v file1.nc file2.nc file3.nc output.nc

Avec tous les fichiers fileXXX.nc disponibles ncrcat file*.nc output.nc

=== - Comment tracer les champs et divers diagnostics en moyenne zonale? ===

<ul>
<li>mettre à jour planetoplot et planets</li>
<li>utiliser precast.py dans planetoplot/examples/ppclass_additional/dynanalysis</li>
<li>changer les options au début (il y a des exemples), exemple pour Saturne
<source lang="python">fileAP="Xhistins_42.nc"
p_upper,p_lower,nlev = 4.0e2,2.5e5,40
targetp1d = np.logspace(np.log10(p_lower),np.log10(p_upper),nlev)
myp = planets.Saturn
day_per_year = 24430.
short = False
includels = False
charx = "0,360"
ispressure = False
vartemp = "temperature"
outfile = "precast.nc"
nopole = True</source></li>
<li>Il faut que apbp.txt soit présent !</li>
<li>Le résultat se trouve dans le fichier indiqué dans outfile.</li></ul>

== F) Debug ==

=== - Pour debugger ===

Réglez <code>info_level</code> à <code>100</code> dans le fichier <code>iodef.xml</code> .

=== - Erreur en début de run, SEGMENTATION FAULT dans la routine <code>advect.f90</code> (ICOSAGCM/ppsrc/transport) au niveau de l'appel à <code>cross_product2</code> ===

La routine <code>cross_product2</code> effectue des produits vectoriels. Cette erreur vient du fait qu'en essayant de faire des optimisations de calcul, le compilateur se prend "les pieds dans le tapis" et ne sait plus faire le produit vectoriel. Pour empêcher le compilateur de faire des optimisations trop poussées, il faut ajouter l'option <code>-fp-model strict</code> au niveau du mot-clé <code>%BASE_FFLAGS</code> dans tous les fichiers <code>.fcm</code> du modèle (LMDZ.COMMON / ICOSAGCM / ICOSA_LMDZ).

== G) Pas de radiatif et perturbations ==

<ul>
<li>faire descendre le modèle</li>
<li>modifier le code pour perturbations vitesse et pas température</li>
<li>compiler (voir README.md)</li>
<li>changer dans callphys.def (dans saturn/makestart)
<pre>surfalbedo = 1.0
surfemis = 0.0</pre></li>
<li>changer dans callphys.def (dans saturn '''et''' saturn/makestart)
<pre>corrk = .false.
enertest = .true.
randompert = 1</pre></li>
<li>voir si on garde le flux interne ou non (<code>intheat</code>)</li>
<li>refaire les états initiaux (dans makestart)</li>
<li>faire un run court pour voir sur les spectres l'injection d'Ek</li></ul>

== H) PROFILING with DYNAMICO-Giant ==

<code>-pg</code>: Generate extra code to write profile information suitable for analysis program <code>gprof</code>. You must use this option when compiling the sources files you want data about and you must '''also use it when linking'''. 1. For DYNAMICO-Giant: in the arch.fcm file of each part of the model (ICOSAGCM, ICOSA_LMDZ, LMDZ.COMMON and XIOS), you have to add the option <code>-pg</code> to <code>%PROD_FFLAGS</code> and <code>%BASE_LD</code> 2. Compile as usual 3. Execute your code as usual 4. Run <code>gprof</code> tool:

<pre> gprof icosa_lmdz.exe goon.out > analysis.txt</pre>

= OLD WIKI (DEPRECATED, JUST FOR REFERENCE/INFORMATION) =

``Welcome to the dynamico-giant wiki!

== pour debugger ==

régler <code>info_level</code> à <code>100</code> dans le fichier <code>iodef.xml</code>

== pourquoi diurnal=.false. ==

le temps radiatif est beaucoup plus long sur les géantes

== comment changer la fréquence d'appel à la physique? ==

il faut changer à ''deux endroits'' en réglant la même valeur - dans run_icosa.def, changer itau_physics - dans run.def, changer iphysiq ces paramètres sont exprimés en pas de temps dynamique. pour appeler la physique à chaque pas de temps dynamique, régler ces paramètres à 1

== comment changer la résolution ==

[https://github.com/aymeric-spiga/dynamico-giant/commit/f39e3b651c19cefb26da72ab77933520ff8f27f5 voir lien ici]

== comment choisir le nombre de processeurs pour la résolution ==

voir commentaires sur run_icosa.def

<pre>###########################################
## There must be less MPIxOpenMP processes than the 10 x nsplit_i x nsplit_j tiles
## typically for pure MPI runs, let nproc = 10 x nsplit_i x nsplit_j
## it is better to have nbp/split >~ 10
###########################################
#### 40 noeuds de 24 processeurs = 960 procs
nsplit_i=12
nsplit_j=8</pre>
Ehouarn: ''Une règle à suivre est de ne pas faire des tuiles "trop petites" et donc de garder nbsplit_i et nbsplit_j tels que nbp/nbsplit >=10-15 (plutôt 15). Car il y a des calculs "redondants" faits sur les bords du domaine. Or si ton domaine est 10''10, le bord du domaine c'est 38% du domaine (38 points sur le bord pour 100 points en tout) alors que si ton domaine est 15''15, le bord du domaine n'est plus que 56/(15''15) = 25.7%. 56=2''15+2''13 (pour ne pas compter 2 fois les coins). On voudrait du coup des domaine les plus grand possible, clairement, mais il faut aussi voir que chaque domaine c'est nbp''nbp colonnes à résoudre sur ce même proc... Et c'est là ou il faut expérimenter un peu pour trouver l'optimum entre le nombre total de proc à employer et le gain effectif (temps total "facturé" au vu du nombre de procs monopolisés pour une simu donnée).''

_Pour les noeuds, il faut tenir compte qu'un noeud sur Occigen, c'est 24 (ou 28) procs. Disons 24. Quand on demande N procs, le système te donne M noeuds, soit M*24 procs (les noeuds ne sont pas partagés avec d'autres applications). Et bien sûr on te facturera ces M noeuds, même si tu n'utilises pas tous les procs. Donc il faut s'efforcer de tomber juste et demander un multiple de 24 procs. Même raisonement si tu demandes à utiliser des noeuds de 28 procs._

== comment changer la rotation? ==

dans saturn_const.def, changer le taux de rotation jamais testé avec 0, mais pourrait créer des problèmes (ex: beta) ''vérifier que omega dans saturn_const.def n'intervient pas dans la physique LMDZ.GENERIC/libf/phystd/''

== la dissipation c'est où? ==

laplacien itéré: hyperviscosity - tau_graddiv,tau_gradrot,tau_divgrad: periode (plus c'est petit, plus la dissipation est efficace) - nitergdiv,nitergrot,niterdivgrad: ordre du laplacien (plus c'est grand, plus on dissipe sélectivement les petites échelles)

== comment changer la grille cible du remapping on-the-fly ==

dans context_lmdz_physics.xml changer les paramètres ni_glo et nj_glo par exemple pour remapper sur du lat/lon à 360 pts en latitude et 720 pts en longitude (0.5°) <domain id="dom_regular" ni_glo="720" nj_glo="360" type="rectilinear">

== Forcer XIOS à écrire tous les .... ==

Il faut utiliser l'attribut (ici par ex pour forcer l'écriture tous les ts (==time step) ; d'autres délais doivent être possible): <code>sync_freq="1ts"</code> dans le <code><file id=... .... ></code> concerné.

== comment changer le bottom du modele ==

* refaire tourner le 1D
** changer la pression <code>psurf</code> dans rcm1d.def
** changer <code>ichoice=1</code> et changer <code>tref</code> (ex: 10b: 330K)
* refaire tourner le 3D (avec makestart) en changeant <code>preff</code> dans <code>saturn_const.def</code>

== comment changer la fréquence de sortie des Xhistins.nc ==

les commandes XIOS sont appelées depuis la physique (dans la config présente) dans context_lmdz_physics.xml il suffit de changer output_freq ''attention'' ts se comprend comme le pas de temps physique (voir donc dans run_icosa.def les paramètres dt et itau_physics pour le connaître) par exemple, pour des sorties tous les 20 jours Saturne

<source lang="xml"><file id="histins"
name="Xhistins"
output_freq="40ts"
type="one_file"
enabled=".true."></source>
== pas de radiatif et perturbations ==

<ul>
<li>faire descendre le modèle</li>
<li>modifier le code pour perturbations vitesse et pas température</li>
<li>compiler (voir README.md)</li>
<li>changer dans callphys.def (dans saturn/makestart)
<pre>surfalbedo = 1.0
surfemis = 0.0</pre></li>
<li>changer dans callphys.def (dans saturn '''et''' saturn/makestart)
<pre>corrk = .false.
enertest = .true.
randompert = 1</pre></li>
<li>voir si on garde le flux interne ou non (<code>intheat</code>)</li>
<li>refaire les états initiaux (dans makestart)</li>
<li>faire un run court pour voir sur les spectres l'injection d'Ek</li></ul>

== comment comparer deux fichiers? file1.nc file2.nc ==

ncdiff file1.nc file2.nc output.nc

== comment coller plusieurs fichiers file1.nc file2.nc file3.nc ==

ncrcat file1.nc file2.nc file3.nc output.nc avec juste la vitesse u et v ncrcat -v u -v v file1.nc file2.nc file3.nc output.nc avec tous les fichiers fileXXX.nc disponibles ncrcat file*.nc output.nc

== comment tracer les champs et divers diagnostics en moyenne zonale? ==

<ul>
<li>mettre à jour planetoplot et planets</li>
<li>utiliser precast.py dans planetoplot/examples/ppclass_additional/dynanalysis</li>
<li>changer les options au début (il y a des exemples), exemple pour Saturne
<source lang="python">fileAP="Xhistins_42.nc"
p_upper,p_lower,nlev = 4.0e2,2.5e5,40
targetp1d = np.logspace(np.log10(p_lower),np.log10(p_upper),nlev)
myp = planets.Saturn
day_per_year = 24430.
short = False
includels = False
charx = "0,360"
ispressure = False
vartemp = "temperature"
outfile = "precast.nc"
nopole = True</source></li>
<li>Il faut que apbp.txt soit présent !</li>
<li>le résultat se trouve dans le fichier indiqué dans outfile</li></ul>

== bilan de moment cinétique ==

dans run_icosa.def

<pre>check_conservation = detailed
itau_check_conserv = 160</pre>
le itau_check_conserv est comme itau_physics pour que la contribution de la physique à AAM ne soit pas nulle. de plus il est coûteux d'appeler les diagnostics trop souvent.

S. Lebonnois, C. Covey, A. Grossman, H. Parish, G. Schubert, R. Walterscheid, P. Lauritzen, and C. Jablonowski. Angular momentum budget in General Circulation Models of superrotating atmospheres: A critical diagnostic. Journal of Geophysical Research (Planets), 117:E12004, 2012.

P. H. Lauritzen, J. T. Bacmeister, T. Dubos, S. Lebonnois, and M. A. Taylor. Held-Suarez simulations with the Community Atmosphere Model Spectral Element (CAM-SE) dynamical core: A global axial angular momentum analysis using Eulerian and floating Lagrangian vertical coordinates. Journal of Advances in Modeling Earth Systems, 6:129-140, 2014.

== XIOS server ou client ==

# XIOS client <code>use_server=False</code> Broadwell 24 processeurs sur 28
# XIOS server <code>use_server=true</code> Broadwell 24 processurs sur 28 + 4 processeurs pour XIOS

cas 1 est 3 min plus lent que cas 2 -- sur 2h20...! parce qu'on fait peu de sorties

== Optimisation ==

# 40 noeuds au lieu de 50 noeuds, même temps de calcul...! Testé dans la branche https://github.com/aymeric-spiga/dynamico-giant/tree/work_61levels
# <code>sync_freq=40ts</code> (ou multiple de, si 40ts est la fréquence d'écriture) dans context_lmdz_physics.xml doit être inclus sinon il n'est pas réglé et XIOS transfère tout à la fin du run ce qui prend du temps (et peut occasionner des problèmes de mémoire)

== Versions fonctionnelles de dynamico-giant ==

=== Jupiter ===

* XIOS revision 1583 

* ICOSAGCM revision 756 

* IOIPSL revision 339 

* FCM_V1.2 revision 12 

* Physics revision 2228 (LMDZ.COMMON, LMDZ.GENERIC, ICOSA_LMDZ, ARCH)

but check Physics version 2142 and 2180

=== Saturne (référence Spiga 2020, à vérifier) ===

* DYNAMICO --> Revision: 756
* PHYSICS --> Revision: 2005
* XIOS --> Revision: 1583
* IOIPSL --> Revision: 310?

=== Saturne (simulation de référence sur 61 niveaux, Bardet et al. Icarus 2021) ===

* XIOS revision 1583 

* ICOSAGCM revision 765 

* IOIPSL revision 310 

* FCM_V1.2 revision 12 

* Physics revision 2005 (LMDZ.COMMON, LMDZ.GENERIC, ICOSA_LMDZ, ARCH)

=== Saturne (simulation avec la GWD paramétrisation sur 61 niveaux, chapitre 6 PhD Bardet) ===

* XIOS revision 1583 

* ICOSAGCM revision 765 

* IOIPSL revision 310 

* FCM_V1.2 revision 12 

* Physics revision 2213 (LMDZ.COMMON, LMDZ.GENERIC, ICOSA_LMDZ, ARCH)

=== Saturne (simulation sur 96 niveaux, Bardet et al. Nature Astronomy 2022) ===

* XIOS revision 1583 

* ICOSAGCM revision 765 

* IOIPSL revision 431 

* FCM_V1.2 revision 12 

* Physics revision 2305 (LMDZ.COMMON, LMDZ.GENERIC, ICOSA_LMDZ, ARCH)

=== Saturne (simulation avec la GWD paramétrisation sur 96 niveaux, chapitre 6 PhD Bardet) ===

* XIOS revision 1583 

* ICOSAGCM revision 765 

* IOIPSL revision 431 

* FCM_V1.2 revision 12 

* Physics revision 2403 (LMDZ.COMMON, LMDZ.GENERIC, ICOSA_LMDZ, ARCH)

=== Uranus & Neptune (old version) ===

* XIOS revision 1944
* ICOSAGCM revision 765
* IOIPSL revision 431
* FCM_V1.2 revision 12 

* Physics revision 2413 (LMDZ.COMMON, LMDZ.GENERIC, ICOSA_LMDZ, ARCH)

=== Uranus & Neptune (new version) ===

* XIOS revision 2203
* ICOSAGCM revision (20/08/2021)
* IOIPSL revision 450
* FCM_V1.2 revision 12 

* Physics revision 2555 (LMDZ.COMMON, LMDZ.GENERIC, ICOSA_LMDZ, ARCH)

works on 03/03/2022 for all HEAD

== Les fichiers d'architecture (Installation sur un nouveau cluster) ==

En pratique LMDZ.COMMON, ICOSA_LMDZ et IOIPSL peuvent utiliser exactement le même fichier arch.fcm ; mais celui pour ICOSAGCM est légèrement différent (les %FPP_DEF diffèrent, peut-être aussi le %FPP).

make_icosa_lmdz doit être lancé avec -arch_path ../ARCH puisque les arch.env et arch.path communs se trouvent dans ../ARCH (l'option -arch_path ne concerne d'ailleurs que les fichiers arch.env et arch.path; le fichier arch.fcm recherché sera toujours celui dans le "arch" de chacun des modèles. Donc il faut bien mettre pour chacun des quatre modèles (LMDZ.COMMON, ICOSA_LMDZ,IOIPSL et ICOSAGCM) le arch.fcm dans le sous-dossier arch/ du modèle correspondant.

XIOS est écrit en C++, et pas en Fortran. Le fichier arch.fcm correspondant est donc nécessairement différent de celui des quatre autres modèles. Pour créér ce arch.fcm, prendre exemple sur les fichiers .fcm déjà présent dans XIOS/arch/, avec une architecture similaire à celle du nouveau cluster. En particulier, il faut utiliser des versions de gcc/gfortran > 6+. Il faut absolument avoir une bibliothèque HDF5 compilée en parallèle, ainsi que netcdf-C et netcdf-fortran (et peut être aussi netcdf-Cxx) compilé en parallèle. Sans ça, il sera peut-être possible de compiler le modèle, mais pas de lancer une simulation en utilisant XIOS.

== Commencer une nouvelle simulation en utilisant le schéma des panaches thermiques ==

# Comme pour toutes nouvelles simulations, il faut commencer par un run 1D de plusieurs décennies permettant d'obtenir un profil de température (<code>temp_profile.txt</code>) et les coefficients ap et bp (<code>apbp.txt</code>) en équilibre radiatif-convectif pour la planète que l'on étudie. Ce run 1D s'effectue dans les dossiers jupiter1d, saturn1d, neptune1d et/ou uranus1d, avec ces options pour le callphys.def :

callrad = true calladj = true tous les autres mots clés à false (y compris <code>calltherm</code>).

<ol start="2" style="list-style-type: decimal;">
<li>Faire un "makestart" run : permet d'obtenir les fichiers restart à partir du profil de température initial. Il s'agit du run dont l'état initial est le profil de température créé par run 1D (<code>temp_profile.txt</code> est appliquer à chaque point de grille horizontale du modèle). L'état initial est ainsi une planète isotherme horizontalement mais qui varie verticalement. Pour ce run, aucun traceur ne va être utilisé dans le modèle. Néanmoins, il faut renseigner au modèle le nombre de traceurs que nous souhaitons utiliser avec le schéma des panaches thermiques afin qu'il puisse créer la dimension nq et le champ q dans les fichiers <code>restart_icosa.nc</code> et <code>restartfi.nc</code>. dans '''callphys.def''' :
<pre> traceur = true</pre>
dans '''run_icosa.def''' : <code>nqtot = 2</code> (par exemple, h2o_vap et h2o_ice) dans '''traceur.def''' :
<pre> 2
h2o_vap
h2o_ice</pre></li>
<li>Ajout des quantités pour chaque traceur. Ici, nous ajoutons les profils pour chacun des traceurs dans les fichiers <code>restart_icosa.nc</code> et <code>restartfi.nc</code> "à la main" en utilisant le programme python <code>/processing_codes/tracer_settings.py</code> pour obtenir des fichiers restart avec la bonne abondance d'eau dans le cas présent.</li>
<li>Commencer la simulation : Il ne reste plus qu'à lancer la simulation 3D avec les nouveaux fichiers restart et les réglages suivant : dans '''callphys.def''' :
<pre> traceur = true
calltherm = true
# thermal plume model options:
divmpl = true
r_aspect_thermals = 2.0
tau_thermals = 0.0
betalpha = 0.9
afact = 0.7
fact_epsilon = 2.e-4
alpha_max = 0.7
fomass_max = 0.5
pres_limit = 2.e5
water = true
watercond = true
waterain = true
evap_prec = true</pre></li></ol>

dans '''run_icosa.def''' : <code>nqtot = 2</code> dans '''traceur.def''' :

<pre> 2
h2o_vap
h2o_ice</pre>
Enfin, ajouter la déclaration et l'écriture des variables relatives à l'utilisation du schéma des thermiques (<code>h2o_vap</code>, <code>h2o_ice</code>, <code>w_plm</code>) dans les fichiers XML de la physique.

== Erreur en début de run, SEGMENTATION FAULT dans la routine <code>advect.f90</code> (ICOSAGCM/ppsrc/transport) au niveau de l'appel à <code>cross_product2</code>. ==

La routine <code>cross_product2</code> effectue des produits vectoriels. Cette erreur vient du fait qu'en essayant de faire des optimisations de calcul, le compilateur se prend "les pieds dans le tapis" et ne sait plus faire le produit vectoriel. Pour empêcher le compilateur de faire des optimisations trop poussées, il faut ajouter l'option <code>-fp-model strict</code> au niveau du mot-clé <code>%BASE_FFLAGS</code> dans tous les fichiers <code>.fcm</code> du modèle (LMDZ.COMMON / ICOSAGCM / ICOSA_LMDZ)

== Branche master version du 04/04/2022 dans le dossier Jupiter ==

La version actuelle des fichiers de réglage dans le dossier Jupiter sont réglés pour utiliser des traceurs. Si vous souhaiter utiliser le modèle sans traceur, il vous faut modifier les fichiers suivants :

run_icosa.def

<pre> nqtot = 0</pre>
traceur.def

<pre> 0</pre>
context_dynamico.xml (ligne 125)

<pre> <field id="q_start" name="q" grid_ref="grid_q_start" prec="8"/> </pre>
== PROFILING with DYNAMICO-Giant: ==

<code>-pg</code>: Generate extra code to write profile information suitable for analysis program <code>gprof</code>. You must use this option when compiling the sources files you want data about and you must '''also use it when linking'''. 1. For DYNAMICO-Giant: in the arch.fcm file of each part of the model (ICOSAGCM, ICOSA_LMDZ, LMDZ.COMMON and XIOS), you have to add the option <code>-pg</code> to <code>%PROD_FFLAGS</code> and <code>%BASE_LD</code> 2. Compile as usual 3. Execute your code as usual 4. Run <code>gprof</code> tool:

<pre> gprof icosa_lmdz.exe goon.out > analysis.txt</pre>

[[Category:Generic-DYNAMICO]]

Advanced Use of the GCM

2024-02-27T09:24:03Z

LTeinturier: /* The generic condensable scheme */

== Running in parallel ==

For large simulation (long run, high resolution etc...), the computational cost can be huge and hence the run time very long.
To overcome this issue, the model can be run in parallel. This however requires a few extra steps (compared to compiling and running the serial version of the code).
For all the details see [[Parallelism | the dedicated page]].

== Disambiguation between ifort, mpif90, etc. ==

For users not used to compilers and/or compiling and running codes in parallel, namely in MPI mode, there is often some confusion which hopefully the following paragraph might help clarify:
* the compiler (typically gfortran, ifort, pgfortran, etc.) is the required tool to compile the Fortran source code and generate an executable. It is strongly recommended that libraries used by a program are also compiled using the same compiler. Thus if you plan to use different compilers to compile the model, note that you should also have at hand versions of the libraries it uses also compiled with these compilers.
* the MPI (Message Passing Interface) library is a library used to solve problems using multiple processes by enabling message-passing between the otherwise independent processes. There are a number of available MPI libraries out there, e.g. OpenMPI, MPICH or IntelMPI to name a few (you can check out the [[Building an MPI library]] page for some information about installing an MPI library). The important point here is that on a given machine the MPI library is related to a given compiler and that it provides related wrappers to compile and run with. Typically (but not always) the compiler wrapper is '''mpif90''' and the execution wrapper is '''mpirun'''. If you want to know which compiler is wrapped in the '''mpif90''' compiler wrapper, check out the output of:
<pre>
mpif90 --version
</pre>
* In addition a second type of parallelism, shared memory parallelism known as OpenMP, is also implemented in the code. In contradistinction to MPI, OpenMP does not require an external library but is instead implemented as a compiler feature. At run time one must then specify some dedicated environment variables (such as OMP_NUM_THREADS and OMP_STACKSIZE) to specify the number of threads to use per process.
* In practice one should favor compiling and running with both MPI and OPenMP enabled.
* For much more detailed information about compiling and running in parallel, check out the [[Parallelism | the page dedicated to Parallelism]].

== A word about the IOIPSL and XIOS libraries ==
* The IOIPSL (Input Output IPSL) library is a library that has developed by the IPSL community to handle input and outputs of (mostly terrestrial) climate models. For the Generic PCM only a small part of this library is actually used, related to reading and processing the input [[The_run.def_Input_File | run.def]] file. For more details check out the [[The IOIPSL Library]] page.
* The [https://forge.ipsl.jussieu.fr/ioserver/wiki XIOS] (Xml I/O Server) library is based on client-server principles where the server manages the outputs asynchronously from the client (the climate model) so that the bottleneck of writing data in a parallel environment is alleviated. All aspects of the outputs (name, units, file, post-processing operations, etc.) are then controlled by dedicated XML files which are read at run-time. Using XIOS is currently optional (and requires compiling the GCM with the XIOS library). More about the XIOS library, how to install and use it, etc. [[The XIOS Library| here]].

== Playing with the output files ==

=== Changing the output temporal resolution and time duration ===

* To change the total time of a simulation, you need to open the 'For all the details see [[The_run.def_Input_File | run.def]]. file and change the variable 'nday':

<pre>
nday = 1000 # this means the simulation will run for 1000 days ; and that the associated output files will also be computed for a total duration of 1000 days
</pre>

Note: in the example, they are not necessarily 1000 Earth days, because it depends on the definition of the day duration that has been taken in the start files.

* To change the temporal resolution of the output files, you need to open the [[The_run.def_Input_File | run.def]] file and change the variable 'ecritphy':

<pre>
ecritphy = 200 # this means the simulation will write variables in the output files every 200 time steps of the simulation.
</pre>

Note: The output temporal resolution of the output files then depends also on the number of timestep per day ('day_step' variable in [[The_run.def_Input_File | run.def]] file). In this example:

<pre>
nday = 1000
daystep = 800
ecritphy = 200
</pre>

The output file will provide results every 0.25 days (800/200), and for a total duration of 1000 days (so 4000 time values in total).

=== Changing the output variable ===

To select the variable provided in the output file diagfi.nc, you simply need to add the list of variables needed in the [[The_diagfi.def_Input_File | diagfi.def]].

Note for experts: Some technical variables need to be de-commented in 'physiq_mod.F90' file to be written in the output files.

=== Spectral outputs ===

It is possible to provide spectral outputs such as the OLR (Outgoing Longwave Radiation, i.e. the thermal emission of the planet at the top of the atmosphere), the OSR (Outgoing Stellar Radiation, i.e. the light reflected by the planet at the top of the atmosphere), or the GSR (Ground Stellar Radiation, i.e. the light emitted by the star that reaches the surface of the planet).

For this, you need to activate the option 'specOLR' in the [[The_callphys.def_Input_File | callphys.def]] file, as follows:

<pre>
specOLR = .true.
</pre>

The simulations will then create diagspec_VI.nc and diagspec_IR.nc files (along with the standard diagfi.nc file), which contain the spectra of OLR, OSR, GSR, etc.

Note: The resolution of the spectra is defined by that of the correlated-k (opacity) files used for the simulation.

=== Statistical outputs ===

TBD (explain how to compute stats.nc files as well as what is inside)

== How to Change Vertical and Horizontal Resolutions ==

=== When you are using the regular longitude/latitude horizontal grid ===
To run at a different grid resolution than available initial conditions files, one needs to use the tools ''newstart.e'' and ''start2archive.e''

For example, to create initial states at grid resolution 32×24×25 from NetCDF files start and startfi at grid resolution 64×48×32 :

* Create file ''start_archive.nc'' with ''start2archive.e'' compiled at grid resolution 64×48×32 using old file ''z2sig.def'' used previously
* Create files ''restart.nc'' and ''restartfi.nc'' with ''newstart.e'' compiled at grid resolution 32×24×25, using a new file ''z2sig.def'' (more details below on the choice of the ''z2sig.def'').
* While executing ''newstart.e'', you need to choose the answer '0 - from a file start_archive' and then press enter to all other requests.

==== What you need to ''know'' about the ''z2sig.def'' file ====

For a model with Nlay layers, the [[The_z2sig.def_Input_File | z2sig.def]] file must contain at least Nlay+1 lines (the other not being read).

The first line is a scale height ($$H$$). The following lines are the target pseudo-altitudes for the model from the bottom up ($$z_i$$).
The units do not matter as long as you use the same ones for both.

The model will use these altitudes to compute a target pressure grid ($$p_i$$ ) as follows:
\begin{align}
\label{def:pseudoalt}
p_i &= p_s \exp(-z_i/H),
\end{align}
where $$p_s$$ is the surface pressure.

As you can see, the scale height and pseudo altitudes enter the equation only through their ratio. So they do not have to to be the real scale-height and altitudes of the atmosphere you are simulating.
So you can use the same [[The_z2sig.def_Input_File | z2sig.def]].def for different planets.

There is no hard rule to follow to determine the altitude/pressure levels you should use. As a rule of thumb, layers should be thiner near the surface to properly resolve the surface boundary layer. Then they should gradually increase in size over a couple scale heights and transition to constant thickness above that. Of course, some specific applications may require thinner layers in some specific parts of the atmospheres.

A little trick for those who prefer to think in terms of (log)pressure: if you use $$H= 1/\ln 10 \approx 0.43429448$$, then $$z_i=x$$ corresponds to a pressure difference with the surface of exactly x pressure decades (i.e. at $$z=1$$, $$p=0.1p_s$$). This is particularly useful for giant-planet applications.



'''WE SHOULD WRITE A PAGE (LINK HERE) ABOUT HYBRID VERTICAL COORDINATES'''

=== When you are using the DYNAMICO icosahedral horizontal grid ===

The horizontal resolution for the DYNAMICO dynamical core is managed from several setting files, online during the execution.
To this purpose, each part of the GCM managing the in/output fields ('''ICOSAGCM''', '''ICOSA_LMDZ''', '''XIOS''') requires to know the input and output grids:

'''1. ''context_lmdz_physics.xml'':'''

You can find several grid setup already defined:

<syntaxhighlight lang="xml" line>
<domain_definition>
<domain id="dom_96_95" ni_glo="96" nj_glo="95" type="rectilinear" >
<generate_rectilinear_domain/>
<interpolate_domain order="1"/>
</domain>

<domain id="dom_144_142" ni_glo="144" nj_glo="142" type="rectilinear" >
<generate_rectilinear_domain/>
<interpolate_domain order="1"/>
</domain>

<domain id="dom_512_360" ni_glo="512" nj_glo="360" type="rectilinear" >
<generate_rectilinear_domain/>
<interpolate_domain order="1"/>
</domain>

<domain id="dom_720_360" ni_glo="720" nj_glo="360" type="rectilinear">
<generate_rectilinear_domain/>
<interpolate_domain order="1"/>
</domain>

<domain id="dom_out" domain_ref="dom_720_360"/>
</domain_definition>
</syntaxhighlight>

In this example, the output grid for the physics fields is defined by

<syntaxhighlight lang="xml">
<domain id="dom_out" domain_ref="dom_720_360"/>
</syntaxhighlight>

which is an half-degree horizontal resolution. To change this resolution, you have to change name of the '''domain_ref''' grid, for instance:

<syntaxhighlight lang="xml">
<domain id="dom_out" domain_ref="dom_96_95"/>
</syntaxhighlight>

'''2. ''run_icosa.def'': setting file to execute a simulation'''

In this file, regarding of the horizontal resolution intended, you have to set the number of subdivision on the main triangle.
For reminder, each hexagonal mesh is divided in several main triangles and each main triangles are divided in suitable number of sub-triangles according the horizontal resolution

<syntaxhighlight lang="bash" line>
#nbp --> number of subdivision on a main triangle: integer (default=40)
# nbp = sqrt((nbr_lat x nbr_lon)/10)
# nbp: 20 40 80 160
# T-edge length (km): 500 250 120 60
# Example: nbp(128x96) = 35 -> 40
# nbp(256x192)= 70 -> 80
# nbp(360x720)= 160 -> 160
nbp = 160
</syntaxhighlight>

If you have chosen the 96_95 output grid in ''context_lmdz_physics.xml'', you have to calculate $$nbp = \sqrt(96x95) / 10 = 10$$ and in this case

<syntaxhighlight lang="bash">
nbp = 20
</syntaxhighlight>

After the number of subdivision of the main triangle, you have to define the number subdivision over each direction. At this stage you need to be careful as the number of subdivisions on each direction:
* needs to be set according to the number of subdivisions on the main triangle '''nbp'''
* will determine the number of processors on which the GCM will be most effective

<syntaxhighlight lang="bash" line>
## sub splitting of main rhombus : integer (default=1)
#nsplit_i=1
#nsplit_j=1
#omp_level_size=1
###############################################################
## There must be less MPIxOpenMP processes than the 10 x nsplit_i x nsplit_j tiles
## typically for pure MPI runs, let nproc = 10 x nsplit_i x nsplit_j
## it is better to have nbp/nsplit_i > 10 and nbp/nplit_j > 10
###############################################################
#### 40 noeuds de 24 processeurs = 960 procs
nsplit_i=12
nsplit_j=8

#### 50 noeuds de 24 processeurs = 1200 procs
#nsplit_i=10
#nsplit_j=12
</syntaxhighlight>

With the same example as above, the 96_95 output grid requires:
$$nsplit_i < 2$$ and $$nsplit_j < 2$$
We advise you to select:

<syntaxhighlight lang="bash">
## sub splitting of main rhombus : integer (default=1)
nsplit_i=1
nsplit_j=1
</syntaxhighlight>

and using 10 processors.

== How to Change the Topography (or remove it) ==

The generic model can use in principle any type of surface topography, provided that the topographic data file is available in the right format, and put in the right place. The information content on the surface topography is contained in the ''startfi.nc'', and we do have developed tools (see below) to modify the ''startfi.nc'' to account for a new surface topography.

To change the surface topography of a simulation, we recommend to follow the procedure detailed below:

* Create file ''start_archive.nc'' with ''start2archive.e'' compiled at the same (horizontal and vertical) resolution than the ''start.nc'' and ''startfi.nc'' files.
* Create files ''restart.nc'' and ''restartfi.nc'' with ''newstart.e'' compiled again at the same (horizontal and vertical) resolution.
* While executing ''newstart.e'', you need to choose the answer '0 - from a file start_archive' and then press enter to all other requests.
* At some point, the script ''newstart.e'' asks you to chose the surface topography you want from the list of files available in your 'datagcm/surface_data/' directory.

We do have a repository of for Venus, Earth and Mars through time available at https://web.lmd.jussieu.fr/~lmdz/planets/LMDZ.GENERIC/datagcm/surface_data/. You can download the surface topography files and place them in your 'datagcm/surface_data/' directory.

We also offer a tutorial to design new topography maps here: https://lmdz-forge.lmd.jussieu.fr/mediawiki/Planets/index.php/Building_Surface_Topography_Files

----

Special note: To remove the topography, you can simply add the following tag in callphys.def:

<pre>
nosurf = .true.
</pre>

== How to Change the Stellar Spectrum ==

To simulate the effect of the star's radiation on a given planetary atmosphere, it is necessary to accurately represent the stellar spectrum (spectral shape and total bolometric flux) at the top of this atmosphere. In the model, we have set up two different options to model the stellar spectra of any star.

=== Black Body Stellar Spectra ===

First, it is possible to simply use a black body. In this case, the stellar spectrum depends only on the effective temperature of the star which is provided to the model.

For this, you need to activate the option 'stelbbody' in the [[The_callphys.def_Input_File | callphys.def]] file, as follows:

<pre>
stelbbody = .true.
</pre>

and then add, also in the [[The_callphys.def_Input_File | callphys.def]] file, the following line:

<pre>
stelTbb = 3500.0
</pre>

to specify the effective temperature of the host star. (in this example, we have chosen a M-star with an effective temperature of 3500K)

=== Pre-Tabulated spectra ===

Second, the model can read a file containing any pre-computed stellar spectrum. Traditionally, we have used synthetic spectra from the PHOENIX database, that we adapt to the Generic PCM by decreasing the spectral resolution (use 10000 points with a fixed spectral resolution of 0.001 micron) and by adapting the units (in W/m2/micron). This is the option that is generally preferred to better represent the effect of the star (whose real spectrum can strongly deviate from the black body approximation).

For this, you need to make sure the option 'stelbbody' in the [[The_callphys.def_Input_File | callphys.def]] file is equal to .false. (it not specified, by default stelbbody is assumed to be .false.).

Then you need to add in the [[The_callphys.def_Input_File | callphys.def]] file, the following line:

<pre>
startype = 1
</pre>

and change the value of startype depending on the star you want to model. (here 1 means we use the solar spectrum)

To know which stellar spectra are available, you need to open the file LMDZ.GENERIC/libf/phystd/ave_stelspec.F90 and adapt the value of startype accordingly. You also need to make sure the spectra are available in your /datadir/stellar_spectra (or /datagcm/stellar_spectra) directory.

To calculate the true stellar spectrum at the top of the atmosphere, the Generic PCM renormalizes the stellar spectrum by the bolometric flux at 1 Astronomical Unit (AU) provided by the user, which it then converts into the true stellar spectrum by using the star-planet distance.

To specify the flux at 1 AU, you need to add in the [[The_callphys.def_Input_File | callphys.def]] file, the following line:

<pre>
Fat1AU = 1366.0
</pre>

(here 1366W/m2 corresponds to the flux at 1AU for the Sun)

1st NOTE: We will improve this second part (Martin and Mathilde) by the end of 2022.

2nd NOTE: The Generic PCM eventually has the capability to run without any stellar flux. To do that, you can simply put Fat1AU = 0. (@LUCAS_TEINTURIER, could you check that?)

== How to Change the Opacity Tables ==

'''SECTION ET PAGES EN COURS D'ÉCRITURE'''

The model uses opacity tables to compute heating rates throughout the atmosphere. These opacity tables are generated "offline" for a given set of pressures, temperatures, a given composition and a specific spectral decomposition.

=== Getting opacity tables for your desired atmospheric composition ===

* You should first check our common repository here (https://web.lmd.jussieu.fr/~lmdz/planets/LMDZ.GENERIC/datagcm/corrk_data/) to check whether your desired opacity table is not already included.
We have a dedicated page on how to build new opacity tables here: https://lmdz-forge.lmd.jussieu.fr/mediawiki/Planets/index.php/Building_Opacity_Tables ; N

=== Implementing your opacity tables in the Generic PCM ===

, follow these steps:

* copy your directory containing your opacity tables in .... It has to contain... [[Building_Opacity_Tables | you need to build opacity tables]]
* change corrkdir = ... in [[The_callphys.def_Input_File | callphys.def]] with the name of that directory.
* change gases.def: it has to be consistent with Q.dat. A special case concerns...
* change -b option when compiling the model with makelmdz_fcm: it has to correspond to the number of bands (in the IR x in the visible) of the new opacity tables. For instance, compile with -b 38x26 if you used 38 bands in the infrared and 26 in the visible to generate the opacity tables.

== How to Change the Aerosols Optical Properties ==

Aerosol optical properties are represented using three distinct properties: the extinction coefficient (Q_ext), the single scattering albedo (omega) and the asymmetry factor (g).

The Generic PCM can compute the radiative effects of any aerosol, provided that they optical properties (Q_ext, omega, g) are tabulated and provided in the right format.

=== Getting optical properties for your aerosols ===

* You should first check our common repository here (https://web.lmd.jussieu.fr/~lmdz/planets/LMDZ.GENERIC/datagcm/aerosol_properties/) to check whether your favorite aerosol is not already included. The optical properties of each aerosol is built using two distinct files: one in the 'visible' (which is used in the visible part of the radiative transfer, to compute the fate of stellar radiation) and one in the 'infrared' (which is used in the thermal infrared part of the radiative transfer, to compute the fate of thermal emission by the surface and atmosphere).

For instance, if you want to include the radiative effect of CO2 ice clouds, then you just need the files:

- optprop_co2ice_vis_n50.dat

- optprop_co2ice_ir_n50.dat

* Otherwise, you can create your own tables of optical properties, using existing databases. Check this page to learn how to do this: https://lmdz-forge.lmd.jussieu.fr/mediawiki/Planets/index.php/Building_Tables_Of_Aerosol_Optical_Properties

=== Implementing your aerosols in the Generic PCM ===

There are two available aerosol schemes you can use to add a new aerosol in the Generic PCM:

====The n-layer aerosol scheme====

You can use the n-layer scheme (implemented by Jan Vatant) to easily prescribe an aerosol vertical distribution. The scheme is activated by adding aeronlay=.true. in callphys.def.

In this scheme, each layer can have different optical properties, particle sizes, etc. You can use different options (e.g. to fix the aerosol distribution between two atmospheric pressures) of the scheme with aeronlay_choice = 1, 2, etc.

You will need to provide the name of your aerosols optical properties tables in callphys.def. For instance:
* optprop_aeronlay_vis = 'optprop_saturn_vis_n20.dat'
* optprop_aeronlay_ir = 'optprop_saturn_ir_n20.dat'

(here, to use the optical properties of aerosols for Saturn)

We encourage you to search for the keyword "aeronlay" in the source code to use more specific options of the scheme.

====The generic condensable scheme====

You can use the generic condensable scheme (implemented by Lucas Teinturier) to easily compute the radiative effect of cloud particles formed by condensation.
The scheme has to be used conjointly to the generic condensation scheme.
To activate it, one needs to add the '''aerogeneric=n'''( with n>0 the number of condensable species handled in the scheme) in callphys.def. On top of that, one needs to add in traceur.def, on the solid/ice traceur the option '''is_rgcs = 1'''.

In the model itself, one needs then to manually add your aerosol by modifying the suaer_corrk.F90 routine (if your favorite aerosol is not there yet already) to specify the correct names for your condensing species. A more dynamical/flexible approach will be added at some point, so one won't need to directly modify the code.

More information on how to use the scheme is provided here: https://lmdz-forge.lmd.jussieu.fr/mediawiki/Planets/index.php/Radiative_Generic_Condensable_Specie

====Last Comments====

Don’t forget to add your new aerosol species in traceur.def and adapt the -s option (with the correct number of radiatively active aerosols) at the compilation stage.

== How to Manage Tracers ==

Tracers are managed thanks to the [https://lmdz-forge.lmd.jussieu.fr/mediawiki/Planets/index.php/The_traceur.def_Input_File ''traceur.def''] file.

Specific treatment of some tracers (e.g., water vapor cycle) can be added directly in the model and an option added in [https://lmdz-forge.lmd.jussieu.fr/mediawiki/Planets/index.php/The_callphys.def_Input_File ''callphys.def''] file.

== Use the Z of LMDZ : Zoomed version ==

Do we need this? Has anyone already made use of the zoom module?

[[Category:Generic-Model]]

Tool Box

2023-12-07T10:52:26Z

LTeinturier: /* XARRAY python library (more modern) */

== Pre-processing Tools ==
=== newstart: a fortran program to modify start files ===

Newstart is an interactive tool to modify the start files (''start.nc'' and ''startfi.nc'').

To be usable, ''newstart'' should be compile in the ''LMDZ.COMMON'' directory by using the following command line:
<syntaxhighlight lang="bash">
./makelmdz_fcm -arch my_arch_file -p std -d 64x48x30 newstart
</syntaxhighlight>
In the example, my_arch_file is the name the arch files (see [https://lmdz-forge.lmd.jussieu.fr/mediawiki/Planets/index.php/The_Target_Architecture_(%22arch%22)_Files arch] ) and 64x48x30 is the resolution of the physical grid.
Then copy the executable from the ''LMDZ.COMMON/bin'' directory to your bench directory.

When you execute newstart, you can use both a ''start2archive'' file or the start files (''start.nc'' and ''startfi.nc''). Then the interactive interface will propose to modify several physical quantities such as the gravity, the surface pressure or the rotation of the planet. At the end of the procedure, two files are created: '' '''re'''start.nc'' and '' '''re'''startfi.nc''. They can be renamed and used as start files to initialize a new simulation.

We have prepared a simple tutorial to learn how to modify ''start.nc'' and ''startfi.nc'' files (= the the initial conditions) for the Generic PCM: https://lmdz-forge.lmd.jussieu.fr/mediawiki/Planets/index.php/Modify_start_Files

=== start2archive ===

The start2archive tool is similar to newstart in the sense that it can be used to modify the start files. But start2archive can modify the resolution of the physical grid, the topography and the surface thermal inertia while newstart cannot. It is also useful to create an archive of different starting states, then extractable as start files.
The command line to compile start2archive is similar to the one used for newstart:
<syntaxhighlight lang="bash">
./makelmdz_fcm -arch my_arch_file -p std -d 64x48x30 start2archive
</syntaxhighlight>
To modify the resolution, you should first create a start_archive (by using start2archive) file at the used resolution, then compile a newstart file at the new resolution. Newstart will interpolate all the physical quantities on the new grid.

=== other third party scripts and tools ===

You can easily modify start.nc and startfi.nc netcdf files with the xarray Python library. Below you can find an easy example where we modify the surface temperature field.

<syntaxhighlight lang="python">
from numpy import *
import numpy as np
import matplotlib.pyplot as mpl
import math
import xarray as xr

#FIRST WE GET THE DATA FROM THE GCM SIMULATION
nc = xr.open_dataset('startfi.nc',decode_times=False) # can be any netcdf file (e.g. start/startfi.nc files)

# BELOW PHYSICAL VARIABLES
physical_points=nc['physical_points']
lat=nc['latitude']
lon=nc['longitude']
aire_GCM=nc['area']

# BELOW THE VARIABLE WE WANT TO UPDATE
tsurf=nc['tsurf']
new_tsurf = np.empty(len(physical_points))

# LOOP TO MODIFY THE VARIABLE
for i in range(0,len(physical_points),1):
new_tsurf[i]=300. # here you put whatever you want ; in this exemple, we assume an isothermal temperature distribution

nc['tsurf'].values = new_tsurf
nc.to_netcdf('restartfi.nc')

# SANITY CHECK PLOTS

fig = mpl.figure(1)
mpl.plot(physical_points,tsurf)
mpl.plot(physical_points,new_tsurf)
mpl.xlabel('GCM Physical Points')
mpl.ylabel('Tsurf (K)')
mpl.show()
</syntaxhighlight>

== Post-processing tools ==
=== zrecast ===

With this program you can recast atmospheric (i.e.: 4D-dimentional longitude-latitude-altitude-time) data from
GCM outputs (e.g. as given in diagfi.nc files) onto either ''pressure'' or ''altitude'' above ''areoid vertical'' coordinates.
Since integrating the hydrostatic equation is required to recast the data, the input file must contain surface pressure
and atmospheric temperature, as well as the ground geopotential.
If recasting data onto ''pressure'' coordinates, then the output file name is given by the input file name to which ''_P.nc'' will be appened. If recasting data onto altitude above areoid coordinates, then a ''_A.nc'' will be appened.

=== mass stream function ===

The mass stream function (and the total angular momentum) can be computed from a diagfi.nc or a stats.nc, using the '''streamfunction.F90''' script. The script is located at

<syntaxhighlight lang="bash">
trunk/LMDZ.GENERIC/utilities
</syntaxhighlight>

To compile the script, open the ''compile'' file in the same directory and do the following:
* Replace "pgf90" with your favorite fortran compiler
* replace "/distrib/local/netcdf/pgi_7.1-6_32/lib" with the lib address and directory that contains your NetCDF library (file ''libnetcdf.a'').
* Replace "/distrib/local/netcdf/pgi_7.1-6_32/include" with the address of the directory that contains the NetCDF include file (''netcdf.inc'').
* You can mess with the compiling options but it is not mandatory.

Once the script is compiled, copy it in the same directory as your '''.nc''' file and run

<syntaxhighlight lang="bash">
./streamfunction.e
</syntaxhighlight>
The script will ask you for the name of your '''.nc''' file, and will run and produce a new '''nameofyourfile_stream.nc''' file.

'''Be careful''' : In this new file, all fields are temporally and zonally averaged.

If you want to use '''python''' instead of '''fortran''', you can take a look at this [https://github.com/aymeric-spiga/dynanalysis repo]. It hosts a tool to perform dynamical analysis of GCM simulations (and therefore, it computes the mass stream function and a lot of other stuff), but it is tailored for Dynamico only. This repo also takes care of recasting (it does the job of both ''zrecast.F90'' and ''streamfunction.F90'')

== Continuing Simulations ==

=== manually ===

At the end of a simulation, the model generates restart files (files 'restart.nc' and 'restartfi.nc') which contain the final state of the model. The 'restart.nc' and 'restartfi.nc' files have the same format as the 'start.nc' and 'startfi.nc' files, respectively.

These files can in fact be used as initial states to continue the simulation, using the following renaming command lines:

<syntaxhighlight lang="bash">
mv restart.nc start.nc
mv restartfi.nc startfi.nc
</syntaxhighlight>

Running a simulation with these start files will in fact resume the simulation from where the previous run ended.

=== with bash scripts ===

We have set up very simple bash scripts to automatize the launching of chain simulations. Here is an example of bash script that does the job:

<syntaxhighlight lang="bash">
#!/bin/bash
###########################################################################
# Script to perform several chained LMD Mars GCM simulations
# SET HERE the maximum total number of simulations

nummax=100

###########################################################################

echo "---------------------------------------------------------"
echo "STARTING LOOP RUN"
echo "---------------------------------------------------------"

dir=`pwd`
machine=`hostname`
address=`whoami`

# Look for file "num_run" which should contain
# the value of the previously computed season
# (defaults to 0 if file "num_run" does not exist)
if [[ -r num_run ]] ; then
echo "found file num_run"
numold=`cat num_run`
else
numold=0
fi
echo "numold is set to" ${numold}

# Set value of current season
(( numnew = ${numold} + 1 ))
echo "numnew is set to" ${numnew}

# Look for initialization data files (exit if none found)
if [[ ( -r start${numold}.nc && -r startfi${numold}.nc ) ]] ; then
\cp -f start${numold}.nc start.nc
\cp -f startfi${numold}.nc startfi.nc
else
if (( ${numold} == 99999 )) ; then
echo "No run because previous run crashed ! (99999 in num_run)"
exit
else
echo "Where is file start"${numold}".nc??"
exit
fi
fi

# Run GCM -- THIS LINE NEEDS TO BE MODIFIED WITH THE CORRECT GCM EXECUTION COMMAND
mpirun -np 8 gcm_64x48x26_phystd_para.e < diagfi.def > lrun${numnew}

# Check if run ended normaly and copy datafiles
if [[ ( -r restartfi.nc && -r restart.nc ) ]] ; then
echo "Run seems to have ended normaly"

\mv -f restart.nc start${numnew}.nc
\mv -f restartfi.nc startfi${numnew}.nc

else
if [[ -r num_run ]] ; then
\mv -f num_run num_run.crash
else
echo "No file num_run to build num_run.crash from !!"
# Impose a default value of 0 for num_run
echo 0 > num_run.crash
fi
echo 99999 > num_run
exit
fi

# Copy other datafiles that may have been generated
if [[ -r diagfi.nc ]] ; then
\mv -f diagfi.nc diagfi${numnew}.nc
fi
if [[ -r diagsoil.nc ]] ; then
\mv -f diagsoil.nc diagsoil${numnew}.nc
fi
if [[ -r stats.nc ]] ; then
\mv -f stats.nc stats${numnew}.nc
fi
if [[ -f profiles.dat ]] ; then
\mv -f profiles.dat profiles${numnew}.dat
\mv -f profiles.hdr profiles${numnew}.hdr
fi

# Prepare things for upcoming runs by writing
# value of computed season in file num_run
echo ${numnew} > num_run

# If we are over nummax : stop
if (( $numnew + 1 > $nummax )) ; then
exit
# If not : restart the loop (copy the executable and run the copy)
else
\cp -f run_gnome exe_mars
./exe_mars
fi

</syntaxhighlight>

'''Summary of what this bash script does''':

* It reads the file 'num_run' which contains the step of the simulation.
If num_run is
<syntaxhighlight lang="bash">
5
</syntaxhighlight>
then the script expects to read start5.nc and startfi5.nc.
* It modifies start5.nc and startfi5.nc into start.nc and startfi.nc, respectively.
* It runs the GCM.
* It modifies restart.nc and restartfi.nc into start6.nc and startfi6.nc
* It rewrite num_run as follows:
<syntaxhighlight lang="bash">
6
</syntaxhighlight>
* It restarts the loop until num_run reaches the value (defined in nummax):
<syntaxhighlight lang="bash">
100
</syntaxhighlight>

== Processing Output Files with NCOs ==

NCOs (NetCdf OperatorS) are a set of powerful command-line utilities – available on Linux, Mac and PC – that allow to perform useful (and very fast!) post-processing operations on netCDF GCM output files. Full documentation can be found on http://research.jisao.washington.edu/data_sets/nco/, but we provide below a few examples of command lines.

* How to calculate a time mean of a netCDF 'diagfi.nc' file

<syntaxhighlight lang="bash">
ncra -F -d Time,1,,1 diagfi.nc diagfi_MEAN.nc # format is "-d dimension,minimum,maximum,stride"
</syntaxhighlight>

* Subsetting time in a netCDF 'diagfi.nc' file

<syntaxhighlight lang="bash">
ncea -F -d Time,first,last diagfi.nc diagfi_subset.nc # format is "-d dimension,minimum,maximum" ; we recall you can type "ncdump -v time diagfi.nc" to see the Time values in the netCDF file.
</syntaxhighlight>

* Decimating a netCDF 'diagfi.nc' file in time

<syntaxhighlight lang="bash">
ncks -F -d Time,1,,8 diagfi.nc diagfi_decimated.nc # format is "-d dimension,minimum,maximum,stride" ; In this example, this means that data is extracted 1 time every 8 time steps, starting from the first time step (number 1), ending at the last time step).
</syntaxhighlight>

* Extract a variable from a netCDF 'diagfi.nc' file

<syntaxhighlight lang="bash">
ncks -v tsurf,temp,p diagfi.nc diagfi_out.nc # Here we created a new file named 'diagfi_out.nc' in which we only kept variables named 'tsurf' (surface temperatures), 'temp' (atmospheric temperatures) and p (atmospheric pressures).
</syntaxhighlight>

Again, more examples can be found on http://research.jisao.washington.edu/data_sets/nco/ .

== Data Handling and Visualization Software ==

There are several data handling and visualization tools that can be used to analyse and plot the results from GCM simulations (using the diagfi.nc NetCDF files). We provide below a panorama of most widely used solutions.

=== panoply ===
Panoply is a user-friendly tool for viewing raw NetCDF data, available here: https://www.giss.nasa.gov/tools/panoply/ . It is very convenient to make pretty visuals (see an example for the exoplanet TRAPPIST-1e). There are many options that can be used (map projections, masks, colorbars, shadows, etc.) to make your plots fancy. However, the tool is not very well suited for manipulating data (compute averages, statistics, etc.).

* Installation on Linux:
You simply need to download and untar the Package from the Panoply website. Note that to work it requires that Java and related Java Runtime environment (JRE) be installed on your system (otherwise it will simply look as if "nothing is happening" when you try to launch Panoply via the "panoply.sh" script), which on Ubuntu simply requires something like:
<syntaxhighlight lang="bash">
sudo apt install java
sudo apt install default-jre
</syntaxhighlight>

* Run on Linux (assuming the panoply.sh script is in a directory included in your PATH environment variable):
<syntaxhighlight lang="bash">
panoply.sh
</syntaxhighlight>

[[File:Example panoply.png|thumb|Screenshot of panoply showing here Generic PCM results for the exoplanet TRAPPIST-1e (surface temperatures)]]

=== ncview ===
ncview is another useful user-friendly tool for viewing raw NetCDF data. This is kind of a very archaic version of panoply, but it is convenient because it allows to have a very quick first look at netCDF data files.

Command line tool to visualize NetCDF data:
* Installation on Linux-Ubuntu:
<syntaxhighlight lang="bash">
sudo apt install ncview
</syntaxhighlight>
* Run on Linux:
<syntaxhighlight lang="bash">
ncview diagfi.nc
</syntaxhighlight>

[[File:Example ncview.png|thumb|Screenshot of ncview showing here Generic PCM results for the exoplanet Proxima b (OLR - Thermal Emission)]]

=== python scripts ===

Python scripts provide a very useful mean to analyse and visualize netCDF files.

==== NETCDF4 python library (old school) ====

You can use the netCDF4 python library to open a netCDF file and put data in tables that can then be manipulated and plotted.

Here is an exemple of how to open and read a netCDF file with Python:

<syntaxhighlight lang="python" line>
import numpy
from netCDF4 import Dataset

# HERE WE OPEN THE NETCDF FILE
nc = Dataset('diagfi.nc')

# HERE WE READ THE VARIABLES (1D OUTPUT)
Time=nc.variables['Time'][:]
lat=nc.variables['latitude'][:]
lon=nc.variables['longitude'][:]
al=nc.variables['altitude'][:]

# HERE WE READ THE AREA (2D OUTPUT)
aire_GCM=nc.variables['aire'][:][:]

# HERE WE READ 3D OUTPUTS
tsurf=nc.variables['tsurf'][:][:][:] # this is the surface temperature 3D field (time, latitude, longitude, altitude)

# HERE WE READ 4D OUTPUTS
temp=nc.variables['temp'][:][:][:][:] # this is the atmospheric temperature 4D field (time, latitude, longitude, altitude)

</syntaxhighlight>

And here is an exemple of how to manipulate the netCDF data (here to compute the time averaged surface temperatures):

<syntaxhighlight lang="python" line>
from numpy import *
import numpy as np

mean_tsurf=np.zeros((len(lat),len(lon)),dtype='f')

for i in range(0,len(Time)):
for j in range(0,len(lat)):
for k in range(0,len(lon)):
mean_tsurf[j,k]=mean_tsurf[j,k]+tsurf[i,j,k]*(1./len(Time))

</syntaxhighlight>

And here is a last exemple of how to plot the data (using matplotlib):

<syntaxhighlight lang="python" line>
import matplotlib.pyplot as plt

plt.figure(1)
plt.contourf(lon_GCM,lat_GCM,mean_tsurf)
plt.colorbar(label='Surface Temperature (K)')
plt.xlabel('Longitude ($^{\circ}$)')
plt.ylabel('Latitude ($^{\circ}$)')
plt.show()

</syntaxhighlight>

==== XARRAY python library (more modern) ====

Another useful library to deal with netcdf files is ''xarray''. We provide a code snippet below, doing the same thing as the snippets above.
<syntaxhighlight lang="python" line>
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt

# HERE WE OPEN THE NETCDF FILE
data = xr.open_dataset('diagfi.nc',
decode_times=False)

# HERE WE READ THE VARIABLES (1D OUTPUT)
Time=data['Time']
lat=data['latitude']
lon=data['longitude']
al=data['altitude']

# HERE WE READ THE AREA (2D OUTPUT)
aire_GCM=data['aire']

# HERE WE READ 3D OUTPUTS
tsurf=data['tsurf'] # this is the surface temperature 3D field (time, latitude, longitude, altitude)

# HERE WE READ 4D OUTPUTS
temp=data['temp'] # this is the atmospheric temperature 4D field (time, latitude, longitude, altitude)

## let's take the time-averaged surface temperature
mean_tsurf = np.mean(tsurf,axis=0)

##Let's plot a lon-lat map
fig = plt.figure()
plt.contourf(lon,lat,mean_tsurf)
plt.colorbar(label='Surface Temperature (K)')
plt.xlabel('Longitude ($^{\circ}$)')
plt.ylabel('Latitude ($^{\circ}$)')
plt.show()
</syntaxhighlight>
Don't hesitate to use the function called ''.values'' to transform any ''xarray'' into a numpy array, especially in case of calculation time problems.
For more examples on how to use ''xarray'', take a look at the [https://docs.xarray.dev/en/stable/index.html documentation].
Here is another example of how one can use xarray with multiples netcdfiles.

<syntaxhighlight lang="python" line>

import xarray as xr
import os

# your folder where output files are stored
FOLDER = './your_folder_with_output_files/'

# take back the files from your FOLDER
list_files_folder=os.listdir(FOLDER)

# If there are several files.
# Sort your simulation files by date,
# so beginning of simulation will be top of the list
# and end of simulation will be end of the list.
list_files_folder.sort()

files = [FOLDER+str(f) for f in list_files_folder]
# if you want to keep only files of special_year you can add this option :
# files = [FOLDER+str(f) for f in list_files_folder if f.startswith("special_year")]

# xarray will magically concatenate your outfiles by 'Time' (or any other 'concat_dime' you want)
nc=xr.open_mfdataset(files,decode_times=False, concat_dim='Time', combine='nested')

# to check your keys
for key in nc.keys():
print(key)

# to load keys (example here with keys for a mesoscale simulation)
Times = nc['Times'][:]
PTOT = nc['PTOT'][:]
T = nc['T'][:]
W = nc['W'][:]

# you can use some functions to make averages etc

T_moy = T.mean(dim=['Time','south_north','west_east'])

# other functions
# .cumsum (cumulative sum)
# .rename (change the name of the object)

</syntaxhighlight>

One cool thing about xarray is that it is well optimized, and can do whatever you want to do on your data, but better than you. See for instance, the example below to plot a temperature lon-lat map. Xarray handles it in 5 lines of code, where you would need a lot more to set-up you plot in traditional matplotlib. And the results look almost good enough for a paper plot.

<syntaxhighlight lang="python" line>
import xarray as xr
import matplotlib.pyplot as plt
##Load your data and print it
file = '/home/lteinturier/Documents/PhD/wasp43b/chemistry_project/input_5xsolar.nc'
data = xr.open_dataset(file,decode_times=False)
print(data)
##extracting the altitude level #20 for the whole file
data = data.isel(altitude=20)
##let's assume that data hold a time-series of the temperature. Let's average it in time
temp = data['temp'].mean("Time",keep_attrs=True) #we keep the attribute when averaging to conserve the DataArray structure
##now we plot
fig = temp.plot.contourf(cmap='gnuplot',levels=50) #choose the colormap and the number of contourf levels
fig.ax.set_title("P = {:.2e} mbar".format(data.p.mean().values)) ##set-up your title. If you don't change it, the title will be the altitude in km of your atmospheric level
plt.show()

</syntaxhighlight>
This is only a fraction of what Xarray can do. Check the [https://docs.xarray.dev/en/stable/user-guide/index.html documentation] for more.

==== Python tutorials to make pretty visuals ====

We provide a tutorial on how to make pretty visuals using Generic PCM 3-D simulations [https://lmdz-forge.lmd.jussieu.fr/mediawiki/Planets/index.php/Plots_With_PyVista here].

=== Planetoplot ===

Planetoplot is a in-house, python based library developped to vizualize Generic PCM data.

The code and documentation can be found at: https://nbviewer.org/github/aymeric-spiga/planetoplot/blob/master/tutorial/planetoplot_tutorial.ipynb

[[Category:Generic-Model]]
[[Category:Generic-LMDZ]]
[[Category:Generic-DYNAMICO]]
[[Category:Generic-WRF]]
[[Category:Generic-1D]]

Tool Box

2023-12-07T10:51:28Z

LTeinturier: /* XARRAY python library (more modern) */

== Pre-processing Tools ==
=== newstart: a fortran program to modify start files ===

Newstart is an interactive tool to modify the start files (''start.nc'' and ''startfi.nc'').

To be usable, ''newstart'' should be compile in the ''LMDZ.COMMON'' directory by using the following command line:
<syntaxhighlight lang="bash">
./makelmdz_fcm -arch my_arch_file -p std -d 64x48x30 newstart
</syntaxhighlight>
In the example, my_arch_file is the name the arch files (see [https://lmdz-forge.lmd.jussieu.fr/mediawiki/Planets/index.php/The_Target_Architecture_(%22arch%22)_Files arch] ) and 64x48x30 is the resolution of the physical grid.
Then copy the executable from the ''LMDZ.COMMON/bin'' directory to your bench directory.

When you execute newstart, you can use both a ''start2archive'' file or the start files (''start.nc'' and ''startfi.nc''). Then the interactive interface will propose to modify several physical quantities such as the gravity, the surface pressure or the rotation of the planet. At the end of the procedure, two files are created: '' '''re'''start.nc'' and '' '''re'''startfi.nc''. They can be renamed and used as start files to initialize a new simulation.

We have prepared a simple tutorial to learn how to modify ''start.nc'' and ''startfi.nc'' files (= the the initial conditions) for the Generic PCM: https://lmdz-forge.lmd.jussieu.fr/mediawiki/Planets/index.php/Modify_start_Files

=== start2archive ===

The start2archive tool is similar to newstart in the sense that it can be used to modify the start files. But start2archive can modify the resolution of the physical grid, the topography and the surface thermal inertia while newstart cannot. It is also useful to create an archive of different starting states, then extractable as start files.
The command line to compile start2archive is similar to the one used for newstart:
<syntaxhighlight lang="bash">
./makelmdz_fcm -arch my_arch_file -p std -d 64x48x30 start2archive
</syntaxhighlight>
To modify the resolution, you should first create a start_archive (by using start2archive) file at the used resolution, then compile a newstart file at the new resolution. Newstart will interpolate all the physical quantities on the new grid.

=== other third party scripts and tools ===

You can easily modify start.nc and startfi.nc netcdf files with the xarray Python library. Below you can find an easy example where we modify the surface temperature field.

<syntaxhighlight lang="python">
from numpy import *
import numpy as np
import matplotlib.pyplot as mpl
import math
import xarray as xr

#FIRST WE GET THE DATA FROM THE GCM SIMULATION
nc = xr.open_dataset('startfi.nc',decode_times=False) # can be any netcdf file (e.g. start/startfi.nc files)

# BELOW PHYSICAL VARIABLES
physical_points=nc['physical_points']
lat=nc['latitude']
lon=nc['longitude']
aire_GCM=nc['area']

# BELOW THE VARIABLE WE WANT TO UPDATE
tsurf=nc['tsurf']
new_tsurf = np.empty(len(physical_points))

# LOOP TO MODIFY THE VARIABLE
for i in range(0,len(physical_points),1):
new_tsurf[i]=300. # here you put whatever you want ; in this exemple, we assume an isothermal temperature distribution

nc['tsurf'].values = new_tsurf
nc.to_netcdf('restartfi.nc')

# SANITY CHECK PLOTS

fig = mpl.figure(1)
mpl.plot(physical_points,tsurf)
mpl.plot(physical_points,new_tsurf)
mpl.xlabel('GCM Physical Points')
mpl.ylabel('Tsurf (K)')
mpl.show()
</syntaxhighlight>

== Post-processing tools ==
=== zrecast ===

With this program you can recast atmospheric (i.e.: 4D-dimentional longitude-latitude-altitude-time) data from
GCM outputs (e.g. as given in diagfi.nc files) onto either ''pressure'' or ''altitude'' above ''areoid vertical'' coordinates.
Since integrating the hydrostatic equation is required to recast the data, the input file must contain surface pressure
and atmospheric temperature, as well as the ground geopotential.
If recasting data onto ''pressure'' coordinates, then the output file name is given by the input file name to which ''_P.nc'' will be appened. If recasting data onto altitude above areoid coordinates, then a ''_A.nc'' will be appened.

=== mass stream function ===

The mass stream function (and the total angular momentum) can be computed from a diagfi.nc or a stats.nc, using the '''streamfunction.F90''' script. The script is located at

<syntaxhighlight lang="bash">
trunk/LMDZ.GENERIC/utilities
</syntaxhighlight>

To compile the script, open the ''compile'' file in the same directory and do the following:
* Replace "pgf90" with your favorite fortran compiler
* replace "/distrib/local/netcdf/pgi_7.1-6_32/lib" with the lib address and directory that contains your NetCDF library (file ''libnetcdf.a'').
* Replace "/distrib/local/netcdf/pgi_7.1-6_32/include" with the address of the directory that contains the NetCDF include file (''netcdf.inc'').
* You can mess with the compiling options but it is not mandatory.

Once the script is compiled, copy it in the same directory as your '''.nc''' file and run

<syntaxhighlight lang="bash">
./streamfunction.e
</syntaxhighlight>
The script will ask you for the name of your '''.nc''' file, and will run and produce a new '''nameofyourfile_stream.nc''' file.

'''Be careful''' : In this new file, all fields are temporally and zonally averaged.

If you want to use '''python''' instead of '''fortran''', you can take a look at this [https://github.com/aymeric-spiga/dynanalysis repo]. It hosts a tool to perform dynamical analysis of GCM simulations (and therefore, it computes the mass stream function and a lot of other stuff), but it is tailored for Dynamico only. This repo also takes care of recasting (it does the job of both ''zrecast.F90'' and ''streamfunction.F90'')

== Continuing Simulations ==

=== manually ===

At the end of a simulation, the model generates restart files (files 'restart.nc' and 'restartfi.nc') which contain the final state of the model. The 'restart.nc' and 'restartfi.nc' files have the same format as the 'start.nc' and 'startfi.nc' files, respectively.

These files can in fact be used as initial states to continue the simulation, using the following renaming command lines:

<syntaxhighlight lang="bash">
mv restart.nc start.nc
mv restartfi.nc startfi.nc
</syntaxhighlight>

Running a simulation with these start files will in fact resume the simulation from where the previous run ended.

=== with bash scripts ===

We have set up very simple bash scripts to automatize the launching of chain simulations. Here is an example of bash script that does the job:

<syntaxhighlight lang="bash">
#!/bin/bash
###########################################################################
# Script to perform several chained LMD Mars GCM simulations
# SET HERE the maximum total number of simulations

nummax=100

###########################################################################

echo "---------------------------------------------------------"
echo "STARTING LOOP RUN"
echo "---------------------------------------------------------"

dir=`pwd`
machine=`hostname`
address=`whoami`

# Look for file "num_run" which should contain
# the value of the previously computed season
# (defaults to 0 if file "num_run" does not exist)
if [[ -r num_run ]] ; then
echo "found file num_run"
numold=`cat num_run`
else
numold=0
fi
echo "numold is set to" ${numold}

# Set value of current season
(( numnew = ${numold} + 1 ))
echo "numnew is set to" ${numnew}

# Look for initialization data files (exit if none found)
if [[ ( -r start${numold}.nc && -r startfi${numold}.nc ) ]] ; then
\cp -f start${numold}.nc start.nc
\cp -f startfi${numold}.nc startfi.nc
else
if (( ${numold} == 99999 )) ; then
echo "No run because previous run crashed ! (99999 in num_run)"
exit
else
echo "Where is file start"${numold}".nc??"
exit
fi
fi

# Run GCM -- THIS LINE NEEDS TO BE MODIFIED WITH THE CORRECT GCM EXECUTION COMMAND
mpirun -np 8 gcm_64x48x26_phystd_para.e < diagfi.def > lrun${numnew}

# Check if run ended normaly and copy datafiles
if [[ ( -r restartfi.nc && -r restart.nc ) ]] ; then
echo "Run seems to have ended normaly"

\mv -f restart.nc start${numnew}.nc
\mv -f restartfi.nc startfi${numnew}.nc

else
if [[ -r num_run ]] ; then
\mv -f num_run num_run.crash
else
echo "No file num_run to build num_run.crash from !!"
# Impose a default value of 0 for num_run
echo 0 > num_run.crash
fi
echo 99999 > num_run
exit
fi

# Copy other datafiles that may have been generated
if [[ -r diagfi.nc ]] ; then
\mv -f diagfi.nc diagfi${numnew}.nc
fi
if [[ -r diagsoil.nc ]] ; then
\mv -f diagsoil.nc diagsoil${numnew}.nc
fi
if [[ -r stats.nc ]] ; then
\mv -f stats.nc stats${numnew}.nc
fi
if [[ -f profiles.dat ]] ; then
\mv -f profiles.dat profiles${numnew}.dat
\mv -f profiles.hdr profiles${numnew}.hdr
fi

# Prepare things for upcoming runs by writing
# value of computed season in file num_run
echo ${numnew} > num_run

# If we are over nummax : stop
if (( $numnew + 1 > $nummax )) ; then
exit
# If not : restart the loop (copy the executable and run the copy)
else
\cp -f run_gnome exe_mars
./exe_mars
fi

</syntaxhighlight>

'''Summary of what this bash script does''':

* It reads the file 'num_run' which contains the step of the simulation.
If num_run is
<syntaxhighlight lang="bash">
5
</syntaxhighlight>
then the script expects to read start5.nc and startfi5.nc.
* It modifies start5.nc and startfi5.nc into start.nc and startfi.nc, respectively.
* It runs the GCM.
* It modifies restart.nc and restartfi.nc into start6.nc and startfi6.nc
* It rewrite num_run as follows:
<syntaxhighlight lang="bash">
6
</syntaxhighlight>
* It restarts the loop until num_run reaches the value (defined in nummax):
<syntaxhighlight lang="bash">
100
</syntaxhighlight>

== Processing Output Files with NCOs ==

NCOs (NetCdf OperatorS) are a set of powerful command-line utilities – available on Linux, Mac and PC – that allow to perform useful (and very fast!) post-processing operations on netCDF GCM output files. Full documentation can be found on http://research.jisao.washington.edu/data_sets/nco/, but we provide below a few examples of command lines.

* How to calculate a time mean of a netCDF 'diagfi.nc' file

<syntaxhighlight lang="bash">
ncra -F -d Time,1,,1 diagfi.nc diagfi_MEAN.nc # format is "-d dimension,minimum,maximum,stride"
</syntaxhighlight>

* Subsetting time in a netCDF 'diagfi.nc' file

<syntaxhighlight lang="bash">
ncea -F -d Time,first,last diagfi.nc diagfi_subset.nc # format is "-d dimension,minimum,maximum" ; we recall you can type "ncdump -v time diagfi.nc" to see the Time values in the netCDF file.
</syntaxhighlight>

* Decimating a netCDF 'diagfi.nc' file in time

<syntaxhighlight lang="bash">
ncks -F -d Time,1,,8 diagfi.nc diagfi_decimated.nc # format is "-d dimension,minimum,maximum,stride" ; In this example, this means that data is extracted 1 time every 8 time steps, starting from the first time step (number 1), ending at the last time step).
</syntaxhighlight>

* Extract a variable from a netCDF 'diagfi.nc' file

<syntaxhighlight lang="bash">
ncks -v tsurf,temp,p diagfi.nc diagfi_out.nc # Here we created a new file named 'diagfi_out.nc' in which we only kept variables named 'tsurf' (surface temperatures), 'temp' (atmospheric temperatures) and p (atmospheric pressures).
</syntaxhighlight>

Again, more examples can be found on http://research.jisao.washington.edu/data_sets/nco/ .

== Data Handling and Visualization Software ==

There are several data handling and visualization tools that can be used to analyse and plot the results from GCM simulations (using the diagfi.nc NetCDF files). We provide below a panorama of most widely used solutions.

=== panoply ===
Panoply is a user-friendly tool for viewing raw NetCDF data, available here: https://www.giss.nasa.gov/tools/panoply/ . It is very convenient to make pretty visuals (see an example for the exoplanet TRAPPIST-1e). There are many options that can be used (map projections, masks, colorbars, shadows, etc.) to make your plots fancy. However, the tool is not very well suited for manipulating data (compute averages, statistics, etc.).

* Installation on Linux:
You simply need to download and untar the Package from the Panoply website. Note that to work it requires that Java and related Java Runtime environment (JRE) be installed on your system (otherwise it will simply look as if "nothing is happening" when you try to launch Panoply via the "panoply.sh" script), which on Ubuntu simply requires something like:
<syntaxhighlight lang="bash">
sudo apt install java
sudo apt install default-jre
</syntaxhighlight>

* Run on Linux (assuming the panoply.sh script is in a directory included in your PATH environment variable):
<syntaxhighlight lang="bash">
panoply.sh
</syntaxhighlight>

[[File:Example panoply.png|thumb|Screenshot of panoply showing here Generic PCM results for the exoplanet TRAPPIST-1e (surface temperatures)]]

=== ncview ===
ncview is another useful user-friendly tool for viewing raw NetCDF data. This is kind of a very archaic version of panoply, but it is convenient because it allows to have a very quick first look at netCDF data files.

Command line tool to visualize NetCDF data:
* Installation on Linux-Ubuntu:
<syntaxhighlight lang="bash">
sudo apt install ncview
</syntaxhighlight>
* Run on Linux:
<syntaxhighlight lang="bash">
ncview diagfi.nc
</syntaxhighlight>

[[File:Example ncview.png|thumb|Screenshot of ncview showing here Generic PCM results for the exoplanet Proxima b (OLR - Thermal Emission)]]

=== python scripts ===

Python scripts provide a very useful mean to analyse and visualize netCDF files.

==== NETCDF4 python library (old school) ====

You can use the netCDF4 python library to open a netCDF file and put data in tables that can then be manipulated and plotted.

Here is an exemple of how to open and read a netCDF file with Python:

<syntaxhighlight lang="python" line>
import numpy
from netCDF4 import Dataset

# HERE WE OPEN THE NETCDF FILE
nc = Dataset('diagfi.nc')

# HERE WE READ THE VARIABLES (1D OUTPUT)
Time=nc.variables['Time'][:]
lat=nc.variables['latitude'][:]
lon=nc.variables['longitude'][:]
al=nc.variables['altitude'][:]

# HERE WE READ THE AREA (2D OUTPUT)
aire_GCM=nc.variables['aire'][:][:]

# HERE WE READ 3D OUTPUTS
tsurf=nc.variables['tsurf'][:][:][:] # this is the surface temperature 3D field (time, latitude, longitude, altitude)

# HERE WE READ 4D OUTPUTS
temp=nc.variables['temp'][:][:][:][:] # this is the atmospheric temperature 4D field (time, latitude, longitude, altitude)

</syntaxhighlight>

And here is an exemple of how to manipulate the netCDF data (here to compute the time averaged surface temperatures):

<syntaxhighlight lang="python" line>
from numpy import *
import numpy as np

mean_tsurf=np.zeros((len(lat),len(lon)),dtype='f')

for i in range(0,len(Time)):
for j in range(0,len(lat)):
for k in range(0,len(lon)):
mean_tsurf[j,k]=mean_tsurf[j,k]+tsurf[i,j,k]*(1./len(Time))

</syntaxhighlight>

And here is a last exemple of how to plot the data (using matplotlib):

<syntaxhighlight lang="python" line>
import matplotlib.pyplot as plt

plt.figure(1)
plt.contourf(lon_GCM,lat_GCM,mean_tsurf)
plt.colorbar(label='Surface Temperature (K)')
plt.xlabel('Longitude ($^{\circ}$)')
plt.ylabel('Latitude ($^{\circ}$)')
plt.show()

</syntaxhighlight>

==== XARRAY python library (more modern) ====

Another useful library to deal with netcdf files is ''xarray''. We provide a code snippet below, doing the same thing as the snippets above.
<syntaxhighlight lang="python" line>
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt

# HERE WE OPEN THE NETCDF FILE
data = xr.open_dataset('diagfi.nc',
decode_times=False)

# HERE WE READ THE VARIABLES (1D OUTPUT)
Time=data['Time']
lat=data['latitude']
lon=data['longitude']
al=data['altitude']

# HERE WE READ THE AREA (2D OUTPUT)
aire_GCM=data['aire']

# HERE WE READ 3D OUTPUTS
tsurf=data['tsurf'] # this is the surface temperature 3D field (time, latitude, longitude, altitude)

# HERE WE READ 4D OUTPUTS
temp=data['temp'] # this is the atmospheric temperature 4D field (time, latitude, longitude, altitude)

## let's take the time-averaged surface temperature
mean_tsurf = np.mean(tsurf,axis=0)

##Let's plot a lon-lat map
fig = plt.figure()
plt.contourf(lon,lat,mean_tsurf)
plt.colorbar(label='Surface Temperature (K)')
plt.xlabel('Longitude ($^{\circ}$)')
plt.ylabel('Latitude ($^{\circ}$)')
plt.show()
</syntaxhighlight>
Don't hesitate to use the function called ''.values'' to transform any ''xarray'' into a numpy array, especially in case of calculation time problems.
For more examples on how to use ''xarray'', take a look at the [https://docs.xarray.dev/en/stable/index.html documentation].
Here is another example of how one can use xarray with multiples netcdfiles.

<syntaxhighlight lang="python" line>

import xarray as xr
import os

# your folder where output files are stored
FOLDER = './your_folder_with_output_files/'

# take back the files from your FOLDER
list_files_folder=os.listdir(FOLDER)

# If there are several files.
# Sort your simulation files by date,
# so beginning of simulation will be top of the list
# and end of simulation will be end of the list.
list_files_folder.sort()

files = [FOLDER+str(f) for f in list_files_folder]
# if you want to keep only files of special_year you can add this option :
# files = [FOLDER+str(f) for f in list_files_folder if f.startswith("special_year")]

# xarray will magically concatenate your outfiles by 'Time' (or any other 'concat_dime' you want)
nc=xr.open_mfdataset(files,decode_times=False, concat_dim='Time', combine='nested')

# to check your keys
for key in nc.keys():
print(key)

# to load keys (example here with keys for a mesoscale simulation)
Times = nc['Times'][:]
PTOT = nc['PTOT'][:]
T = nc['T'][:]
W = nc['W'][:]

# you can use some functions to make averages etc

T_moy = T.mean(dim=['Time','south_north','west_east'])

# other functions
# .cumsum (cumulative sum)
# .rename (change the name of the object)

</syntaxhighlight>

One cool thing about xarray is that it is well optimized, and can do whatever you want to do on your data, but better than you. See for instance, the example below to plot a temperature lon-lat map. Xarray handles it in 5 lines of code, where you would need a lot more to set-up you plot in traditional matplotlib. And the results look almost good enough for a paper plot.

<syntaxhighlight lang="python" line>
import xarray as xr
import matplotlib.pyplot as plt
##Load your data and print it
file = '/home/lteinturier/Documents/PhD/wasp43b/chemistry_project/input_5xsolar.nc'
data = xr.open_dataset(file,decode_times=False)
print(data)
##extracting the altitude level # for the whole file
data = data.isel(altitude=20)
##let's assume that data hold a time-series of the temperature. Let's average it in time
temp = data['temp'].mean("Time",keep_attrs=True) #we keep the attribute when averaging
##now we plot
fig = temp.plot.contourf(cmap='gnuplot',levels=50) #choose the colormap and the number of contourf levels
fig.ax.set_title("P = {:.2e} mbar".format(data.p.mean().values)) ##set-up your title. If you don't change it, the title will be the altitude in km of your atmospheric level
plt.show()

</syntaxhighlight>
This is only a fraction of what Xarray can do. Check the [https://docs.xarray.dev/en/stable/user-guide/index.html documentation] for more.

==== Python tutorials to make pretty visuals ====

We provide a tutorial on how to make pretty visuals using Generic PCM 3-D simulations [https://lmdz-forge.lmd.jussieu.fr/mediawiki/Planets/index.php/Plots_With_PyVista here].

=== Planetoplot ===

Planetoplot is a in-house, python based library developped to vizualize Generic PCM data.

The code and documentation can be found at: https://nbviewer.org/github/aymeric-spiga/planetoplot/blob/master/tutorial/planetoplot_tutorial.ipynb

[[Category:Generic-Model]]
[[Category:Generic-LMDZ]]
[[Category:Generic-DYNAMICO]]
[[Category:Generic-WRF]]
[[Category:Generic-1D]]

Radiative Generic Condensable Specie

2023-03-02T09:50:50Z

LTeinturier: /* Radiative Generic Condensable Specie */

== Radiative Generic Condensable Specie ==
The model can take into account the radiative effects of Generic condensable specie in the radiative transfer. This is only tested for Hot Jupiter simulations for now.
Here is just a copy of the merge request to include the code into the master branch. This will need further editing:

We can activate the scheme by putting aerogeneric = # of aerosols in callphys.def. This is the only needed thing for activating the radiative effects. They must be tracer in modern tracer.def

Commented out the abort if we use more than 4 aerosols

Added reading of optprop files for Radiative Generic Condensable Aerosols

We use the same file for IR and VI channel. For now, only MnS, Cr,Fe and Mg2SiO4 can be read. If you want to add another specie, check the code, it is explained how to quickly do that (right above the Initializations)

Added radii calculation for Radiative Generic Condensable aerosols

Changed the hardcoded size of the totalemit array

The hardcoded size is now 1900 instead of 100 so we don't exceed the array size when working at high spectral resolution (very rare case)

Added opacity computation for Radiative Generic Condensable aerosols

We do this computation in the same fashion as what's performed on water and dust.

switch iniaerosol and initracer order, to prepare for the RGCS scheme

Needed to switch the order of initialization so we can use the RGCS scheme without the assumption that ice and vap tracer of the same specie are following each other in the traceur.def file
All my commits regarding the radiative effects of generic tracers. A lot of changes have been made. I've rebased the branch on the latest revision of master + squashed a lot of commits together to avoid the "avalanche de commits". By default, the scheme is NOT activated. Read above on how to activate it

==Adding stuff on how to really add a specie (because I'm once again stuck even though I coded that) ==

go to suaer_corrk. There, need to add the reading of the optical properties file.

Check that in your datadir, the table_tracers_condensable containe the species you wan. You need to have delta_vapH Tref, Pref, mass and metallicity_coeff.
In you have all that, you need your traceur.def with the mmol (actually we do't since it's in the table ? check that because weird), the radius and is_condensable=1, is_Rgcs=1

Tool Box

2022-11-23T14:12:50Z

LTeinturier: /* python scripts */

== Pre-processing Tools ==
=== newstart: a fortran program to modify start files ===

Newstart is an interactive tool to modify the start files (''start.nc'' and ''startfi.nc'').

To be usable, ''newstart'' should be compile in the ''LMDZ.COMMON'' directory by using the following command line:
<syntaxhighlight lang="bash">
./makelmdz_fcm -arch my_arch_file -p std -d 64x48x30 newstart
</syntaxhighlight>
In the example, my_arch_file is the name the arch files (see [https://lmdz-forge.lmd.jussieu.fr/mediawiki/Planets/index.php/The_Target_Architecture_(%22arch%22)_Files arch] ) and 64x48x30 is the resolution of the physical grid.
Then copy the executable from the ''LMDZ.COMMON/bin'' directory to your bench directory.

When you execute newstart, you can use both a ''start2archive'' file or the start files (''start.nc'' and ''startfi.nc''). Then the interactive interface will propose to modify several physical quantities such as the gravity, the surface pressure or the rotation of the planet. At the end of the procedure, two files are created: ''restart.nc'' and ''restartfi.nc''. They can be renamed and used as start files to initialize a new simulation.

=== start2archive ===

The start2archive tool is similar to newstart in the sense that it can be used to modify the start files. But start2archive can modify the resolution of the physical grid, the topography and the surface thermal inertia while newstart cannot. It is also useful to create an archive of different starting states, then extractable as start files.
The command line to compile start2archive is similar to the one used for newstart:
<syntaxhighlight lang="bash">
./makelmdz_fcm -arch my_arch_file -p std -d 64x48x30 start2archive
</syntaxhighlight>
To modify the resolution, you should first create a start_archive (by using start2archive) file at the used resolution, then compile a newstart file at the new resolution. Newstart will interpolate all the physical quantities on the new grid.

=== other third party scripts and tools ===

TO BE COMPLETED

== Post-processing tools ==
=== zrecast ===

With this program you can recast atmospheric (i.e.: 4D-dimentional longitude-latitude-altitude-time) data from
GCM outputs (e.g. as given in diagfi.nc files) onto either ''pressure'' or ''altitude'' above ''areoid vertical'' coordinates.
Since integrating the hydrostatic equation is required to recast the data, the input file must contain surface pressure
and atmospheric temperature, as well as the ground geopotential.
If recasting data onto ''pressure'' coordinates, then the output file name is given by the input file name to which ''_P.nc'' will be appened. If recasting data onto altitude above areoid coordinates, then a ''_A.nc'' will be appened.

=== mass stream function ===

The mass stream function (and the total angular momentum) can be computed from a diagfi.nc or a stats.nc, using the '''streamfunction.F90''' script. The script is located at

<syntaxhighlight lang="bash">
trunk/LMDZ.GENERIC/utilities
</syntaxhighlight>

To compile the script, open the ''compile'' file in the same directory and do the following:
* Replace "pgf90" with your favorite fortran compiler
* replace "/distrib/local/netcdf/pgi_7.1-6_32/lib" with the lib address and directory that contains your NetCDF library (file ''libnetcdf.a'').
* Replace "/distrib/local/netcdf/pgi_7.1-6_32/include" with the address of the directory that contains the NetCDF include file (''netcdf.inc'').
* You can mess with the compiling options but it is not mandatory.

Once the script is compiled, copy it in the same directory as your '''.nc''' file and run

<syntaxhighlight lang="bash">
./streamfunction.e
</syntaxhighlight>
The script will ask you for the name of your '''.nc''' file, and will run and produce a new '''nameofyourfile_stream.nc''' file.

'''Be careful''' : In this new file, all fields are temporally and zonally averaged.

If you want to use '''python''' instead of '''fortran''', you can take a look at this [https://github.com/aymeric-spiga/dynanalysis repo]. It hosts a tool to perform dynamical analysis of GCM simulations (and therefore, it computes the mass stream function and a lot of other stuff), but it is tailored for Dynamico only. This repo also takes care of recasting (it does the job of both ''zrecast.F90'' and ''streamfunction.F90'')

== Continuing Simulations ==

=== manually ===

At the end of a simulation, the model generates restart files (files 'restart.nc' and 'restartfi.nc') which contain the final state of the model. The 'restart.nc' and 'restartfi.nc' files have the same format as the 'start.nc' and 'startfi.nc' files, respectively.

These files can in fact be used as initial states to continue the simulation, using the following renaming command lines:

<syntaxhighlight lang="bash">
mv restart.nc start.nc
mv restartfi.nc startfi.nc
</syntaxhighlight>

Running a simulation with these start files will in fact resume the simulation from where the previous run ended.

=== with bash scripts ===

We have set up very simple bash scripts to automatize the launching of chain simulations. Here is an example of bash script that does the job:

<syntaxhighlight lang="bash">
#!/bin/bash
###########################################################################
# Script to perform several chained LMD Mars GCM simulations
# SET HERE the maximum total number of simulations

nummax=100

###########################################################################

echo "---------------------------------------------------------"
echo "STARTING LOOP RUN"
echo "---------------------------------------------------------"

dir=`pwd`
machine=`hostname`
address=`whoami`

# Look for file "num_run" which should contain
# the value of the previously computed season
# (defaults to 0 if file "num_run" does not exist)
if [[ -r num_run ]] ; then
echo "found file num_run"
numold=`cat num_run`
else
numold=0
fi
echo "numold is set to" ${numold}

# Set value of current season
(( numnew = ${numold} + 1 ))
echo "numnew is set to" ${numnew}

# Look for initialization data files (exit if none found)
if [[ ( -r start${numold}.nc && -r startfi${numold}.nc ) ]] ; then
\cp -f start${numold}.nc start.nc
\cp -f startfi${numold}.nc startfi.nc
else
if (( ${numold} == 99999 )) ; then
echo "No run because previous run crashed ! (99999 in num_run)"
exit
else
echo "Where is file start"${numold}".nc??"
exit
fi
fi

# Run GCM -- THIS LINE NEEDS TO BE MODIFIED WITH THE CORRECT GCM EXECUTION COMMAND
mpirun -np 8 gcm_64x48x26_phystd_para.e < diagfi.def > lrun${numnew}

# Check if run ended normaly and copy datafiles
if [[ ( -r restartfi.nc && -r restart.nc ) ]] ; then
echo "Run seems to have ended normaly"

\mv -f restart.nc start${numnew}.nc
\mv -f restartfi.nc startfi${numnew}.nc

else
if [[ -r num_run ]] ; then
\mv -f num_run num_run.crash
else
echo "No file num_run to build num_run.crash from !!"
# Impose a default value of 0 for num_run
echo 0 > num_run.crash
fi
echo 99999 > num_run
exit
fi

# Copy other datafiles that may have been generated
if [[ -r diagfi.nc ]] ; then
\mv -f diagfi.nc diagfi${numnew}.nc
fi
if [[ -r diagsoil.nc ]] ; then
\mv -f diagsoil.nc diagsoil${numnew}.nc
fi
if [[ -r stats.nc ]] ; then
\mv -f stats.nc stats${numnew}.nc
fi
if [[ -f profiles.dat ]] ; then
\mv -f profiles.dat profiles${numnew}.dat
\mv -f profiles.hdr profiles${numnew}.hdr
fi

# Prepare things for upcoming runs by writing
# value of computed season in file num_run
echo ${numnew} > num_run

# If we are over nummax : stop
if (( $numnew + 1 > $nummax )) ; then
exit
else
\cp -f run_gnome exe_mars
./exe_mars
fi

</syntaxhighlight>

'''Summary of what this bash script does''':

* It reads the file 'num_run' which contains the step of the simulation.
If num_run is
<syntaxhighlight lang="bash">
5
</syntaxhighlight>
then the script expects to read start5.nc and startfi5.nc.
* It modifies start5.nc and startfi5.nc into start.nc and startfi.nc, respectively.
* It runs the GCM.
* It modifies restart.nc and restartfi.nc into start6.nc and startfi6.nc
* It rewrite num_run as follows:
<syntaxhighlight lang="bash">
6
</syntaxhighlight>
* It restarts the loop until num_run reaches the value (defined in nummax):
<syntaxhighlight lang="bash">
100
</syntaxhighlight>

== Processing Output Files with NCOs ==

NCOs (NetCdf OperatorS) are a set of powerful command-line utilities – available on Linux, Mac and PC – that allow to perform useful (and very fast!) post-processing operations on netCDF GCM output files. Full documentation can be found on http://research.jisao.washington.edu/data_sets/nco/, but we provide below a few examples of command lines.

* How to calculate a time mean of a netCDF 'diagfi.nc' file

<syntaxhighlight lang="bash">
ncra -F -d Time,1,,1 diagfi.nc diagfi_MEAN.nc # format is "-d dimension,minimum,maximum,stride"
</syntaxhighlight>

* Subsetting time in a netCDF 'diagfi.nc' file

<syntaxhighlight lang="bash">
ncea -F -d Time,first,last diagfi.nc diagfi_subset.nc # format is "-d dimension,minimum,maximum" ; we recall you can type "ncdump -v time diagfi.nc" to see the Time values in the netCDF file.
</syntaxhighlight>

* Decimating a netCDF 'diagfi.nc' file in time

<syntaxhighlight lang="bash">
ncks -F -d Time,1,,8 diagfi.nc diagfi_decimated.nc # format is "-d dimension,minimum,maximum,stride" ; In this example, this means that data is extracted 1 time every 8 time steps, starting from the first time step (number 1), ending at the last time step).
</syntaxhighlight>

* Extract a variable from a netCDF 'diagfi.nc' file

<syntaxhighlight lang="bash">
ncks -v tsurf,temp,p diagfi.nc diagfi_out.nc # Here we created a new file named 'diagfi_out.nc' in which we only kept variables named 'tsurf' (surface temperatures), 'temp' (atmospheric temperatures) and p (atmospheric pressures).
</syntaxhighlight>

Again, more examples can be found on http://research.jisao.washington.edu/data_sets/nco/ .

== Data Handling and Visualization Software ==

There are several data handling and visualization tools that can be used to analyse and plot the results from GCM simulations (using the diagfi.nc NetCDF files). We provide below a panorama of most widely used solutions.

=== panoply ===
Panoply is a user-friendly tool for viewing raw NetCDF data, available here: https://www.giss.nasa.gov/tools/panoply/ . It is very convenient to make pretty visuals (see an example for the exoplanet TRAPPIST-1e). There are many options that can be used (map projections, masks, colorbars, shadows, etc.) to make your plots fancy. However, the tool is not very well suited for manipulating data (compute averages, statistics, etc.).

[[File:Example panoply.png|thumb|Screenshot of panoply showing here Generic PCM results for the exoplanet TRAPPIST-1e (surface temperatures)]]

=== ncview ===
ncview is another useful user-friendly tool for viewing raw NetCDF data. This is kind of a very archaic version of panoply, but it is convenient because it allows to have a very quick first look at netCDF data files.

Command line tool to visualize NetCDF data:
* Installation on Linux:
<syntaxhighlight lang="bash">
sudo apt install ncview
</syntaxhighlight>
* Run on Linux:
<syntaxhighlight lang="bash">
ncview diagfi.nc
</syntaxhighlight>

[[File:Example ncview.png|thumb|Screenshot of ncview showing here Generic PCM results for the exoplanet Proxima b (OLR - Thermal Emission)]]

=== python scripts ===

Python scripts provide a very useful mean to analyse and visualize netCDF files.

You can use the netCDF4 python library to open a netCDF file and put data in tables that can then be manipulated and plotted.

Here is an exemple of how to open and read a netCDF file with Python:

<syntaxhighlight lang="python" line>
import numpy
from netCDF4 import Dataset

# HERE WE OPEN THE NETCDF FILE
nc = Dataset('diagfi.nc')

# HERE WE READ THE VARIABLES (1D OUTPUT)
Time=nc.variables['Time'][:]
lat=nc.variables['latitude'][:]
lon=nc.variables['longitude'][:]
al=nc.variables['altitude'][:]

# HERE WE READ THE AREA (2D OUTPUT)
aire_GCM=nc.variables['aire'][:][:]

# HERE WE READ 3D OUTPUTS
tsurf=nc.variables['tsurf'][:][:][:] # this is the surface temperature 3D field (time, latitude, longitude, altitude)

# HERE WE READ 4D OUTPUTS
temp=nc.variables['temp'][:][:][:][:] # this is the atmospheric temperature 4D field (time, latitude, longitude, altitude)

</syntaxhighlight>

And here is an exemple of how to manipulate the netCDF data (here to compute the time averaged surface temperatures):

<syntaxhighlight lang="python" line>
from numpy import *
import numpy as np

mean_tsurf=np.zeros((len(lat_GCM),len(lon_GCM)),dtype='f')

for i in range(0,len(Time)):
for j in range(0,len(lat)):
for k in range(0,len(lon)):
mean_tsurf[j,k]=mean_tsurf[j,k]+tsurf[i,j,k]*(1./len(Time))

</syntaxhighlight>

And here is a last exemple of how to plot the data (using matplotlib):

<syntaxhighlight lang="python" line>
import matplotlib.pyplot as plt

plt.figure(1)
plt.contourf(lon_GCM,lat_GCM,mean_tsurf)
plt.colorbar(label='Surface Temperature (K)')
plt.xlabel('Longitude ($^{\circ}$)')
plt.ylabel('Latitude ($^{\circ}$)')
plt.show()

</syntaxhighlight>

Another useful library to deal with netcdf files is ''xarray''. We provide a code snippet below, doing the same thing as the snippets above.
<syntaxhighlight lang="python" line>
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt

# HERE WE OPEN THE NETCDF FILE
data = xr.open_dataset('diagfi.nc',
decode_times=False)

# HERE WE READ THE VARIABLES (1D OUTPUT)
Time=data['Time']
lat=data['latitude']
lon=data['longitude']
al=data['altitude']

# HERE WE READ THE AREA (2D OUTPUT)
aire_GCM=data['aire']

# HERE WE READ 3D OUTPUTS
tsurf=data['tsurf'] # this is the surface temperature 3D field (time, latitude, longitude, altitude)

# HERE WE READ 4D OUTPUTS
temp=data['temp'] # this is the atmospheric temperature 4D field (time, latitude, longitude, altitude)

## let's take the time-averaged surface temperature
mean_tsurf = np.mean(tsurf,axis=0)

##Let's plot a lon-lat map
fig = plt.figure()
plt.contourf(lon,lat,mean_tsurf)
plt.colorbar(label='Surface Temperature (K)')
plt.xlabel('Longitude ($^{\circ}$)')
plt.ylabel('Latitude ($^{\circ}$)')
plt.show()
</syntaxhighlight>
For more examples on how to use ''xarray'', take a look at the [https://docs.xarray.dev/en/stable/index.html documentation].
Here is another example of how one can use xarray with multiples netcdfiles.

<syntaxhighlight lang="python" line>

import xarray as xr
import os

# your folder where output files are stored
FOLDER = './your_folder_with_output_files/'

# take back the files from your FOLDER
list_files_folder=os.listdir(FOLDER)

# If there are several files.
# Sort your simulation files by date,
# so beginning of simulation will be top of the list
# and end of simulation will be end of the list.
list_files_folder.sort()

files = [FOLDER+str(f) for f in list_files_folder]
# if you want to keep only files of special_year you can add this option :
# files = [FOLDER+str(f) for f in list_files_folder if f.startswith("special_year")]

# xarray will magically concatenate your outfiles by 'Time' (or any other 'concat_dime' you want)
nc=xr.open_mfdataset(files,decode_times=False, concat_dim='Time', combine='nested')

# to check your keys
for key in nc.keys():
print(key)

# to load keys (example here with keys for a mesoscale simulation)
Times = nc['Times'][:]
PTOT = nc['PTOT'][:]
T = nc['T'][:]
W = nc['W'][:]

# you can use some functions to make averages etc

T_moy = T.mean(dim=['Time','south_north','west_east'])

# other functions
# .cumsum (cumulative sum)
# .rename (change the name of the object)

</syntaxhighlight>

NB: ADD SOMETHING HERE ABOUT XARRAY LIBRARY? AYMERIC? (with a simple tutorial?) TBD

We provide a tutorial on how to make pretty visuals using Generic PCM 3-D simulations [https://lmdz-forge.lmd.jussieu.fr/mediawiki/Planets/index.php/Plots_With_PyVista here].

=== Planetoplot ===

Planetoplot is a in-house, python based library developped to vizualize Generic PCM data.

The code and documentation can be found at: https://nbviewer.org/github/aymeric-spiga/planetoplot/blob/master/tutorial/planetoplot_tutorial.ipynb

Radiative Generic Condensable Specie

2022-10-19T11:48:25Z

LTeinturier: Created page with "== Radiative Generic Condensable Specie == The model can take into account the radiative effects of Generic condensable specie in the radiative transfer. This is only tested f..."

Other GCM Configurations worth knowing about

2022-09-02T09:47:04Z

LTeinturier: /* Hot Jupiter with DYNAMICO */

== 3D lon-lat LMDZ setup ==

=== early Mars ===

It is already described in the [https://lmdz-forge.lmd.jussieu.fr/mediawiki/Planets/index.php/Quick_Install_and_Run ''Quick Install and Run''] section.

=== Earth with slab ocean ===

TBD by Martin

=== TRAPPIST-1e with photochemistry ===

A temperate rocky planet in synchronous rotation around a low mass star

TBD by Yassin

=== TRAPPIST-1c in Venus-like conditions ===

A warm rocky planet in synchronous rotation around a low mass star

TBD by Gabriella (waiting for the SVN update by Ehouarn)

=== mini-Neptune GJ1214b ===

A warm mini-Neptune

TBD by Benjamin

== 3D DYNAMICO setup ==

Due to the rich dynamical activities in their atmospheres (banded zonal jets, eddies, vortices, storms, equatorial oscillations,...) resulting from multi-scale dynamic interactions, the Global Climate Modelling of the giant planet requires to resolve eddies arising from hydrodynamical instabilities to correctly establish the planetary-scaled jets regime. To this purpose, their Rossby radius deformation $$L_D$$, which is the length scale at which rotational effects become as important as buoyancy or gravity wave effects in the evolution of the flow about some disturbance, is calculated to determine the most suitable horizontal grid resolution. At mid-latitude range, for the giant planets, $$L_D$$ is of the same order of magnitude as that of the Earth. As the giant planets have a size of roughly 10 times the Earth size (i.e., Jupiter and Saturn), the modelling grid must be of a horizontal resolution of 0.5$$^{\circ}$$ over longitude and latitude (vs 5$$^{\circ}$$ for the Earth), considering 3 grid points to resolved $$L_D$$.
Moreover, to have a chance to model the equatorial oscillation, meridional cell circulations and/or a seasonal inter-hemispheric circulation, a giant planet GCM must also include a high vertical resolution. Indeed, these climate phenomena have been studied for decades for the Earth's atmosphere, and result from small- and large-scale interactions between the troposphere and stratosphere. This implies that the propagation of dynamic instabilities, waves and turbulence should be resolved as far as possible along the vertical. Contrary to horizontal resolution, it doesn't really exist a criterion (similar to $$L_D$$) to determine the most suitable vertical grid resolution and still an adjustable parameter according to the processes to be represented. However, we advise the user to set a vertical resolution of at least 5 grid points per scale height as first stage.
Finally, these atmospheres are cold, with long radiative response time which needs radiative transfer computations over decade-long years of Jupiter (given that a Jupiter year $$\approx$$ 12 Earth years), Saturn ( a Saturn year $$\approx$$ 30 Earth years), Uranus (a Uranus year $$\approx$$ 84 earth years) or Neptune (a Neptune year $$\approx$$ 169 Earth years), depending on the chosen planet.

To be able to deal with these three -- and non-exhaustive -- requirements to build a giant planet GCM, we need massive computational ressources. For this, we use a dynamical core suitable and numerical stable for massive parallel ressource computations: DYNAMICO [Dubos et al,. 2015].

In these two following subsections, we purpose an example of installation for Jupiter and a Hot Jupiter. All the install, compiling, setting and parameters files for each giant planets could be found on:
https://github.com/aymeric-spiga/dynamico-giant

If you have already downloaded '''LMDZ.COMMON''', '''LMDZ.GENERIC''', '''IOIPSL''', '''ARCH''', you only have to download:

'''ICOSAGCM''': the DYNAMICO dynamical core
<syntaxhighlight lang="bash">
git clone https://gitlab.in2p3.fr/ipsl/projets/dynamico/dynamico.git ICOSAGCM
cd ICOSAGCM
git checkout 90f7138a60ebd3644fbbc42bc9dfa22923386385
</syntaxhighlight>

'''ICOSA_LMDZ''': the interface using to link LMDZ.GENERIC physical packages and ICOSAGCM
<syntaxhighlight lang="bash">
svn update -r 2655 -q ICOSA_LMDZ
</syntaxhighlight>

'''XIOS (XML Input Output Server)''': the library to interpolate input/output fields between the icosahedral and longitude/latitude regular grids on fly
<syntaxhighlight lang="bash">
svn co -r 2319 -q http://forge.ipsl.jussieu.fr/ioserver/svn/XIOS/trunk XIOS
</syntaxhighlight>

If you haven't already download '''LMDZ.COMMON''', '''LMDZ.GENERIC''', '''IOIPSL''', '''ARCH''', you can use the '''install.sh''' script provided by the Github repository.

=== Jupiter with DYNAMICO ===

TBD by Deborah

=== Hot Jupiter with DYNAMICO ===

Modelling the atmosphere of Hot Jupiter is challenging because of the extreme temperature conditions, and the fact that these planets are gas giants. Therefore, using a dynamical core such as Dynamico is strongly recommended. Here, we discuss how to perform a cloudless simulation of the Hot Jupiter WASP-43 b, using Dynamico.

'''1st step''': You need to go to the github mentionned previously for Dynamico: https://github.com/aymeric-spiga/dynamico-giant. ''Git clone'' this repo on your favorite cluster, and ''checkout'' to the "hot_jupiter" branch.

'''2nd step''': Now, run the install.sh script. This script will install '''all''' the required models ('''LMDZ.COMMON''', '''LMDZ.GENERIC''','''ICOSA_LMDZ''','''XIOS''','''FCM''','''ICOSAGCM'''). At this point, you only miss '''IOIPSL'''. To install it, go to
<syntaxhighlight lang="bash">
dynamico-giant/code/LMDZ.COMMON/ioipsl/
</syntaxhighlight>

There, you will find some examples of installations script. You need to create one that will work on your cluster, with your own arch files.
During the installation of '''IOIPSL''', you might be asked for a login/password. Contact TGCC computing center to get access.

'''3rd step''': Great, now we have all we need to get started. Navigate to the ''hot_jupiter'' folder. You will find a ''compile_mesopsl.sh'' and a ''compile_occigen.sh'' script. Use them as examples to create the compile script adapted to your own cluster, then run it.
While running, I suggest that you take a look at the ''log_compile'' file. The compilation can take a while (~ 10minutes, especially because of XIOS). On quick trick to make sure that everything went right is to check the number of ''Build command finished'' messages in ''log_compile''. If everything worked out, there should be 6 of them.

'''4th step''': Okay, the model compiled, good job ! Now we need to create the initial condition for our run. In the hot_jupiter1d folder, you already have a ''temp_profile.txt'' computed with the 1D version of the LMDZ.GENERIC (see rcm1d on this page). Thus, no need to recompute a 1D model but it will be needed if you want to model another Hot Jupiter.
Navigate to the 'makestart' folder, located at
<syntaxhighlight lang="bash">
dynamico-giant/hot_jupiter/makestart/
</syntaxhighlight>
To generate the initial conditions for the 3D run, we're gonna start the model using the temperature profile from the 1D run. to do that, you will find a "job_mpi" script. Open it, and adapt it to your cluster and launch the job. This job is using 20 procs, and it runs 5 days of simulations.
If everything goes well, you should see few netcdf files appear. The important ones are '''start_icosa0.nc''', '''startfi0.nc''' and '''Xhistins.nc'''.
If you see these files, you're all set to launch a real simulation !

'''5th step''': Go back to ''hot_jupiter'' folder. There are a bunch of script to launch your simulation. Take a look at the ''astro_fat_mpi'' script, and adapt it to your cluster. Then you can launch your simulation by doing
<syntaxhighlight lang="bash">
./run_astro_fat
</syntaxhighlight>
This will start the simulation, using 90 procs. In the same folder, check if the icosa_lmdz.out file is created. This is the logfile of the simulation, while it is running. You can check there that everything is going well.

'''Important side note''': When using the ''run_astro_fat'' script to run a simulation, it will run a chained simulation, restarting the simulation from the previous state after 100 days of simulations and generating ''Xhistins.nc'' files. This is your results file, where you will find all the variables that controls your atmosphere (temperature field, wind fields, etc..).

Good luck and enjoy the generic PCM Dynamico for Hot Jupiter !

'''2nd important side note''': These 5 steps are the basic needed steps to run a simulation. If you want to tune simulations to another planet, or change other stuff, you need to take a look at '''*.def''' and '''*.xml''' files. If you're lost in all of this, take a look at the different pages of this website and/or contact us !
Also, you might want to check the wiki on the [https://github.com/aymeric-spiga/dynamico-giant ''Github''], that explains a lot of settings for Dynamico

== 3D LES setup ==

=== Proxima b with LES ===

TBD by Maxence

== 1D setup ==

=== rcm1d test case ===

Our 1-D forward model

TBD by Gwenael ? (you can have a look at the Generic GCM User Manual for inspiration)

=== kcm1d test case ===

Our 1-D inverse model

TBD by Guillaume or Martin

Other GCM Configurations worth knowing about

2022-09-02T09:27:13Z

LTeinturier: /* Hot Jupiter with DYNAMICO */

== 3D lon-lat LMDZ setup ==

=== early Mars ===

It is already described in the [https://lmdz-forge.lmd.jussieu.fr/mediawiki/Planets/index.php/Quick_Install_and_Run ''Quick Install and Run''] section.

=== Earth with slab ocean ===

TBD by Martin

=== TRAPPIST-1e with photochemistry ===

A temperate rocky planet in synchronous rotation around a low mass star

TBD by Yassin

=== TRAPPIST-1c in Venus-like conditions ===

A warm rocky planet in synchronous rotation around a low mass star

TBD by Gabriella (waiting for the SVN update by Ehouarn)

=== mini-Neptune GJ1214b ===

A warm mini-Neptune

TBD by Benjamin

== 3D DYNAMICO setup ==

Due to the rich dynamical activities in their atmospheres (banded zonal jets, eddies, vortices, storms, equatorial oscillations,...) resulting from multi-scale dynamic interactions, the Global Climate Modelling of the giant planet requires to resolve eddies arising from hydrodynamical instabilities to correctly establish the planetary-scaled jets regime. To this purpose, their Rossby radius deformation $$L_D$$, which is the length scale at which rotational effects become as important as buoyancy or gravity wave effects in the evolution of the flow about some disturbance, is calculated to determine the most suitable horizontal grid resolution. At mid-latitude range, for the giant planets, $$L_D$$ is of the same order of magnitude as that of the Earth. As the giant planets have a size of roughly 10 times the Earth size (i.e., Jupiter and Saturn), the modelling grid must be of a horizontal resolution of 0.5$$^{\circ}$$ over longitude and latitude (vs 5$$^{\circ}$$ for the Earth), considering 3 grid points to resolved $$L_D$$.
Moreover, to have a chance to model the equatorial oscillation, meridional cell circulations and/or a seasonal inter-hemispheric circulation, a giant planet GCM must also include a high vertical resolution. Indeed, these climate phenomena have been studied for decades for the Earth's atmosphere, and result from small- and large-scale interactions between the troposphere and stratosphere. This implies that the propagation of dynamic instabilities, waves and turbulence should be resolved as far as possible along the vertical. Contrary to horizontal resolution, it doesn't really exist a criterion (similar to $$L_D$$) to determine the most suitable vertical grid resolution and still an adjustable parameter according to the processes to be represented. However, we advise the user to set a vertical resolution of at least 5 grid points per scale height as first stage.
Finally, these atmospheres are cold, with long radiative response time which needs radiative transfer computations over decade-long years of Jupiter (given that a Jupiter year $$\approx$$ 12 Earth years), Saturn ( a Saturn year $$\approx$$ 30 Earth years), Uranus (a Uranus year $$\approx$$ 84 earth years) or Neptune (a Neptune year $$\approx$$ 169 Earth years), depending on the chosen planet.

To be able to deal with these three -- and non-exhaustive -- requirements to build a giant planet GCM, we need massive computational ressources. For this, we use a dynamical core suitable and numerical stable for massive parallel ressource computations: DYNAMICO [Dubos et al,. 2015].

In these two following subsections, we purpose an example of installation for Jupiter and a Hot Jupiter. All the install, compiling, setting and parameters files for each giant planets could be found on:
https://github.com/aymeric-spiga/dynamico-giant

If you have already downloaded '''LMDZ.COMMON''', '''LMDZ.GENERIC''', '''IOIPSL''', '''ARCH''', you only have to download:

'''ICOSAGCM''': the DYNAMICO dynamical core
<syntaxhighlight lang="bash">
git clone https://gitlab.in2p3.fr/ipsl/projets/dynamico/dynamico.git ICOSAGCM
cd ICOSAGCM
git checkout 90f7138a60ebd3644fbbc42bc9dfa22923386385
</syntaxhighlight>

'''ICOSA_LMDZ''': the interface using to link LMDZ.GENERIC physical packages and ICOSAGCM
<syntaxhighlight lang="bash">
svn update -r 2655 -q ICOSA_LMDZ
</syntaxhighlight>

'''XIOS (XML Input Output Server)''': the library to interpolate input/output fields between the icosahedral and longitude/latitude regular grids on fly
<syntaxhighlight lang="bash">
svn co -r 2319 -q http://forge.ipsl.jussieu.fr/ioserver/svn/XIOS/trunk XIOS
</syntaxhighlight>

If you haven't already download '''LMDZ.COMMON''', '''LMDZ.GENERIC''', '''IOIPSL''', '''ARCH''', you can use the '''install.sh''' script provided by the Github repository.

=== Jupiter with DYNAMICO ===

TBD by Deborah

=== Hot Jupiter with DYNAMICO ===

Modelling the atmosphere of Hot Jupiter is challenging because of the extreme temperature conditions, and the fact that these planets are gas giants. Therefore, using a dynamical core such as Dynamico is strongly recommended. Here, we discuss how to perform a cloudless simulation of the Hot Jupiter WASP-43 b, using Dynamico.

'''1st step''': You need to go to the github mentionned previously for Dynamico: https://github.com/aymeric-spiga/dynamico-giant. ''Git clone'' this repo on your favorite cluster, and ''checkout'' to the "hot_jupiter" branch.

'''2nd step''': Now, run the install.sh script. This script will install '''all''' the required models ('''LMDZ.COMMON''', '''LMDZ.GENERIC''','''ICOSA_LMDZ''','''XIOS''','''FCM''','''ICOSAGCM'''). At this point, you only miss '''IOIPSL'''. To install it, go to
<syntaxhighlight lang="bash">
dynamico-giant/code/LMDZ.COMMON/ioipsl/
</syntaxhighlight>

There, you will find some examples of installations script. You need to create one that will work on your cluster, with your own arch files.
During the installation of '''IOIPSL''', you might be asked for a login/password. Contact TGCC computing center to get access.

'''3rd step''': Great, now we have all we need to get started. Navigate to the ''hot_jupiter'' folder. You will find a ''compile_mesopsl.sh'' and a ''compile_occigen.sh'' script. Use them as examples to create the compile script adapted to your own cluster, then run it.
While running, I suggest that you take a look at the ''log_compile'' file. The compilation can take a while (~ 10minutes, especially because of XIOS). On quick trick to make sure that everything went right is to check the number of ''Build command finished'' messages in ''log_compile''. If everything worked out, there should be 6 of them.

'''4th step''': Okay, the model compiled, good job ! Now we need to create the initial condition for our run. In the hot_jupiter1d folder, you already have a ''temp_profile.txt'' computed with the 1D version of the LMDZ.GENERIC (see rcm1d on this page). Thus, no need to recompute a 1D model but it will be needed if you want to model another Hot Jupiter.
Navigate to the 'makestart' folder, located at
<syntaxhighlight lang="bash">
dynamico-giant/hot_jupiter/makestart/
</syntaxhighlight>
To generate the initial conditions for the 3D run, we're gonna start the model using the temperature profile from the 1D run. to do that, you will find a "job_mpi" script. Open it, and adapt it to your cluster and launch the job. This job is using 20 procs, and it runs 5 days of simulations.
If everything goes well, you should see few netcdf files appear. The important ones are '''start_icosa0.nc''', '''startfi0.nc''' and '''Xhistins.nc'''.
If you see these files, you're all set to launch a real simulation !

'''5th step''': Go back to ''hot_jupiter'' folder. There are a bunch of script to launch your simulation. Take a look at the ''astro_fat_mpi'' script, and adapt it to your cluster. Then you can launch your simulation by doing
<syntaxhighlight lang="bash">
./run_astro_fat
</syntaxhighlight>
This will start the simulation, using 90 procs. In the same folder, check if the icosa_lmdz.out file is created. This is the logfile of the simulation, while it is running. You can check there that everything is going well.

'''Important side note''': When using the ''run_astro_fat'' script to run a simulation, it will run a chained simulation, restarting the simulation from the previous state after 100 days of simulations and generating ''Xhistins.nc'' files. This is your results file, where you will find all the variables that controls your atmosphere (temperature field, wind fields, etc..).

Good luck and enjoy the generic PCM Dynamico for Hot Jupiter !

'''2nd important side note''': These 5 steps are the basic needed steps to run a simulation. If you want to tune simulations to another planet, or change other stuff, you need to take a look at '''*.def''' and '''*.xml''' files. If you're lost in all of this, take a look at the different pages of this website and/or contact us !

== 3D LES setup ==

=== Proxima b with LES ===

TBD by Maxence

== 1D setup ==

=== rcm1d test case ===

Our 1-D forward model

TBD by Gwenael ? (you can have a look at the Generic GCM User Manual for inspiration)

=== kcm1d test case ===

Our 1-D inverse model

TBD by Guillaume or Martin

Other GCM Configurations worth knowing about

2022-09-02T09:24:58Z

LTeinturier: /* Hot Jupiter with DYNAMICO */

== 3D lon-lat LMDZ setup ==

=== early Mars ===

It is already described in the [https://lmdz-forge.lmd.jussieu.fr/mediawiki/Planets/index.php/Quick_Install_and_Run ''Quick Install and Run''] section.

=== Earth with slab ocean ===

TBD by Martin

=== TRAPPIST-1e with photochemistry ===

A temperate rocky planet in synchronous rotation around a low mass star

TBD by Yassin

=== TRAPPIST-1c in Venus-like conditions ===

A warm rocky planet in synchronous rotation around a low mass star

TBD by Gabriella (waiting for the SVN update by Ehouarn)

=== mini-Neptune GJ1214b ===

A warm mini-Neptune

TBD by Benjamin

== 3D DYNAMICO setup ==

Due to the rich dynamical activities in their atmospheres (banded zonal jets, eddies, vortices, storms, equatorial oscillations,...) resulting from multi-scale dynamic interactions, the Global Climate Modelling of the giant planet requires to resolve eddies arising from hydrodynamical instabilities to correctly establish the planetary-scaled jets regime. To this purpose, their Rossby radius deformation $$L_D$$, which is the length scale at which rotational effects become as important as buoyancy or gravity wave effects in the evolution of the flow about some disturbance, is calculated to determine the most suitable horizontal grid resolution. At mid-latitude range, for the giant planets, $$L_D$$ is of the same order of magnitude as that of the Earth. As the giant planets have a size of roughly 10 times the Earth size (i.e., Jupiter and Saturn), the modelling grid must be of a horizontal resolution of 0.5$$^{\circ}$$ over longitude and latitude (vs 5$$^{\circ}$$ for the Earth), considering 3 grid points to resolved $$L_D$$.
Moreover, to have a chance to model the equatorial oscillation, meridional cell circulations and/or a seasonal inter-hemispheric circulation, a giant planet GCM must also include a high vertical resolution. Indeed, these climate phenomena have been studied for decades for the Earth's atmosphere, and result from small- and large-scale interactions between the troposphere and stratosphere. This implies that the propagation of dynamic instabilities, waves and turbulence should be resolved as far as possible along the vertical. Contrary to horizontal resolution, it doesn't really exist a criterion (similar to $$L_D$$) to determine the most suitable vertical grid resolution and still an adjustable parameter according to the processes to be represented. However, we advise the user to set a vertical resolution of at least 5 grid points per scale height as first stage.
Finally, these atmospheres are cold, with long radiative response time which needs radiative transfer computations over decade-long years of Jupiter (given that a Jupiter year $$\approx$$ 12 Earth years), Saturn ( a Saturn year $$\approx$$ 30 Earth years), Uranus (a Uranus year $$\approx$$ 84 earth years) or Neptune (a Neptune year $$\approx$$ 169 Earth years), depending on the chosen planet.

To be able to deal with these three -- and non-exhaustive -- requirements to build a giant planet GCM, we need massive computational ressources. For this, we use a dynamical core suitable and numerical stable for massive parallel ressource computations: DYNAMICO [Dubos et al,. 2015].

In these two following subsections, we purpose an example of installation for Jupiter and a Hot Jupiter. All the install, compiling, setting and parameters files for each giant planets could be found on:
https://github.com/aymeric-spiga/dynamico-giant

If you have already downloaded '''LMDZ.COMMON''', '''LMDZ.GENERIC''', '''IOIPSL''', '''ARCH''', you only have to download:

'''ICOSAGCM''': the DYNAMICO dynamical core
<syntaxhighlight lang="bash">
git clone https://gitlab.in2p3.fr/ipsl/projets/dynamico/dynamico.git ICOSAGCM
cd ICOSAGCM
git checkout 90f7138a60ebd3644fbbc42bc9dfa22923386385
</syntaxhighlight>

'''ICOSA_LMDZ''': the interface using to link LMDZ.GENERIC physical packages and ICOSAGCM
<syntaxhighlight lang="bash">
svn update -r 2655 -q ICOSA_LMDZ
</syntaxhighlight>

'''XIOS (XML Input Output Server)''': the library to interpolate input/output fields between the icosahedral and longitude/latitude regular grids on fly
<syntaxhighlight lang="bash">
svn co -r 2319 -q http://forge.ipsl.jussieu.fr/ioserver/svn/XIOS/trunk XIOS
</syntaxhighlight>

If you haven't already download '''LMDZ.COMMON''', '''LMDZ.GENERIC''', '''IOIPSL''', '''ARCH''', you can use the '''install.sh''' script provided by the Github repository.

=== Jupiter with DYNAMICO ===

TBD by Deborah

=== Hot Jupiter with DYNAMICO ===

Modelling the atmosphere of Hot Jupiter is challenging because of the extreme temperature conditions, and the fact that these planets are gas giants. Therefore, using a dynamical core such as Dynamico is strongly recommended. Here, we discuss how to perform a cloudless simulation of the Hot Jupiter WASP-43 b, using Dynamico.

'''1st step''': You need to go to the github mentionned previously for Dynamico: https://github.com/aymeric-spiga/dynamico-giant. ''Git clone'' this repo on your favorite cluster, and ''checkout'' to the "hot_jupiter" branch.

'''2nd step''': Now, run the install.sh script. This script will install '''all''' the required models ('''LMDZ.COMMON''', '''LMDZ.GENERIC''','''ICOSA_LMDZ''','''XIOS''','''FCM''','''ICOSAGCM'''). At this point, you only miss '''IOIPSL'''. To install it, go to
<syntaxhighlight lang="bash">
dynamico-giant/code/LMDZ.COMMON/ioipsl/
</syntaxhighlight>

There, you will find some examples of installations script. You need to create one that will work on your cluster, with your own arch files.
During the installation of '''IOIPSL''', you might be asked for a login/password. Contact TGCC computing center to get access.

'''3rd step''': Great, now we have all we need to get started. Navigate to the ''hot_jupiter'' folder. You will find a ''compile_mesopsl.sh'' and a ''compile_occigen.sh'' script. Use them as examples to create the compile script adapted to your own cluster, then run it.
While running, I suggest that you take a look at the ''log_compile'' file. The compilation can take a while (~ 10minutes, especially because of XIOS). On quick trick to make sure that everything went right is to check the number of ''Build command finished'' messages in ''log_compile''. If everything worked out, there should be 6 of them.

'''4th step''': Okay, the model compiled, good job ! Now we need to create the initial condition for our run. In the hot_jupiter1d folder, you already have a ''temp_profile.txt'' computed with the 1D version of the LMDZ.GENERIC (see rcm1d on this page). Thus, no need to recompute a 1D model but it will be needed if you want to model another Hot Jupiter.
Navigate to the 'makestart' folder, located at
<syntaxhighlight lang="bash">
dynamico-giant/hot_jupiter/makestart/
</syntaxhighlight>
To generate the initial conditions for the 3D run, we're gonna start the model using the temperature profile from the 1D run. to do that, you will find a "job_mpi" script. Open it, and adapt it to your cluster and launch the job. This job is using 20 procs, and it runs 5 days of simulations.
If everything goes well, you should see few netcdf files appear. The important ones are '''start_icosa0.nc''', '''startfi0.nc''' and '''Xhistins.nc'''.
If you see these files, you're all set to launch a real simulation !

'''5th step''': Go back to ''hot_jupiter'' folder. There are a bunch of script to launch your simulation. Take a look at the ''astro_fat_mpi'' script, and adapt it to your cluster. Then you can launch your simulation by doing
<syntaxhighlight lang="bash">
./run_astro_fat
</syntaxhighlight>
This will start the simulation, using 90 procs. In the same folder, check if the icosa_lmdz.out file is created. This is the logfile of the simulation, while it is running. You can check there that everything is going well.

'''Important side note''': When using the ''run_astro_fat'' script to run a simulation, it will run a chained simulation, restarting the simulation from the previous state after 100 days of simulations and generating ''Xhistins.nc'' files. This is your results file, where you will find all the variables that controls your atmosphere (temperature field, wind fields, etc..).

Good luck and enjoy the generic PCM Dynamico for Hot Jupiter !

== 3D LES setup ==

=== Proxima b with LES ===

TBD by Maxence

== 1D setup ==

=== rcm1d test case ===

Our 1-D forward model

TBD by Gwenael ? (you can have a look at the Generic GCM User Manual for inspiration)

=== kcm1d test case ===

Our 1-D inverse model

TBD by Guillaume or Martin

Tool Box

2022-05-11T09:41:54Z

LTeinturier: /* mass stream function */

== Pre-processing Tools ==
=== newstart ===
=== start2archive ===
=== visualization tools ===
=== other third party scripts and tools ===

TO BE COMPLETED

== Post-processing tools ==
=== zrecast ===

TO BE COMPLETED

=== mass stream function ===

The mass stream function (and the total angular momentum) can be computed from a diagfi.nc or a stats.nc, using the '''streamfunction.F90''' script. The script is located at

<syntaxhighlight lang="bash">
trunk/LMDZ.GENERIC/utilities
</syntaxhighlight>

To compile the script, open the ''compile'' file in the same directory and do the following:
* Replace "pgf90" with your favorite fortran compiler
* replace "/distrib/local/netcdf/pgi_7.1-6_32/lib" with the lib address and directory that contains your NetCDF library (file ''libnetcdf.a'').
* Replace "/distrib/local/netcdf/pgi_7.1-6_32/include" with the address of the directory that contains the NetCDF include file (''netcdf.inc'').
* You can mess with the compiling options but it is not mandatory.

Once the script is compiled, copy it in the same directory as your '''.nc''' file and run

<syntaxhighlight lang="bash">
./streamfunction.e
</syntaxhighlight>
The script will ask you for the name of your '''.nc''' file, and will run and produce a new '''nameofyourfile_stream.nc''' file.

'''Be careful''' : In this new file, all fields are temporally and zonally averaged.

If you want to use '''python''' instead of '''fortran''', you can take a look at this [https://github.com/aymeric-spiga/dynanalysis repo]. It hosts a tool to perform dynamical analysis of GCM simulations (and therefore, it computes the mass stream function and a lot of other stuff), but it is tailored for Dynamico only. This repo also takes care of recasting (it does the job of both ''zrecast.F90'' and ''streamfunction.F90'')

== Continuing Simulations ==

=== manually ===

At the end of a simulation, the model generates restart files (files 'restart.nc' and 'restartfi.nc') which contain the final state of the model. The 'restart.nc' and 'restartfi.nc' files have the same format as the 'start.nc' and 'startfi.nc' files, respectively.

These files can in fact be used as initial states to continue the simulation, using the following renaming command lines:

<syntaxhighlight lang="bash">
mv restart.nc start.nc
mv restartfi.nc startfi.nc
</syntaxhighlight>

Running a simulation with these start files will in fact resume the simulation from where the previous run ended.

=== with bash scripts ===

We have set up very simple bash scripts to automatize the launching of chain simulations. Here is an example of bash script that does the job:

<syntaxhighlight lang="bash">
#!/bin/bash
###########################################################################
# Script to perform several chained LMD Mars GCM simulations
# SET HERE the maximum total number of simulations

nummax=100

###########################################################################

echo "---------------------------------------------------------"
echo "STARTING LOOP RUN"
echo "---------------------------------------------------------"

dir=`pwd`
machine=`hostname`
address=`whoami`

# Look for file "num_run" which should contain
# the value of the previously computed season
# (defaults to 0 if file "num_run" does not exist)
if [[ -r num_run ]] ; then
echo "found file num_run"
numold=`cat num_run`
else
numold=0
fi
echo "numold is set to" ${numold}

# Set value of current season
(( numnew = ${numold} + 1 ))
echo "numnew is set to" ${numnew}

# Look for initialization data files (exit if none found)
if [[ ( -r start${numold}.nc && -r startfi${numold}.nc ) ]] ; then
\cp -f start${numold}.nc start.nc
\cp -f startfi${numold}.nc startfi.nc
else
if (( ${numold} == 99999 )) ; then
echo "No run because previous run crashed ! (99999 in num_run)"
exit
else
echo "Where is file start"${numold}".nc??"
exit
fi
fi

# Run GCM -- THIS LINE NEEDS TO BE MODIFIED WITH THE CORRECT GCM EXECUTION COMMAND
mpirun -np 8 gcm_64x48x26_phystd_para.e < diagfi.def > lrun${numnew}

# Check if run ended normaly and copy datafiles
if [[ ( -r restartfi.nc && -r restart.nc ) ]] ; then
echo "Run seems to have ended normaly"

\mv -f restart.nc start${numnew}.nc
\mv -f restartfi.nc startfi${numnew}.nc

else
if [[ -r num_run ]] ; then
\mv -f num_run num_run.crash
else
echo "No file num_run to build num_run.crash from !!"
# Impose a default value of 0 for num_run
echo 0 > num_run.crash
fi
echo 99999 > num_run
exit
fi

# Copy other datafiles that may have been generated
if [[ -r diagfi.nc ]] ; then
\mv -f diagfi.nc diagfi${numnew}.nc
fi
if [[ -r diagsoil.nc ]] ; then
\mv -f diagsoil.nc diagsoil${numnew}.nc
fi
if [[ -r stats.nc ]] ; then
\mv -f stats.nc stats${numnew}.nc
fi
if [[ -f profiles.dat ]] ; then
\mv -f profiles.dat profiles${numnew}.dat
\mv -f profiles.hdr profiles${numnew}.hdr
fi

# Prepare things for upcoming runs by writing
# value of computed season in file num_run
echo ${numnew} > num_run

# If we are over nummax : stop
if (( $numnew + 1 > $nummax )) ; then
exit
else
\cp -f run_gnome exe_mars
./exe_mars
fi

</syntaxhighlight>

'''Summary of what this bash script does''':

* It reads the file 'num_run' which contains the step of the simulation.
If num_run is
<syntaxhighlight lang="bash">
5
</syntaxhighlight>
then the script expects to read start5.nc and startfi5.nc.
* It modifies start5.nc and startfi5.nc into start.nc and startfi.nc, respectively.
* It runs the GCM.
* It modifies restart.nc and restartfi.nc into start6.nc and startfi6.nc
* It rewrite num_run as follows:
<syntaxhighlight lang="bash">
6
</syntaxhighlight>
* It restarts the loop until num_run reaches the value (defined in nummax):
<syntaxhighlight lang="bash">
100
</syntaxhighlight>

== Visualization software ==
=== Panoply ===
Panoply is a user-friendly tool for viewing raw NetCDF data, available here: https://www.giss.nasa.gov/tools/panoply/

[[File:Example panoply.png|thumb|Screenshot of panoply showing here LMD Generic results for the exoplanet TRAPPIST-1e]]

=== ncview ===
=== paraview ===
=== planetoplot ===
=== python scripts ===

TO BE COMPLETED

Dissipation

2022-05-11T09:08:34Z

LTeinturier: /* How to change it in the model */

== Description==
In the LMD 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 disspation 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 = 2, n = 1,and n = 2.

== How to change it in the model ==

In practise, the values of $$n$$ and $$\tau_{diss}$$ are prescribed in the [https://lmdz-forge.lmd.jussieu.fr/mediawiki/Planets/index.php/The_run.def_Input_File run.def] with the keys:
*nitergdiv
*nitergrot
*niterh
for the values of $$n$$ on each field, and the associated $$\tau$$:
*tetagdiv
*tetagrot
*tetatemp

In [https://lmdz-forge.lmd.jussieu.fr/mediawiki/Planets/index.php/The_run.def_Input_File run.def], there is also a key ''idissip'' which is the frequency (in dynamical steps) at which to apply the dissipation.

== 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$$.
* Optimal values for the dissipation timescales depends on the resolution of the horizontal grid. The higher the resolution, the more dissipation we need.

Dissipation

2022-05-11T09:04:16Z

LTeinturier: /* Description */

== Description==
In the LMD 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 disspation 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 = 2, n = 1,and n = 2.

== How to change it in the model ==

In practise, the values of $$n$$ and $$\tau_{diss}$$ are prescribed in the ''run.def'' with the keys:
*nitergdiv
*nitergrot
*niterh
for the values of $$n$$ on each field, and the associated $$\tau$$:
*tetagdiv
*tetagrot
*tetatemp

In ''run.def'', there is also a key ''idissip'' which is the frequency (in dynamical steps) at which to apply the dissipation.
== 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$$.
* Optimal values for the dissipation timescales depends on the resolution of the horizontal grid. The higher the resolution, the more dissipation we need.

Dissipation

2022-05-11T08:44:44Z

LTeinturier:

== Description==
In the LMD 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 $$\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 = 2, n = 1,and n = 2.
== How to change it in the model ==

In practise, the values of $$n$$ and $$\tau_{diss}$$ are prescribed in the ''run.def'' with the keys:
*nitergdiv
*nitergrot
*niterh
for the values of $$n$$ on each field, and the associated $$\tau$$:
*tetagdiv
*tetagrot
*tetatemp

In ''run.def'', there is also a key ''idissip'' which is the frequency (in dynamical steps) at which to apply the dissipation.
== 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$$.
* Optimal values for the dissipation timescales depends on the resolution of the horizontal grid. The higher the resolution, the more dissipation we need.

Dissipation

2022-05-11T08:37:19Z

LTeinturier: Created page with "In the LMD grid point model, nonlinear interactions between explicitly resolved scales and subgrid-scale processes are parameterized by applying a scale-selective horizontal d..."

In the LMD grid point model, nonlinear interactions between explicitly resolved scales
and subgrid-scale processes are parameterized by applying a scale-selective horizontal dis-
sipation 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 $$\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 = 2, n = 1,and n = 2.

In practise, the values of $$n$$ and $$\tau_{diss}$$ are prescribed in the ''run.def'' with the keys:
*nitergdiv
*nitergrot
*niterh
for the values of $$n$$ on each field, and the associated $$\tau$$:
*tetagdiv
*tetagrot
*tetatemp

Tool Box

2022-05-11T08:14:24Z

LTeinturier: /* mass stream function */

Tool Box

2022-05-11T08:14:12Z

LTeinturier: /* mass stream function */

Tool Box

2022-05-11T08:13:39Z

LTeinturier: /* mass stream function */