View source for PEM (Planetary Evolution Model)

== Overview of the PEM (Planetary Evolution Model) ==

The objective of the “Planetary Evolution Model” is to simulate and understand the past environments of a planet and it's long term evolution. It's also usefull to explain the observed landforms etc.

[[File:PEM.png|thumb|PEM schematic view ; credit: Francois Forget]]

To accurately simulate the climate and the fate of volatiles for thousands to millions of years we must couple
physical processes with very different timescale, ranging from clouds microphysics and atmospheric
dynamics (represented in the GCM) to the evolution of lakes, glacier accumulation, and subsurface ice
evolution.

Given the diversity and the complexity of the Martian paleoclimates, the model use an “asynchronous coupling” between
the slow ice and water reservoirs models and the GCM. In practice the Planetary evolution model
use a horizontal grid identical to that of the PCM, and include the same representation of the microclimate on slopes. 
In our case, we run the Planetary Evolution Model with an adaptative time step,
depending upon the dynamics of the modeled system (smaller timesteps must first be used so that the
different volatile reservoirs reach a quasi-equilibrium, then the timestep will depends on the evolution of the
forcing, which is slow in the case of obliquity, for instance) . At each timestep, the inputs from the
atmosphere (e.g. mean precipitation, sublimation and evaporation, temperatures, dust deposition) is obtained through a multi-annual run of the Global Climate model using the outcome of the Mars Evolution Model as initial state.

== The two modes of simulations ==

The PEM can works in two modes:
* M1 : Search of the steady state given an external forcing and an initial situation
* M2 : Realistic following of orbital forcing

For M1, we always apply the same external forcing (the GCM is always called at the same obliquity for example). The PEM will run to find a steady state. Be carefull that it can oscillate around the steady state so all the stopping criterions needs to be choosen wisely.
For M1, be also carefull that the uniquity of the steady state is not proven for this model. This means that for a given external forcing it is possible that the steady state is not unique and could depend on the initial condition and other parameters.

For M2, we follow a predefined orbit. It is read in the file ob_ex_lsp.asc. Stopping criterion specific to orbital change can be specified in run_pem.def. The steady state is not necessary reached for all the orbit configuration as it could take more time to reach it than the time spent in this orbital configuration.

== Stopping criterion ==

For both modes, the PEM will run and extrapolate tendencies 'as long as it is meaningfull'. This means that for the different physical processes modelled, a criterion is computed. This criterion express how much the system has changed. The PEM compares it to an acceptance rate and decides to stop or continue.
The criterions are the following: 
* '''Water ice ''' : Ratio of the surface of water ice that is sublimating at the time step t over the surface of water ice that is sublimating from the GCM run. 
* '''CO2_ice ''' : Ratio of the surface of co2 ice that is sublimating at the time step t over the surface of co2 ice that is sublimating from the GCM run. 
* '''Pressure''' : Ratio of the planet averaged surface pressure at iteration t over the planet averaged surface pressure at the end of the last GCM run.
* '''Orbital parameter''' : Multiple possibility; Percentage of change of a each orbital parameter, maximum absolute value of change, maximum number of years...

== Overview of a PEM run ==

In either mode, a run works like this:
* Run 2 years of GCM: From this, we extract data such as daily average pressure, daily averaged vmr of atmospheric components, daily minimum of ice etc... You thus need to use '''XIOS''' 
* The GCM restart file are renamed to be read by the PEM. They will be modified by the PEM to be able to start a new GCM simulation.
* If we are in a M2 configuration mode, a beginning year needs to be specified in the run_PEM.def, and the file ob_ex_lsp.asc needs to be present.
* The PEM is launched.
* GCM tendencies are computed from the data of the first item. They are applied to the start files and can be adapted to the situation (change of surface pressure etc...). All the physical processes also take place.
* The stopping criterion are checked and it either stop the time iteration or apply again the tendencies and run the physical processes
* GCM like start files are written.
* GCM is run again for 2 years and it loops again from the beginning.

== Physical processes modelled ==

The PEM is still in development, so far, the physical processes modeled are quite specific to Mars at low obliquity.

'''Current processes that have been validated:'''

✅ Ice accumulation 

✅ CO2 / H_2O adsorption 

✅ Soil temperature from the surface to ten's of kms.

✅  Glacier flow for CO2
 

'''Current processes that are in validation:'''

∼ Subsurface Ground Ice at equilibrium 

∼ Orbital Forcing (RV & LL)

'''Future works:'''

❌ Exchange between subsurface reservoirs and atmosphere / Formation of a lag deposit(will be done by Simon Neviere, Futur intern of Lucas)

❌ Dynamical Ice Table

❌ Hydrology: lake, river, ocean, etc. (Alexandre Gauvain)

❌ Photochemistry

== PEM inputs and outputs ==

Inputs:
* '''start_evol.nc''' A Dynamical GCM start file
* '''startfi_evol.nc''' A physic GCM start file
* '''startfi_PEM.nc''' A start file specific to the PEM (optional)
* '''run.def''' and '''run_PEM.def''' txt files specifying run parameters such as the initial year, the acceptance rate for each criterion etc...
* '''z2sig.def''' A txt file specifying the vertical grid
* '''tracer.def''' A txt file specifyng which tracers we are modelling.
* '''ob_ex_lsp.asc''' A txt files specifying orbit parameters (optional - only needed for M2)
* '''data_GCM_Y1.nc''' A netcdf file containing information about the first typical year run by the GCM
* '''data_GCM_Y2.nc''' A netcdf file containing information about the second typical year run by the GCM

Outputs:
* '''restart_evol.nc''' A Dynamical GCM start file to start the next GCM iteration run
* '''restartfi_evol.nc''' A physic GCM start file to start the next GCM iteration run
* '''restartfi_PEM.nc''' A start file specific to the PEM to start the next PCM run
* '''diagfi_PEM.nc''' (optional) A netcdf file containing physical information about the PCM run

== Naming convention of variables ==

Some variables comes from the GCM run. There are 3 cases:
* '''Constant variables''' : Some variable are not changed during the PEM run, they are named var_cGCM. ex:???
* '''Adapted variables''' : Some variables are adapted to the current state and recomputed through the PEM run, they are named var_aGCM. ex:???
* '''Tendencies from GCM''' : Some variables are tendencies obtained from the GCM run, they are named var_tGCM. ex:???

Some variables coming from the GCM can be shaped as the dynamical grid but needs to be used in the physical grid. The change is done via the subroutine grid_dyn_phys. In this routine, the name is changed:
* '''Dynamical to physical grid''' the variable is renamed from var to var_phys

Variables that exists in the GCM but that here belong to the PEM are named:
* '''Exist in the GCM but belong to the PEM ''' var_PEM

Variables that are specific to the PEM are not named following a special convention.

Index name for loops through:
* '''physical grid''' : ig (max ngrid)
* '''tracer''' : iq (max nq)
* '''subslopes''' : islope (max nslope)

TO DO:
* Rename variables
* Tendencies from the GCM : qsurf_GCM_tend

* Remove _slope

*ini_qsurf

*min_qsurf

== Choices of parameters (run_PEM.def) ==

Parameters of the PEM run can be specified in run_PEM.def. This file is included at the end of the run.def. For now, these are the possible parameters:

* evol_orbit_pem: Boolean. Do you want to follow an orbital forcing predefined (read in ob_ex_lsp.asc for example)? (default=false)
* year_bp_ini: Integer. Number of (Earth) year before present to start the pem run if evol_orbit_pem=.true. , default=0
* Max_iter_pem: Integer. Maximal number of iteration if none of the stopping criterion is reached and if evol_orbit_pem=.false., default=99999999
* dt_pem: Integer. Time step of the PEM in year, default=1
* alpha_criterion: Real. Acceptance rate of sublimating ice surface change, default=0.2
* soil_pem: Boolean. Do you want to run with subsurface physical processes in the PEM? default=.true.

== How to compile the model ==

We compile the model the same way as the GCM (see section Quick Install and Run).
So far it can be compiled only in sequential mode, we qlso have to specify the grid, the planet, the output utilitary etc... 
An example of compilation line command:

The Fortran code of the PEM is in the following directory : trunk/LMDZ.COMMON/libf/evolution/
The main program is in the file pem.F90

For the PEM, the output of the GCM are usually done using XIOS. This allow the possibility to do averages and minimum of variable easily. Example of .xml files used for PEM run can be found in : trunk/???

== How to run the model ==

An example of script to run the PEM in a cluster can be found in the following directory : trunk/???