diff --git a/examples/reactmicp/unsaturated/drying.cpp b/examples/reactmicp/unsaturated/drying.cpp
index 457e8a8..a9abacc 100644
--- a/examples/reactmicp/unsaturated/drying.cpp
+++ b/examples/reactmicp/unsaturated/drying.cpp
@@ -1,578 +1,578 @@
 /* =============================================================================
 
  Copyright (c) 2014 - 2016
  F. Georget <fabieng@princeton.edu> Princeton University
  All rights reserved.
 
 
 Redistribution and use in source and binary forms, with or without modification,
 are permitted provided that the following conditions are met:
 
     1. Redistributions of source code must retain the above copyright notice,
     this list of conditions and the following disclaimer.
 
     2. Redistributions in binary form must reproduce the above copyright notice,
     this list of conditions and the following disclaimer in the documentation
     and/or other materials provided with the distribution.
 
     3. Neither the name of the copyright holder nor the names of its
     contributors may be used to endorse or promote products derived from this
     software without specific prior written permission.
 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
 ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
 ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. *
 
 ============================================================================= */
 
 
 // ====================================================== //
 //                                                        //
 //                     Drying example                     //
 //                     ==============                     //
 //                                                        //
 // ====================================================== //
 //
 //
 // This is a simple example to show how the unsaturated
 // system of ReactMiCP is able to solve a complex drying
 // equation.
 //
 // It implements a special chemistry stagger (DryingStagger)
 // which only computes the vapor pressure.
 
 
 // Macro to define the problem
 // ===========================
 
 #define NB_NODES 50
 #define DX 0.05*1e-2                 // cm
 #define CROSS_SECTION 5.0*1e-4      // cm^2
 
 #define GAMMA_WATER_GLASS 0.1  // N/m
 #define GAMMA_WATER_AIR 71.97e-3  // N/m = kg/s^2
 #define R_SMALL_BEAD 0.015*1e-2     // cm
 #define R_BIG_BEAD 0.04*1e-2        // cm
 
 #define M_v 18.015e-3          // kg/mol
 
 // Includes
 // ========
 
 #include "specmicp_database/database.hpp"
 
 #include "reactmicp/systems/unsaturated/variables.hpp"
 #include "reactmicp/systems/unsaturated/variables_sub.hpp"
 #include "reactmicp/systems/unsaturated/transport_stagger.hpp"
 
 #include "reactmicp/solver/staggers_base/chemistry_stagger_base.hpp"
 #include "reactmicp/solver/staggers_base/upscaling_stagger_base.hpp"
 #include "reactmicp/solver/staggers_base/stagger_structs.hpp"
 
 #include "reactmicp/solver/runner.hpp"
 
 #include "dfpm/meshes/generic_mesh1d.hpp"
 
 #include "specmicp_common/physics/laws.hpp"
 #include "specmicp_common/physics/units.hpp"
 
 #include "specmicp_common/log.hpp"
 #include "specmicp_common/io/csv_formatter.hpp"
 
 #include <functional>
 #include <iostream>
 
 // Namespace
 // =========
 using namespace specmicp;
 using namespace specmicp::reactmicp;
 using namespace specmicp::reactmicp::systems;
 using namespace specmicp::reactmicp::systems::unsaturated;
 
 // Declarations
 // ============
 
 // Constants
 // ==========
 static const scalar_t rho_l_si = constants::water_density_25;
 static const scalar_t rho_l = 1e-6*constants::water_density_25;
 static const scalar_t cw_tilde_si = rho_l_si / M_v;
 static const scalar_t cw_tilde = 1e-6*cw_tilde_si;
 
 int main();
 
 database::RawDatabasePtr get_database();
 mesh::Mesh1DPtr get_mesh();
 
 scalar_t leverett_function(scalar_t s_eff);
 
 scalar_t big_bead_cap_pressure(scalar_t saturation);
 scalar_t small_bead_cap_pressure(scalar_t saturation);
 scalar_t vapor_pressure(scalar_t pc);
 
 // Drying Stagger
 // ==============
 //
 // This stagger computes the exchange between water vapor and liquid water
 class DryingStagger: public solver::ChemistryStaggerBase
 {
 public:
 
     DryingStagger() {}
 
     virtual void initialize(solver::VariablesBase * const var) override {}
 
      //! \brief Initialize the stagger at the beginning of an iteration
     //!
     //! This is where the predictor can be saved, the first trivial iteration done, ...
     //!
     //! \param dt the duration of the timestep
     //! \param var a shared_ptr to the variables
     virtual void initialize_timestep(
             scalar_t dt,
             solver::VariablesBase * const var
             ) override {
         m_dt = dt;
     }
 
     //! \brief Solve the equation for the timestep
     //!
     //! \param var a shared_ptr to the variables
     virtual solver::StaggerReturnCode restart_timestep(
             solver::VariablesBase * const var) override;
 
     void compute_one_node(index_t node, UnsaturatedVariables * const vars);
 
 private:
     scalar_t m_dt {-1.0};
 };
 
 // Upscaling Stagger
 // =================
 //
 // Define how we compute capillary pressure, vapor pressure, ...
 class UpscalingDryingStagger: public solver::UpscalingStaggerBase
 {
 public:
     UpscalingDryingStagger(std::vector<bool>& is_big_bead):
         m_is_big_bead(is_big_bead)
     {}
 
     //! \brief Initialize the stagger at the beginning of the computation
     //!
     //! \param var a shared_ptr to the variables
     virtual void initialize(solver::VariablesBase * const var) override;
 
     //! \brief Initialize the stagger at the beginning of an iteration
     //!
     //! This is where the predictor can be saved, the first trivial iteration done, ...
     //!
     //! \param dt the duration of the timestep
     //! \param var a shared_ptr to the variables
     virtual void initialize_timestep(
             scalar_t dt,
             solver::VariablesBase * const var
             ) override {}
 
     //! \brief Solve the equation for the timestep
     //!
     //! \param var a shared_ptr to the variables
     virtual solver::StaggerReturnCode restart_timestep(
             solver::VariablesBase * const var) override
     {
         return solver::StaggerReturnCode::ResidualMinimized;
     }
 
 
     scalar_t capillary_pressure_model(index_t node, scalar_t saturation);
 
     scalar_t vapor_pressure_model(index_t node, scalar_t saturation);
 
     scalar_t relative_liquid_permeability_model(index_t node, scalar_t saturation);
 
     scalar_t relative_liquid_diffusivity_model(index_t node, scalar_t saturation);
 
     scalar_t relative_gas_diffusivity_model(index_t node, scalar_t saturation);
 
 private:
     std::vector<bool> m_is_big_bead;
 };
 
 
 
 // Implementation
 // ==============
 
 scalar_t UpscalingDryingStagger::capillary_pressure_model(index_t node, scalar_t saturation)
 {
     scalar_t pc = 0.0;
     if (m_is_big_bead[node])
     {
         pc = big_bead_cap_pressure(saturation);
     }
     else
         pc = small_bead_cap_pressure(saturation);
     return pc;
 }
 
 scalar_t UpscalingDryingStagger::vapor_pressure_model(index_t node, scalar_t saturation)
 {
     return vapor_pressure(capillary_pressure_model(node, saturation));
 }
 
 
 scalar_t UpscalingDryingStagger::relative_liquid_permeability_model(index_t node, scalar_t saturation)
 {
     return std::pow(saturation, 3);
 }
 
 scalar_t UpscalingDryingStagger::relative_liquid_diffusivity_model(index_t node, scalar_t saturation)
 {
     return saturation;
 }
 
 scalar_t UpscalingDryingStagger::relative_gas_diffusivity_model(index_t node, scalar_t saturation)
 {
     return (1.0 - saturation);
 }
 
 void UpscalingDryingStagger::initialize(solver::VariablesBase * const var)
 {
     UnsaturatedVariables * const true_vars = static_cast<UnsaturatedVariables * const>(var);
 
     true_vars->set_capillary_pressure_model(
                 [this](index_t node, scalar_t sat) -> scalar_t {
         return this->capillary_pressure_model(node, sat);});
 
     true_vars->set_vapor_pressure_model(std::bind(std::mem_fn(
                             &UpscalingDryingStagger::vapor_pressure_model
                             ), this, std::placeholders::_1, std::placeholders::_2));
 
     true_vars->set_relative_liquid_permeability_model(std::bind(std::mem_fn(
                             &UpscalingDryingStagger::relative_liquid_permeability_model
                             ), this, std::placeholders::_1, std::placeholders::_2));
 
     true_vars->set_relative_liquid_diffusivity_model(std::bind(std::mem_fn(
                             &UpscalingDryingStagger::relative_liquid_diffusivity_model
                             ), this, std::placeholders::_1, std::placeholders::_2));
 
     true_vars->set_relative_gas_diffusivity_model(std::bind(std::mem_fn(
                             &UpscalingDryingStagger::relative_gas_diffusivity_model
                             ), this, std::placeholders::_1, std::placeholders::_2));
 
     MainVariable& vapor_pressure = true_vars->get_pressure_main_variables(0);
     MainVariable& saturation = true_vars->get_liquid_saturation();
 
     SecondaryTransientVariable& porosity   = true_vars->get_porosity();
     SecondaryVariable& liquid_diffusivity  = true_vars->get_liquid_diffusivity();
     SecondaryVariable& liquid_permeability = true_vars->get_liquid_permeability();
     SecondaryVariable& resistance_gas_diffusivity = true_vars->get_resistance_gas_diffusivity();
 
     true_vars->get_binary_gas_diffusivity(0) =  0.282*1e-4;
 
     porosity(0) = 1;
     resistance_gas_diffusivity(0) = 1.0;
 
 
     for (index_t node=1; node<true_vars->get_mesh()->nb_nodes(); ++node)
     {
         vapor_pressure(node) = vapor_pressure_model(node, saturation(node));
 
         if (m_is_big_bead[node]) {
             porosity(node) = 0.371;
             liquid_diffusivity(node) = 0;
             liquid_permeability(node) = 3.52e-11;
             resistance_gas_diffusivity(node) = 2*0.371/(3-0.371);
         } else {
             porosity(node) = 0.387;
             liquid_diffusivity(node) = 0;
             liquid_permeability(node) = 8.41e-12;
             resistance_gas_diffusivity(node) = 2*0.387/(3-0.387);
         }
     }
 
     porosity.predictor = porosity.variable;
 
     true_vars->set_relative_variables();
 }
 
 
 
 solver::StaggerReturnCode DryingStagger::restart_timestep(
         solver::VariablesBase * const var)
 {
     UnsaturatedVariables * const true_vars = static_cast<UnsaturatedVariables * const>(var);
 
     for (index_t node=1; node<true_vars->get_mesh()->nb_nodes(); ++node) {
         compute_one_node(node, true_vars);
     }
 
     return solver::StaggerReturnCode::ResidualMinimized;
 }
 
 void DryingStagger::compute_one_node(
         index_t node,
         UnsaturatedVariables * const vars)
 {
     MainVariable& satvars = vars->get_liquid_saturation();
     MainVariable& presvars = vars->get_pressure_main_variables(0);
 
 
     scalar_t rt = vars->get_rt();
     scalar_t saturation = satvars(node);
     scalar_t pressure = presvars(node);
     scalar_t porosity = vars->get_porosity()(node);
     auto vapor_pressure_f = vars->get_vapor_pressure_model();
 
     const scalar_t tot_conc = cw_tilde_si*saturation + (1.0-saturation)*pressure/rt;
 
     pressure = vapor_pressure_f(node, saturation);
     scalar_t cw_hat = pressure/rt;
     saturation = (tot_conc - cw_hat) / (cw_tilde_si-cw_hat);
 
     pressure = vapor_pressure_f(node, saturation);
     cw_hat = pressure/rt;
 
     scalar_t res = tot_conc - cw_tilde_si*saturation - (1.0-saturation)*cw_hat;
 
     while (std::abs(res)/tot_conc > 1e-12) {
         saturation =  (tot_conc - cw_hat) / (cw_tilde_si-cw_hat);
         pressure = vapor_pressure_f(node, saturation);
         cw_hat = pressure/rt;
 
         res = tot_conc - cw_tilde_si*saturation - (1.0-saturation)*cw_hat;
     }
 
     satvars(node) = saturation;
     satvars.velocity(node) = (saturation - satvars.predictor(node))/m_dt;
 
     presvars(node) = pressure;
     presvars.velocity(node) = (pressure - presvars.predictor(node))/m_dt;
     presvars.chemistry_rate(node) =  - porosity/(rt*m_dt) * (
                 (pressure*(1.0-saturation)) -
                 (presvars.predictor(node)*(1.0-satvars.predictor(node)))
                 )
             + presvars.transport_fluxes(node);
 }
 
 
 database::RawDatabasePtr get_database()
 {
 
     specmicp::database::Database thedatabase(CEMDATA_PATH);
     std::vector<std::string> to_keep = {"H2O", "H[+]"};
     thedatabase.keep_only_components(to_keep);
     thedatabase.remove_gas_phases();
     database::RawDatabasePtr raw_data = thedatabase.get_database();
     raw_data->freeze_db();
     return raw_data;
 }
 
 mesh::Mesh1DPtr get_mesh()
 {
     mesh::Uniform1DMeshGeometry geom;
     geom.dx = DX;
     geom.nb_nodes = NB_NODES;
     geom.section = CROSS_SECTION;
     return mesh::uniform_mesh1d(geom);
 }
 
 scalar_t big_bead_cap_pressure(scalar_t saturation) {
     return GAMMA_WATER_AIR/std::sqrt(3.52e-11/0.371)*leverett_function(saturation);
 }
 
 scalar_t small_bead_cap_pressure(scalar_t saturation) {
     return GAMMA_WATER_AIR/std::sqrt(8.41e-12/0.387)*leverett_function(saturation);
 }
 
 scalar_t vapor_pressure(scalar_t pc) {
     return 3e3*std::exp(-M_v * pc / (rho_l_si * (8.314*(25+273.15))));
 }
 
 
 scalar_t leverett_function(scalar_t s_eff) {
     return 0.325*std::pow((1.0/s_eff - 1), 0.217);
 }
 
 
 class Printer
 {
 public:
     Printer(UnsaturatedVariables* vars, std::size_t reserve_size);
     
     
     void register_timestep(scalar_t time, solver::VariablesBasePtr base_vars);
     
     void print(const std::string& name);
     
 private:
     mesh::Mesh1DPtr m_mesh;
 
     std::vector<std::vector<scalar_t>> m_pressure;
     std::vector<std::vector<scalar_t>> m_saturation;
     std::vector<scalar_t> m_timestep;
 };
 
 Printer::Printer(UnsaturatedVariables *vars, std::size_t reserve_size):
     m_mesh(vars->get_mesh()),
     m_pressure(m_mesh->nb_nodes()),
     m_saturation(m_mesh->nb_nodes())
 {
     m_timestep.reserve(reserve_size);
     for (auto node: m_mesh->range_nodes())
     {
         m_pressure[node].reserve(reserve_size);
         m_saturation[node].reserve(reserve_size);
     }
 }
 
 void Printer::register_timestep(scalar_t time, solver::VariablesBasePtr base_vars)
 {
     UnsaturatedVariables* vars = static_cast<UnsaturatedVariables*>(base_vars.get());
 
     m_timestep.push_back(time);
     MainVariable& saturation = vars->get_liquid_saturation();
     MainVariable& pressure = vars->get_pressure_main_variables(0);
 
     for (auto node: m_mesh->range_nodes())
     {
         m_pressure[node].push_back(pressure(node));
         m_saturation[node].push_back(saturation(node));
     }
 }
 
 void Printer::print(const std::string& name)
 {
 
     io::CSVFile pressure_file(name+"_pressure.dat");
     io::CSVFile saturation_file(name+"_saturation.dat");
 
     pressure_file << "Node";
     saturation_file << "Node";
     for (auto& time: m_timestep)
     {
         pressure_file.separator();
         pressure_file << time;
         saturation_file.separator();
         saturation_file << time;
     }
     pressure_file.eol();
     saturation_file.eol();
 
     for (auto node: m_mesh->range_nodes())
     {
         pressure_file << m_mesh->get_position(node);
         saturation_file << m_mesh->get_position(node);
         for (std::size_t ind=0; ind<m_timestep.size(); ++ind)
         {
             pressure_file.separator();
             pressure_file << m_pressure[node][ind];
             saturation_file.separator();
             saturation_file << m_saturation[node][ind];
         }
         pressure_file.eol();
         saturation_file.eol();
     }
 
     pressure_file.flush();
     saturation_file.flush();
 }
 
 
 
 
 int main()
 {
     init_logger(&std::cerr, logger::Warning);
 
     mesh::Mesh1DPtr the_mesh = get_mesh();
     database::RawDatabasePtr raw_data = get_database();
     std::vector<bool> has_gas = {true, false, false};
 
     std::vector<bool> is_big_bead(the_mesh->nb_nodes(), true);
     //for (int node=25; node<NB_NODES; ++node)
     //    is_big_bead[node] = true;
 
     UnsaturatedVariablesPtr vars = std::make_shared<UnsaturatedVariables>(
-                the_mesh, raw_data, has_gas, units::LengthUnit::meter);
+                the_mesh, raw_data, has_gas, units::SI_units);
 
 
     auto bcs = BoundaryConditions::make(the_mesh->nb_nodes(), raw_data->nb_component());
     bcs->add_gas_node(0);
     bcs->fork_constraint(0, "gas_node");
 
 
     // variables initialisation
     MainVariable& saturation = vars->get_liquid_saturation();
     MainVariable& vapor_pressure = vars->get_pressure_main_variables(0);
     SecondaryTransientVariable& water_aqueous_concentration = vars->get_water_aqueous_concentration();
     water_aqueous_concentration.predictor = water_aqueous_concentration.variable;
 
     SecondaryVariable& rel_gas_d = vars->get_relative_gas_diffusivity();
     rel_gas_d(0) = 1.0;
 
     vapor_pressure(0) = 100;
 
     for (index_t node = 1; node < the_mesh->nb_nodes(); ++node) {
         saturation(node) = 0.9;
     }
 
     water_aqueous_concentration.set_constant(cw_tilde_si);
 
 
 
     std::shared_ptr<UpscalingDryingStagger> upscaling_stagger = std::make_shared<UpscalingDryingStagger>(is_big_bead);
     upscaling_stagger->initialize(vars.get());
 
     std::shared_ptr<DryingStagger> drying_stagger = std::make_shared<DryingStagger>();
     // init vapor pressure
     drying_stagger->initialize_timestep(1.0, vars.get());
     drying_stagger->restart_timestep(vars.get());
 
 
     std::shared_ptr<UnsaturatedTransportStagger> transport_stagger =
             std::make_shared<UnsaturatedTransportStagger>(vars, bcs);
 
     solver::ReactiveTransportSolver react_solver(transport_stagger, drying_stagger, upscaling_stagger);
     solver::ReactiveTransportOptions& opts = react_solver.get_options();
     opts.residuals_tolerance = 1e-2;
     opts.good_enough_tolerance = 0.99;
 
     vapor_pressure.chemistry_rate.setZero();
 
     //for (auto it=0; it<100; ++it)
     //{
 //    auto retcode = react_solver.solve_timestep(0.1, vars);
 //    std::cout << (int) retcode << std::endl;
     //}
 
     //std::cout << saturation.variable << std::endl;
     //std::cout << vapor_pressure.variable << std::endl;
     //std::cout << vars->get_gas_diffusivity().variable << std::endl;
     //std::cout << vars->get_relative_gas_diffusivity().variable << std::endl;
 
     int nb_hours = 9;
     Printer printer(vars.get(), 4*nb_hours);
     solver::SimulationInformation simul_info("drying", 900);
     solver::ReactiveTransportRunner runner(react_solver, 0.01, 10, simul_info);
 
     solver::output_f out_func = std::bind(&Printer::register_timestep,
                              &printer,
                              std::placeholders::_1, std::placeholders::_2);
 
     runner.set_output_policy(out_func);
 
     runner.run_until(nb_hours*3600, vars);
 
     printer.print("out_drying");
 
     std::cout << saturation.variable << std::endl;
     std::cout << vapor_pressure.variable << std::endl;
 }
diff --git a/src/bin/reactmicp/unsaturated.cpp b/src/bin/reactmicp/unsaturated.cpp
index f8d4068..1554e6d 100644
--- a/src/bin/reactmicp/unsaturated.cpp
+++ b/src/bin/reactmicp/unsaturated.cpp
@@ -1,857 +1,858 @@
 /* =============================================================================
 
  Copyright (c) 2014 - 2016
  F. Georget <fabieng@princeton.edu> Princeton University
  All rights reserved.
 
 
 Redistribution and use in source and binary forms, with or without modification,
 are permitted provided that the following conditions are met:
 
     1. Redistributions of source code must retain the above copyright notice,
     this list of conditions and the following disclaimer.
 
     2. Redistributions in binary form must reproduce the above copyright notice,
     this list of conditions and the following disclaimer in the documentation
     and/or other materials provided with the distribution.
 
     3. Neither the name of the copyright holder nor the names of its
     contributors may be used to endorse or promote products derived from this
     software without specific prior written permission.
 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
 ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
 ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. *
 
 ============================================================================= */
 
 //! \file bin/reactmicp/unsaturated.cpp
 //! \brief Unsaturated
 
 #include "specmicp_common/types.hpp"
 
 #include "specmicp_common/cli/parser.hpp"
 #include "specmicp_common/plugins/plugin_manager.hpp"
 #include "specmicp_common/io/yaml.hpp"
 #include "specmicp_common/io/safe_config.hpp"
 #include "specmicp_common/compat.hpp"
 #include "specmicp_common/log.hpp"
 #include "specmicp_common/filesystem.hpp"
 #include "specmicp_common/string_algorithms.hpp"
 #include "specmicp_common/openmp_support.hpp"
 
 #include "specmicp_common/physics/io/configuration.hpp"
 #include "specmicp_database/io/configuration.hpp"
 #include "dfpm/io/configuration.hpp"
 #include "specmicp_common/io/configuration.hpp"
 
 #include "specmicp/io/print.hpp"
 
 #include "reactmicp/reactmicp_unsaturated.hpp"
 #include "reactmicp/systems/unsaturated/variables.hpp"
 #include "dfpm/solver/parabolic_structs.hpp"
 
 #include "reactmicp/io/hdf5_unsaturated.hpp"
 
 #include "specmicp_common/config.h"
 
 #include <string>
 #include <memory>
 #include <stdexcept>
 #include <iostream>
 
 #include "specmicp_common/io/config_yaml_sections.h"
 
 #define PLUGIN_MODULE_UPSCALING "upscaling_factory"
 #define PLUGIN_MODULE_VARIABLES "variables_factory"
 
 using namespace specmicp;
 using namespace specmicp::reactmicp;
 using namespace specmicp::reactmicp::systems::unsaturated;
 
 const char* welcome =
 R"(===============================
 =    Reactmicp - unsaturated  =
 ===============================
 (C) copyright 2014-2016 F.Georget, Princeton University
 
 version: not good enough yet
 )";
 
 int main (int argc, char* argv[]);
 
 plugins::PluginManager& initialize_plugin_manager(
         io::YAMLConfigHandle& config
         );
 
 SimulationData initialize_simul_data(
         io::YAMLConfigHandle& safe_config,
         const std::vector<std::string>& db_dirs
         );
 
 std::shared_ptr<UnsaturatedVariables> initialize_variables(
         std::shared_ptr<database::DataContainer> raw_data,
         std::shared_ptr<mesh::Mesh1D> the_mesh,
         const units::UnitsSet& the_units,
         io::YAMLConfigHandle&& configuration
         );
 
 struct Staggers
 {
     std::shared_ptr<solver::UpscalingStaggerBase> upscaling;
     std::shared_ptr<solver::ChemistryStaggerBase>  chemistry;
     std::shared_ptr<solver::TransportStaggerBase>  transport;
 };
 
 Staggers initialize_staggers(
         SimulationData& simul_data,
         io::YAMLConfigHandle&& configuration
         );
 
 std::shared_ptr<solver::UpscalingStaggerBase>
 initialize_upscaling_stagger(
         SimulationData& simul_data,
         io::YAMLConfigHandle& conf_staggers
         );
 
 std::shared_ptr<solver::ChemistryStaggerBase>
 initialize_chemistry_stagger(
         SimulationData& simul_data,
         io::YAMLConfigHandle& conf_staggers
         );
 
 std::shared_ptr<solver::TransportStaggerBase>
 initialize_transport_stagger(
         SimulationData& simul_data,
         io::YAMLConfigHandle& conf_staggers
         );
 
 
 std::string initialize_output(
         SimulationData& simul_data,
         const std::string& work_dir,
         io::YAMLConfigHandle&& conf
         );
 
 void print_db_to_conf_log(database::DataContainer* raw_data);
 
 int main(int argc, char *argv[])
 {
     std::cout << welcome << std::endl;
 
     init_logger(&std::cout, logger::Warning);
 
     // CLI options
     cli::CommandLineParser parser;
     parser.add_option('i', "input",
                       cli::ValueType::string,
                       "input_file (yaml format)");
     parser.add_option('d', "database_dir",
                       std::string(""),
                       "directory where database are stored");
     parser.add_option('w', "workdir",
                       utils::get_current_directory(),
                       "working directory");
     parser.set_help_message("ReactMiCP Unsaturated System");
     const auto ret_parse = parser.parse(argc, argv);
     if (ret_parse > 0) {
         return EXIT_SUCCESS;
     }
 
 
     std::string input_filepath =
             parser.get_option<cli::ValueType::string>("input");
     std::string database_dir =
             parser.get_option<cli::ValueType::string>("database_dir");
     std::string working_dir =
             parser.get_option<cli::ValueType::string>("workdir");
 
     std::vector<std::string> db_dirs;
     const auto env = utils::get_env("SPECMICP_DATABASE_PATH");
     if (not env.empty()) {
         db_dirs = utils::split(env, ';');
     }
     if (not database_dir.empty()) {db_dirs.push_back(database_dir);}
     if (not working_dir.empty()) {db_dirs.push_back(working_dir);}
 
     // Configuration File
     // ------------------
     auto config = io::YAMLConfigFile::make(input_filepath);
 
     std::unique_ptr<std::ostream> out_conf_logger;
     if (config->has_section(SPC_CF_S_CONF_LOGS)) {
         out_conf_logger = io::configure_conf_log(
                     config->get_section(SPC_CF_S_CONF_LOGS),
                     working_dir
                     );
     }
     SPC_CONF_LOG << "Running ReactMiCP unsaturated driver";
 
     std::unique_ptr<std::ostream> out_logger;
     if (config->has_section(SPC_CF_S_LOGS)) {
         out_logger = io::configure_log(
                     config->get_section(SPC_CF_S_LOGS),
                     working_dir
                     );
     }
 
     initialize_plugin_manager(*(config.get()));
 
     // openmp initialization
     // ---------------------
     utils::initialize_parallel();
 
     // Simulation Data
     // ---------------
     SimulationData simul_data = initialize_simul_data(*config.get(), db_dirs);
 
     // Staggers
     // -------
     Staggers staggers = initialize_staggers(
                 simul_data,
                 config->get_section(SPC_CF_S_STAGGERS)
                 );
 
     // Check the variables
     // -------------------
 
     auto res = VariablesInterface(simul_data.variables).check_variables();
     if (res != VariablesValidity::good) {
         const auto msg = "Variables were not initialized correctly."
                          "Please see log for errors.";
         ERROR << msg;
         throw std::invalid_argument(msg);
 
     }
 
     // Solver
     // ------
 
     solver::ReactiveTransportSolver react_solver(
                 staggers.transport, staggers.chemistry, staggers.upscaling);
     io::configure_reactmicp_options(
                 react_solver.get_options(),
                 config->get_section(SPC_CF_S_REACTMICP)
                 );
     SPC_CONF_LOG << "\n\n++ ReactMiCP ++\n"
                  << SPC_CONF_LOG_SECTION
                  << "\n ReactMiCP options\n"
                  << SPC_CONF_LOG_HLINE;
     std::ostringstream msg;
     io::print_reactmicp_options(msg, react_solver.get_options());
     SPC_CONF_LOG << msg.str();
 
     // Initialize output
     std::unique_ptr<io::UnsaturatedHDF5Saver> saver;
 
     solver::SimulationInformation simul_info =
             io::configure_simulation_information(
                 config->get_section(SPC_CF_S_SIMULINFO));
     simul_info.working_dir = working_dir;
 
     solver::ReactiveTransportRunner runner(react_solver,
                                            1.0, 10.0, // dummy value, set later
                                            simul_info);
     io::configure_reactmicp_timestepper(
                 runner.get_timestepper_options(),
                 config->get_section(SPC_CF_S_TIMESTEPPER)
                 );
 
     //
     if (config->has_section(SPC_CF_S_REACTOUTPUT)) {
         auto filepath = initialize_output(
                     simul_data,
                     working_dir,
                     config->get_section(SPC_CF_S_REACTOUTPUT));
         saver = make_unique<io::UnsaturatedHDF5Saver>(
                     filepath, simul_data.variables, simul_data.units);
         auto* saver_ptr = saver.get();
         auto out_pol = [saver_ptr](
                        scalar_t timestep,
                 solver::VariablesBasePtr _) {
             saver_ptr->save_timestep(timestep);
             };
         saver->save_timestep(0.0);
         runner.set_output_policy(out_pol);
     }
     //  run info
     auto run_section = config->get_section(SPC_CF_S_RUN);
     auto run_until = run_section.get_required_attribute<scalar_t>(
                 SPC_CF_S_RUN_A_RUNUTNIL);
 
     // we destroy config to get information
     config.reset(nullptr);
 
     runner.run_until(run_until, simul_data.variables);
 
 
     return EXIT_SUCCESS;
 }
 
 // ==================== //
 //    Simulation data   //
 // ==================== //
 
 SimulationData initialize_simul_data(
         io::YAMLConfigHandle& safe_config,
         const std::vector<std::string>& db_dirs
         )
 {
     // Package the information
     SimulationData simul_data;
 
     // The units
     // ---------
     if (safe_config.has_section(SPC_CF_S_UNITS))
     {
         simul_data.units = io::configure_units(
                     safe_config.get_section(SPC_CF_S_UNITS));
     }
     SPC_CONF_LOG << SPC_CONF_LOG_HLINE
                  << "\nUnits :"
-                 << "\n  - length : " << io::to_string(simul_data.units.length)
-                 << "\n  - mass   : " << io::to_string(simul_data.units.mass)
+                 << "\n  - length   : " << io::to_string(simul_data.units.length)
+                 << "\n  - mass     : " << io::to_string(simul_data.units.mass)
+                 << "\n  - quantity : " << io::to_string(simul_data.units.quantity)
                  << "\n" << SPC_CONF_LOG_HLINE;
 
     // The database
     // ------------
     SPC_CONF_LOG << "Database initialization\n";
     simul_data.raw_data = io::configure_database(
                 safe_config.get_section(SPC_CF_S_DATABASE),
                 db_dirs
                 );
     SPC_CONF_LOG << SPC_CONF_LOG_HLINE;
     if (not simul_data.raw_data->is_valid()) {
         CRITICAL << "Invalid database from the configuration.";
         throw std::runtime_error("Database parsed from the configuration has"
                                  " been detected to be invalid. Abort.");
     }
     // The mesh
     // --------
     SPC_CONF_LOG << "Mesh initialization\n"
                  << SPC_CONF_LOG_SECTION;
     simul_data.mesh1d =
             io::configure_mesh(safe_config.get_section(SPC_CF_S_MESH));
     SPC_CONF_LOG << SPC_CONF_LOG_HLINE;
 
     // The variables
     // -------------
     SPC_CONF_LOG << "Variables initialization\n"
                  << SPC_CONF_LOG_SECTION;
     simul_data.variables = initialize_variables(
                                simul_data.raw_data,
                                simul_data.mesh1d,
                                simul_data.units,
                                safe_config.get_section(SPC_CF_S_REACTMICPVARIABLES)
                                );
 
     if (not simul_data.raw_data->is_valid()) {
         CRITICAL << "Invalid database after variable initialization.";
         throw std::runtime_error("Database has been detected to be invalid"
                                  " after the variables initialization. Abort.");
     }
     SPC_CONF_LOG << SPC_CONF_LOG_HLINE;
 
     // The database is supposed set after this point !
     simul_data.raw_data->freeze_db();
     print_db_to_conf_log(simul_data.raw_data.get());
 
     // The boundary conditions
     // -----------------------
     SPC_CONF_LOG << "Boundary conditions initialization\n"
                  << SPC_CONF_LOG_SECTION;
     simul_data.bcs =
             io::configure_unsaturated_boundary_conditions(
                 simul_data.mesh1d->nb_nodes(),
                 simul_data.raw_data.get(),
                 safe_config.get_section(SPC_CF_S_BOUNDARYCONDITIONS)
                 );
     SPC_CONF_LOG << SPC_CONF_LOG_HLINE;
 
     // Simulation data ready
     return simul_data;
 }
 
 
 std::shared_ptr<UnsaturatedVariables> initialize_variables(
         std::shared_ptr<database::DataContainer> raw_data,
         std::shared_ptr<mesh::Mesh1D> the_mesh,
         const units::UnitsSet& the_units,
         io::YAMLConfigHandle&& configuration
         )
 {
     auto type = configuration.get_optional_attribute<std::string>(
                 SPC_CF_S_REACTMICPVARIABLES_A_TYPE, SPC_CF_V_PLUGIN);
     if (type == SPC_CF_V_DEFAULT) type = SPC_CF_V_PLUGIN;
     if (type != SPC_CF_V_PLUGIN) {
         configuration.report_error(
                     io::YAMLConfigError::InvalidArgument,
                     "Only accepted value for attribute '"
                     SPC_CF_S_REACTMICPVARIABLES_A_TYPE
                     "' is '" SPC_CF_V_PLUGIN "'."
                     );
     }
     auto plugin_name = configuration.get_required_attribute<std::string>(
                 SPC_CF_A_PLUGIN_NAME
                 );
     SPC_CONF_LOG << "Plugin name : " << plugin_name;
     plugins::PluginManager& plugin_manager = plugins::get_plugin_manager();
     auto extra_plugin_file = configuration.get_optional_attribute<std::string>(
                                  SPC_CF_A_PLUGIN_FILE, "");
 
     if (extra_plugin_file != "") {
         SPC_CONF_LOG << "Load plugin file : " << extra_plugin_file;
         auto retcode = plugin_manager.load_plugin(extra_plugin_file);
         if (not retcode) {
             CRITICAL << "Failed to load plugin file : "
                      << extra_plugin_file << std::endl;
             std::runtime_error("Fail to load plugin file : "
                                + extra_plugin_file + ".");
         }
     }
     auto initializer = plugin_manager.get_object<InitializeVariables>(
                                PLUGIN_MODULE_VARIABLES, plugin_name);
     if (initializer == nullptr) {
         CRITICAL << "No variable initializer found";
         throw std::runtime_error("No variable initializer found");
     }
     return initializer->initialize_variables(
                 raw_data, the_mesh, the_units, configuration);
 }
 
 // ================ //
 //   Plugin Manager //
 // ================ //
 
 plugins::PluginManager& initialize_plugin_manager(
         io::YAMLConfigHandle& config
         )
 {
     SPC_CONF_LOG << SPC_CONF_LOG_SECTION
                  << "\nInitialize plugin manager\n"
                  << SPC_CONF_LOG_HLINE;
     plugins::PluginManager& plugin_manager = plugins::get_plugin_manager();
     SPC_CONF_LOG << "Registering module " << PLUGIN_MODULE_VARIABLES;
     plugin_manager.register_module(PLUGIN_MODULE_VARIABLES,
                                    make_unique<plugins::ModuleBase>());
     SPC_CONF_LOG << "Registering module " << PLUGIN_MODULE_UPSCALING;
     plugin_manager.register_module(PLUGIN_MODULE_UPSCALING,
                                    make_unique<plugins::ModuleBase>());
 
     if (config.has_section(SPC_CF_S_PLUGINS)) {
         io::configure_plugin_manager(config.get_section(SPC_CF_S_PLUGINS));
     }
     auto dirs = plugin_manager.get_plugin_directories();
     if (not dirs.empty()) {
         SPC_CONF_LOG << "plugin path : ";
         for (auto dir: dirs) {
             SPC_CONF_LOG << "  - " << dir;
         }
     }
     SPC_CONF_LOG << SPC_CONF_LOG_HLINE;
     return plugin_manager;
 }
 
 // ================ //
 //    Staggers      //
 // ================ //
 
 Staggers initialize_staggers(
         SimulationData& simul_data,
         io::YAMLConfigHandle&& conf_staggers
         )
 {
     SPC_CONF_LOG << "\n\n" << SPC_CONF_LOG_SECTION
                  << "\n=    Staggers initialization        =\n"
                  << SPC_CONF_LOG_SECTION;
 
     Staggers staggers;
 
     auto* var_ptr = simul_data.variables.get();
 
     staggers.upscaling =
             initialize_upscaling_stagger(simul_data, conf_staggers);
     staggers.upscaling->initialize(var_ptr);
 
     staggers.chemistry =
             initialize_chemistry_stagger(simul_data, conf_staggers);
     staggers.chemistry->initialize(var_ptr);
 
     staggers.transport =
             initialize_transport_stagger(simul_data, conf_staggers);
     staggers.transport->initialize(var_ptr);
 
     staggers.upscaling->register_transport_stagger(staggers.transport);
     staggers.upscaling->register_chemistry_stagger(staggers.chemistry);
 
     // check the database
     if (not simul_data.raw_data->is_valid()) {
         CRITICAL << "Change in database during stagger initialization !";
         throw std::logic_error("Change in the database has been detected"
                                " during the staggers initialization. Abort.");
     }
 
     return staggers;
 }
 
 std::shared_ptr<solver::UpscalingStaggerBase>
 initialize_upscaling_stagger(
         SimulationData& simul_data,
         io::YAMLConfigHandle& config_staggers
         )
 {
     SPC_CONF_LOG << "\n\n --> Upscaling stagger <-- \n"
                  << SPC_CONF_LOG_SECTION;
 
     auto conf = config_staggers.get_section(SPC_CF_S_STAGGERS_SS_UPSCALINGSTAGGER);
     auto type = conf.get_optional_attribute<std::string>(
                 SPC_CF_S_STAGGERS_SS_UPSCALINGSTAGGER_A_TYPE,
                 SPC_CF_V_PLUGIN);
     // only plugin accepted for now
     if (type == SPC_CF_V_DEFAULT) type = SPC_CF_V_PLUGIN;
     if (type != SPC_CF_V_PLUGIN) {
         conf.report_error(
                     io::YAMLConfigError::InvalidArgument,
                     "Invalid argument for attribute '"
                     SPC_CF_A_TYPE
                     "' only value accepted is '" SPC_CF_V_PLUGIN "'.");
     }
     auto plugin_name = conf.get_required_attribute<std::string>(
                 SPC_CF_A_PLUGIN_NAME);
     SPC_CONF_LOG << "Use plugin " << plugin_name;
 
     plugins::PluginManager& plugin_manager = plugins::get_plugin_manager();
     auto extra_plugin_file = conf.get_optional_attribute<std::string>(
                                  SPC_CF_A_PLUGIN_FILE, "");
     if (extra_plugin_file != "") {
         SPC_CONF_LOG << "Load plugin file " << extra_plugin_file;
         auto retcode = plugin_manager.load_plugin(extra_plugin_file);
         if (not retcode) {
             CRITICAL << "Failed to load plugin file : " << extra_plugin_file << std::endl;
             std::runtime_error("Fail to load plugin file : "
                                + extra_plugin_file + ".");
         }
     }
    auto upscaling_factory = plugin_manager.get_object<UpscalingStaggerFactory>(
                                  PLUGIN_MODULE_UPSCALING, plugin_name);
    if (upscaling_factory == nullptr) {
         CRITICAL << "No upscaling stagger factory found";
         throw std::runtime_error("No upscaling stagger factory found");
    }
    SPC_CONF_LOG << SPC_CONF_LOG_HLINE;
    return upscaling_factory->get_upscaling_stagger(
                 simul_data, std::move(conf));
 }
 
 
 std::shared_ptr<solver::ChemistryStaggerBase>
 initialize_chemistry_stagger(
         SimulationData& simul_data,
         io::YAMLConfigHandle& conf_staggers
         )
 {
     SPC_CONF_LOG << "\n\n --> Chemistry stagger <-- \n"
                  << SPC_CONF_LOG_SECTION;
     // no conf
     // -------
     if (not conf_staggers.has_section(SPC_CF_S_STAGGERS_SS_CHEMISTRYSTAGGER)) {
         auto opts = EquilibriumOptionsVector::make(
                                  simul_data.mesh1d->nb_nodes());
         return EquilibriumStagger::make(
                     simul_data.variables,
                     simul_data.bcs,
                     opts
                     );
     }
     // conf is given
     // --------------
     auto conf_chem = conf_staggers.get_section(SPC_CF_S_STAGGERS_SS_CHEMISTRYSTAGGER);
     auto type = conf_chem.get_optional_attribute<std::string>(
                     SPC_CF_S_STAGGERS_SS_CHEMISTRYSTAGGER_A_TYPE,
                     SPC_CF_S_STAGGERS_SS_CHEMISTRYSTAGGER_A_TYPE_V_EQUILIBRIUM);
     if (type == SPC_CF_V_DEFAULT)
         type = SPC_CF_S_STAGGERS_SS_CHEMISTRYSTAGGER_A_TYPE_V_EQUILIBRIUM;
     if (type !=SPC_CF_S_STAGGERS_SS_CHEMISTRYSTAGGER_A_TYPE_V_EQUILIBRIUM) {
         CRITICAL << "Unknown chemistry stagger";
         throw std::invalid_argument("Only equilibrium stagger supported.");
     }
     SPC_CONF_LOG << "Use equilibrium stagger";
     std::shared_ptr<EquilibriumOptionsVector> opts {nullptr};
     if (conf_chem.has_section(SPC_CF_S_STAGGERS_SS_CHEMISTRYSTAGGER_SS_EQUILIBRIUM_OPTS)) {
         opts = io::configure_unsaturated_equilibrium_options(
                    simul_data.mesh1d->nb_nodes(),
                    simul_data.units,
                    conf_chem.get_section(SPC_CF_S_STAGGERS_SS_CHEMISTRYSTAGGER_SS_EQUILIBRIUM_OPTS)
                    );
     } else {
         SPC_CONF_LOG << " Use default options ";
         opts = EquilibriumOptionsVector::make(simul_data.mesh1d->nb_nodes());
         SPC_CONF_LOG << "Default SpecMiCP options : ";
         std::ostringstream msg;
         io::print_specmicp_options(msg, opts->get("default"));
         SPC_CONF_LOG << msg.str() << "\n" << SPC_CONF_LOG_HLINE;
 
     }
     SPC_CONF_LOG << SPC_CONF_LOG_HLINE;
     return EquilibriumStagger::make(
                 simul_data.variables,
                 simul_data.bcs,
                 opts
                 );
 }
 
 std::shared_ptr<solver::TransportStaggerBase> initialize_transport_stagger(
         SimulationData& simul_data,
         io::YAMLConfigHandle& conf_staggers
         )
 {
     SPC_CONF_LOG << "\n\n --> Transport stagger <-- \n"
                  << SPC_CONF_LOG_SECTION;
     const UnsaturatedVariables* const variables = simul_data.variables.get();
     const database::DataContainer* const raw_data = simul_data.raw_data.get();
 
     // No conf
     // -------
     if (not conf_staggers.has_section(SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER)) {
         return UnsaturatedTransportStagger::make(
                     simul_data.variables,
                     simul_data.bcs,
                     true, true);
     }
     // Conf is given
     // -------------
     auto transport_conf = conf_staggers.get_section(
                 SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER);
     const bool merge_sat = transport_conf.get_optional_attribute<bool>(
                 SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_A_MERGE_SATURATION,
                 true);
     if (merge_sat) {
         SPC_CONF_LOG << "Use merged saturation/vapor pressure equation";
     } else {
         SPC_CONF_LOG << "Use separate saturation/vapor pressure equations";
     }
     const bool merge_aq = transport_conf.get_optional_attribute<bool>(
                 SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_A_MERGE_AQUEOUS,
                 true);
     if (merge_aq) {
         SPC_CONF_LOG << "Use merged aqueous liquid/pressure equation";
     } else {
         SPC_CONF_LOG << "Use separate aqueous liquid/pressure equations";
     }
 
     auto stagger = UnsaturatedTransportStagger::make(
                        simul_data.variables,
                        simul_data.bcs,
                        merge_sat, merge_aq);
 
     // ! scaling => done in variables
 
     // default options
     // ---------------
     dfpmsolver::ParabolicDriverOptions default_opts;
     bool rewrite_default = false;
     if (transport_conf.has_section(SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_SS_DEFAULT_OPTS))
     {
         io::configure_transport_options(
                     default_opts,
                     transport_conf.get_section(SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_SS_DEFAULT_OPTS)
                     );
         rewrite_default = true;
     }
     SPC_CONF_LOG << "\nDefault transport options\n"
                  << SPC_CONF_LOG_HLINE;
     std::ostringstream msg;
     io::print_transport_options(default_opts, msg);
     SPC_CONF_LOG << msg.str() << SPC_CONF_LOG_HLINE;
 
 
     // Saturation option
     // -----------------
     auto& sat_opts = *(stagger->get_saturation_options());
     sat_opts = default_opts;
     if (transport_conf.has_section(SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_SS_SATURATION_OPTS))
     {
         io::configure_transport_options(
                     sat_opts,
                     transport_conf.get_section(SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_SS_SATURATION_OPTS)
                     );
     }
     SPC_CONF_LOG << "\nTransport options for saturation equations \n"
                  << SPC_CONF_LOG_HLINE;
     msg = std::ostringstream();
     io::print_transport_options(sat_opts, msg);
     SPC_CONF_LOG << msg.str() << SPC_CONF_LOG_HLINE;
 
     // Aqueous options
     // ---------------
     if (transport_conf.has_section(SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_SS_AQUEOUS_OPTS))
     {
         auto conf = transport_conf.get_section(SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_SS_AQUEOUS_OPTS);
         if (not conf.is_sequence()) {
             // config for every aqueous equations
             dfpmsolver::ParabolicDriverOptions opts;
             if (rewrite_default) opts = default_opts;
             io::configure_transport_options(
                         opts, std::move(conf));
             for (index_t aq: raw_data->range_aqueous_component()) {
                 *(stagger->get_aqueous_options(aq)) = opts;
             }
             SPC_CONF_LOG << "\nTransport options for aqueous equations \n"
                          << SPC_CONF_LOG_HLINE;
             msg = std::ostringstream();
             io::print_transport_options(opts, msg);
             SPC_CONF_LOG << msg.str() << SPC_CONF_LOG_HLINE;
 
         }
         else {
             // First set default is needed
             if (rewrite_default) {
                 for (index_t aq: simul_data.raw_data->range_aqueous_component()) {
                     *(stagger->get_aqueous_options(aq)) = default_opts;
                 }
             }
             uindex_t size = conf.size();
             for (uindex_t ind=0; ind<size; ++ind)
             {
                 auto subconf = conf.get_section(ind);
                 // personalized options per component
                 auto comp_label = subconf.get_required_attribute<std::string>(
                           SPC_CF_A_COMPONENT);
                 auto comp_id = raw_data->get_id_component(comp_label);
                 io::configure_transport_options(
                             *(stagger->get_aqueous_options(comp_id)),
                             std::move(subconf)
                             );
             }
         }
 
     }
     // Gas options
     // -----------
     if (    transport_conf.has_section(SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_SS_GAS_OPTS)
         and not (merge_aq and merge_sat))  // skip conf if not needed
     {
         auto conf = transport_conf.get_section(SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_SS_GAS_OPTS);
         if (not conf.is_sequence()) {
             // config for every gas equations
             dfpmsolver::ParabolicDriverOptions opts;
             if (rewrite_default) opts = default_opts;
             io::configure_transport_options(
                         opts, std::move(conf));
             if ( not merge_sat and variables->component_has_gas(0)) {
                 *(stagger->get_gas_options(0)) = opts;
             }
             if ( not merge_aq) {
                 for (index_t aq: raw_data->range_aqueous_component()) {
                     if (simul_data.variables->component_has_gas(aq) )
                         *(stagger->get_aqueous_options(aq)) = opts;
                 }
             }
             SPC_CONF_LOG << "\nTransport options for gaseous equations \n"
                          << SPC_CONF_LOG_HLINE;
             msg = std::ostringstream();
             io::print_transport_options(opts, msg);
             SPC_CONF_LOG << msg.str() << SPC_CONF_LOG_HLINE;
         }
         else
         {
             // First rewrite default
             if (rewrite_default) {
                 // H2O
                 if (not merge_sat and variables->component_has_gas(0)) {
                     *(stagger->get_gas_options(0)) = default_opts;
                 }
                 // aqueous components
                 if (not merge_aq)
                 {
                     for (index_t aq: raw_data->range_aqueous_component()) {
                         if (simul_data.variables->component_has_gas(aq))
                             *(stagger->get_aqueous_options(aq)) = default_opts;
                     }
                 }
             }
             // Then go over all conf
             uindex_t size = conf.size();
             for (uindex_t ind=0; ind<size; ++ind)
             {
                 auto subconf = conf.get_section(ind);
                 auto comp_label = subconf.get_required_attribute<std::string>(
                       SPC_CF_A_COMPONENT);
                 auto comp_id = raw_data->get_id_component(comp_label);
                 // H2O
                 if (comp_id == 0) {
                     if (merge_sat) continue;
                     io::configure_transport_options(
                                 *(stagger->get_gas_options(0)),
                                 std::move(subconf)
                                 );
 
                 }
                 // aqueous component
                 else {
                     if (merge_aq) continue;
                     io::configure_transport_options(
                                 *(stagger->get_gas_options(comp_id)),
                                 std::move(subconf)
                                 );
                 }
             }
         }
     }
     return stagger;
 }
 
 
 std::string initialize_output(
         SimulationData& simul_data,
         const std::string& work_dir,
         io::YAMLConfigHandle&& conf
         )
 {
     auto type = conf.get_optional_attribute<std::string>(
                 SPC_CF_S_REACTOUTPUT_A_TYPE, SPC_CF_S_REACTOUTPUT_A_TYPE_V_HDF5);
     if (type == SPC_CF_V_DEFAULT) type = SPC_CF_S_REACTOUTPUT_A_TYPE_V_HDF5;
     if (type != SPC_CF_S_REACTOUTPUT_A_TYPE_V_HDF5) {
         conf.report_error(
                     io::YAMLConfigError::InvalidArgument,
                     "Only hdf5 is accepted for attribute "
                      SPC_CF_S_REACTOUTPUT_A_TYPE " for now.");
     }
     auto filepath = conf.get_required_attribute<std::string>(
                 SPC_CF_S_REACTOUTPUT_A_FILEPATH);
     if (not work_dir.empty() and not utils::is_path_absolute(filepath)) {
         filepath = utils::complete_path(work_dir, filepath);
     }
     if (conf.has_attribute(SPC_CF_S_REACTOUTPUT_A_DATABASE)) {
         database::Database db_manager(simul_data.raw_data);
         auto db_path = conf.get_attribute<std::string>(SPC_CF_S_REACTOUTPUT_A_DATABASE);
         if (not work_dir.empty() and not utils::is_path_absolute(db_path)) {
             db_path = utils::complete_path(work_dir, db_path);
         }
         db_manager.save(db_path);
     }
     return filepath;
 }
 
 
 void print_db_to_conf_log(database::DataContainer* raw_data)
 {
     SPC_CONF_LOG << "Freeze database\n"
                  << SPC_CONF_LOG_SECTION;
     SPC_CONF_LOG <<   " - " << raw_data->nb_component() << " components"
                  << "\n - " << raw_data->nb_aqueous()   << " secondary aqueous species"
                  << "\n - " << raw_data->nb_mineral()   << " solid phases at equilibrium"
                  << "\n - " << raw_data->nb_gas()       << " gas";
     SPC_CONF_LOG << "List component : ";
     for (auto comp: raw_data->range_component()) {
         SPC_CONF_LOG << "   -  " << raw_data->get_label_component(comp);
     }
     SPC_CONF_LOG << SPC_CONF_LOG_SECTION;
 }
diff --git a/src/reactmicp/systems/unsaturated/variables.cpp b/src/reactmicp/systems/unsaturated/variables.cpp
index 7d20dee..6b5b653 100644
--- a/src/reactmicp/systems/unsaturated/variables.cpp
+++ b/src/reactmicp/systems/unsaturated/variables.cpp
@@ -1,408 +1,410 @@
 /* =============================================================================
 
  Copyright (c) 2014 - 2016
  F. Georget <fabieng@princeton.edu> Princeton University
  All rights reserved.
 
 
 Redistribution and use in source and binary forms, with or without modification,
 are permitted provided that the following conditions are met:
 
     1. Redistributions of source code must retain the above copyright notice,
     this list of conditions and the following disclaimer.
 
     2. Redistributions in binary form must reproduce the above copyright notice,
     this list of conditions and the following disclaimer in the documentation
     and/or other materials provided with the distribution.
 
     3. Neither the name of the copyright holder nor the names of its
     contributors may be used to endorse or promote products derived from this
     software without specific prior written permission.
 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
 ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
 ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. *
 
 ============================================================================= */
 
 #include "variables.hpp"
 #include "dfpm/meshes/mesh1d.hpp"
 
 #include "variables_box.hpp"
 
 #include "specmicp_common/log.hpp"
 
 #include <stdexcept>
 
 namespace specmicp {
 namespace reactmicp {
 namespace systems {
 namespace unsaturated {
 
 //  ===================
 //  List main variables
 //  ===================
 
 ListMainVariable::ListMainVariable(index_t nb_vars, index_t nb_nodes)
 {
     variables.reserve(nb_vars);
     variables.emplace_back(make_unique<MainVariable>(nb_nodes));
     variables.emplace_back(nullptr);
     for (index_t i=2; i<nb_vars; ++i)
     {
         variables.emplace_back(make_unique<MainVariable>(nb_nodes));
     }
 }
 
 ListMainVariable::ListMainVariable(
         index_t nb_vars,
         index_t nb_nodes,
         std::vector<bool> to_init
         )
 {
     variables.reserve(nb_vars);
     for (const auto& is_present: to_init) {
         if (is_present) {
             variables.emplace_back(make_unique<MainVariable>(nb_nodes));
         } else {
             variables.emplace_back(nullptr);
         }
     }
 }
 
 //! \brief Reset the variables
 void ListMainVariable::reset()
 {
     //if (variables[0] != nullptr)
     //    variables[0]->reset();
     for (vector_type::size_type i=0; i<variables.size(); ++i)
     {
         if (variables[i] != nullptr)
             variables[i]->reset();
     }
 }
 
 
 void ListMainVariable::hard_reset(index_t component, index_t nb_nodes)
 {
     specmicp_assert(component >= 0 and (unsigned) component < variables.size());
     if (nb_nodes == 0)
         variables[component].reset(nullptr);
     else
         variables[component].reset(new MainVariable(nb_nodes));
 }
 
 
 
 //  ==============
 //  Variables box
 //  ==============
 
 SaturationVariableBox::SaturationVariableBox(UnsaturatedVariables& vars):
     liquid_saturation(vars.m_liquid_var.water()),
     solid_concentration(vars.m_solid_var.water()),
     vapor_pressure(vars.m_gas_var.water()),
     aqueous_concentration(vars.m_water_aq_concentration),
     porosity(vars.m_porosity),
     liquid_permeability(vars.m_liquid_permeability),
     liquid_diffusivity(vars.m_liquid_diffusivity),
     relative_liquid_permeability(vars.m_relative_liquid_permeability),
     relative_liquid_diffusivity(vars.m_relative_liquid_diffusivity),
     capillary_pressure(vars.m_capillary_pressure),
     advection_flux(vars.m_advection_flux),
     constants(vars.m_constants),
     capillary_pressure_f(vars.m_capillary_pressure_f),
     relative_liquid_permeability_f(vars.m_relative_liquid_permeability_f),
     relative_liquid_diffusivity_f(vars.m_relative_liquid_diffusivity_f)
 {}
 
 SaturationPressureVariableBox::SaturationPressureVariableBox(
         UnsaturatedVariables& vars):
     binary_diffusion_coefficient(vars.get_binary_gas_diffusivity(0)),
     liquid_saturation(vars.m_liquid_var.water()),
     partial_pressure(vars.get_pressure_main_variables(0)),
     solid_concentration(vars.m_solid_var.water()),
     aqueous_concentration(vars.m_water_aq_concentration),
     porosity(vars.m_porosity),
     liquid_permeability(vars.m_liquid_permeability),
     liquid_diffusivity(vars.m_liquid_diffusivity),
     relative_liquid_permeability(vars.m_relative_liquid_permeability),
     relative_liquid_diffusivity(vars.m_relative_liquid_diffusivity),
     capillary_pressure(vars.m_capillary_pressure),
     resistance_gas_diffusivity(vars.m_resistance_gas_diffusivity),
     relative_gas_diffusivity(vars.m_relative_gas_diffusivity),
     advection_flux(vars.m_advection_flux),
     constants(vars.m_constants),
     capillary_pressure_f(vars.m_capillary_pressure_f),
     relative_liquid_permeability_f(vars.m_relative_liquid_permeability_f),
     relative_liquid_diffusivity_f(vars.m_relative_liquid_diffusivity_f),
     relative_gas_diffusivity_f(vars.m_relative_gas_diffusivity_f),
     partial_pressure_f(vars.m_vapor_pressure_f)
 {}
 
 LiquidAqueousComponentVariableBox::LiquidAqueousComponentVariableBox(
         UnsaturatedVariables& vars,
         index_t component
         ):
     aqueous_concentration(vars.m_liquid_var.aqueous_component(component)),
     solid_concentration(vars.m_solid_var.aqueous_component(component)),
     partial_pressure(vars.get_pressure_main_variables(component)),
     saturation(vars.m_liquid_var.water()),
     porosity(vars.m_porosity),
     liquid_diffusivity(vars.m_liquid_diffusivity),
     relative_liquid_diffusivity(vars.m_relative_liquid_diffusivity),
     advection_flux(vars.m_advection_flux),
     constants(vars.m_constants)
 {}
 
 
 LiquidGasAqueousVariableBox::LiquidGasAqueousVariableBox(
         UnsaturatedVariables& vars,
         index_t component
         ):
     LiquidAqueousComponentVariableBox(vars, component),
     binary_diffusion_coefficient(vars.get_binary_gas_diffusivity(component)),
     resistance_gas_diffusivity(vars.m_resistance_gas_diffusivity),
     relative_gas_diffusivity(vars.m_relative_gas_diffusivity)
 {}
 
 PressureVariableBox::PressureVariableBox(
         UnsaturatedVariables& vars,
         index_t component
         ):
     binary_diffusion_coefficient(vars.get_binary_gas_diffusivity(component)),
     partial_pressure(vars.get_pressure_main_variables(component)),
     liquid_saturation(vars.m_liquid_var.water()),
     porosity(vars.m_porosity),
     resistance_gas_diffusivity(vars.m_resistance_gas_diffusivity),
     relative_gas_diffusivity(vars.m_relative_gas_diffusivity),
     constants(vars.m_constants)
 {}
 
-void ConstantBox::scale(units::LengthUnit length_unit)
+
+void ConstantBox::scale(const units::UnitsSet& units_set)
 {
-    scalar_t scaling_factor = units::scaling_factor(length_unit);
+    scalar_t scaling_factor_l = units::scaling_factor(units_set.length);
+    scalar_t scaling_factor_m = units::scaling_factor(units_set.quantity);
     //viscosity_liquid_water *= 1.0/scaling_factor;
-    rt /= std::pow(scaling_factor, 3);
+    rt *= scaling_factor_m / std::pow(scaling_factor_l, 3);
     //total_pressure *= 1.0/scaling_factor;
 
 }
 
 //  ======================
 //  Unsaturated variables
 //  ======================
 
 namespace internal {
 
 //! \brief Find the correspond gas id for a component with gas
 static index_t get_id_gas(database::RawDatabasePtr the_database,
                           index_t component)
 {
     index_t id = no_species;
     for (auto gas: the_database->range_gas())
     {
         const scalar_t nu = the_database->nu_gas(gas, component);
         if (nu == 0.0)
             continue;
 
         if (nu != 1.0)
         {
             ERROR << "Database not in the good format"
                   << "; expect : gas <=> component \n"
                   << "component : "
                   << the_database->get_label_component(component) << " \n"
                   << "gas : "
                   << the_database->get_label_gas(gas) << ".";
             throw std::runtime_error("Badly formatted database");
         }
         else
         {
             id = gas;
         }
     }
     if (id == no_species)
     {
         ERROR << "No corresponding gas in the database for component : "
               << the_database->get_label_component(component);
         throw std::runtime_error("Badly formatted database");
     }
     return id;
 }
 
 } //end namespace internal
 
 UnsaturatedVariables::UnsaturatedVariables(
         mesh::Mesh1DPtr the_mesh,
         database::RawDatabasePtr the_database,
         std::vector<bool> has_gas
         ):
     database::DatabaseHolder(the_database),
     m_mesh(the_mesh),
     m_has_gas(has_gas),
     m_id_gas(the_database->nb_component(), no_species),
     // main var
     m_liquid_var(the_database->nb_component(),the_mesh->nb_nodes()),
     m_gas_var(the_database->nb_component(),the_mesh->nb_nodes(), has_gas),
     m_solid_var(the_database->nb_component(),the_mesh->nb_nodes()),
     // chem sol
     m_chem_sol(the_mesh->nb_nodes()),
     // secondary vars
     m_porosity(the_mesh->nb_nodes()),
     m_water_aq_concentration(the_mesh->nb_nodes()),
     m_liquid_permeability(the_mesh->nb_nodes()),
     m_relative_liquid_permeability(the_mesh->nb_nodes()),
     m_liquid_diffusivity(the_mesh->nb_nodes()),
     m_relative_liquid_diffusivity(the_mesh->nb_nodes()),
     m_binary_gas_diffusivity(Vector::Zero(the_database->nb_component())),
     m_resistance_gas_diffusivity(the_mesh->nb_nodes()),
     m_relative_gas_diffusivity(the_mesh->nb_nodes()),
     m_capillary_pressure(the_mesh->nb_nodes()),
     m_advection_flux(the_mesh->nb_nodes()),
     m_aqueous_scaling(Vector::Ones(the_database->nb_component())),
     m_gaseous_scaling(Vector::Ones(the_database->nb_component()))
 {
     specmicp_assert(has_gas.size() == (unsigned) the_database->nb_component());
     //
     // trick to avoid more branching in equations :
     m_gas_var.hard_reset(1, the_mesh->nb_nodes());
 
     for (auto component: the_database->range_aqueous_component())
     {
         if (has_gas[component])
         {
             m_id_gas[component] = internal::get_id_gas(the_database, component);
         }
     }
 }
 
-
 UnsaturatedVariables::UnsaturatedVariables(
         mesh::Mesh1DPtr the_mesh,
         database::RawDatabasePtr the_database,
         std::vector<bool> has_gas,
-        units::LengthUnit length_unit):
+        const units::UnitsSet& units_set):
     UnsaturatedVariables(the_mesh, the_database, has_gas)
 {
-    m_constants.scale(length_unit);
+    m_constants.scale(units_set);
 }
 
+
 void UnsaturatedVariables::reset_main_variables()
 {
     m_liquid_var.reset();
     m_gas_var.reset();
     m_solid_var.reset();
 
     m_porosity.variable = m_porosity.predictor;
     m_porosity.velocity.setZero();
 
     m_water_aq_concentration.variable = m_water_aq_concentration.predictor;
     m_water_aq_concentration.velocity.setZero();
 
     m_advection_flux.set_zero();
 
     set_relative_variables();
 
 }
 
 SaturationVariableBox
 UnsaturatedVariables::get_saturation_variables() {
     return SaturationVariableBox(*this);
 }
 
 SaturationPressureVariableBox
 UnsaturatedVariables::get_saturation_pressure_variables() {
     return SaturationPressureVariableBox(*this);
 }
 
 PressureVariableBox
 UnsaturatedVariables::get_vapor_pressure_variables() {
     return PressureVariableBox(*this, 0);
 }
 
 
 LiquidAqueousComponentVariableBox
 UnsaturatedVariables::get_liquid_aqueous_component_variables(
         index_t component
         ) {
     return LiquidAqueousComponentVariableBox(*this, component);
 }
 
 LiquidGasAqueousVariableBox
 UnsaturatedVariables::get_liquid_gas_aqueous_variables(
         index_t component
         ) {
     return LiquidGasAqueousVariableBox(*this, component);
 }
 
 PressureVariableBox
 UnsaturatedVariables::get_pressure_variables(
         index_t component
         ) {
     return PressureVariableBox(*this, component);
 }
 
 
 MainVariable& UnsaturatedVariables::get_pressure_main_variables(index_t component)
 {
     if (component_has_gas(component))
     {
         specmicp_assert(m_gas_var(component) != nullptr);
         return *m_gas_var(component);
     }
     else
         return *m_gas_var(1); // a zero var..., this is a trick to always get a valid info
 }
 
 index_t UnsaturatedVariables::nb_gas() {
     return std::count(m_has_gas.begin(), m_has_gas.end(), true);
 }
 
 void UnsaturatedVariables::set_relative_variables(
         index_t node
         )
 {
     const scalar_t saturation = m_liquid_var.water().variable(node);
 
     m_relative_liquid_diffusivity(node)  =
             m_relative_liquid_diffusivity_f(node, saturation);
     m_relative_gas_diffusivity(node)     =
             m_relative_gas_diffusivity_f(node, saturation);
     m_relative_liquid_permeability(node) =
             m_relative_liquid_permeability_f(node, saturation);
     m_capillary_pressure(node)           =
             m_capillary_pressure_f(node, saturation);
 
 }
 
 void UnsaturatedVariables::set_relative_variables()
 {
     for (index_t node: m_mesh->range_nodes())
     {
         set_relative_variables(node);
     }
 }
 
 void UnsaturatedVariables::set_vapor_pressure(index_t node)
 {
     if (component_has_gas(0))
     {
         const scalar_t saturation = m_liquid_var.water().variable(node);
         m_gas_var(0)->variable(node) = m_vapor_pressure_f(node, saturation);
     }
 
 
 }
 
 
 
 } //end namespace unsatura
 } //end namespace systems
 } //end namespace reactmicp
 } //end namespace specmicp
diff --git a/src/reactmicp/systems/unsaturated/variables.hpp b/src/reactmicp/systems/unsaturated/variables.hpp
index 3e9d873..db080dc 100644
--- a/src/reactmicp/systems/unsaturated/variables.hpp
+++ b/src/reactmicp/systems/unsaturated/variables.hpp
@@ -1,420 +1,417 @@
 /* =============================================================================
 
  Copyright (c) 2014 - 2016
  F. Georget <fabieng@princeton.edu> Princeton University
  All rights reserved.
 
 
 Redistribution and use in source and binary forms, with or without modification,
 are permitted provided that the following conditions are met:
 
     1. Redistributions of source code must retain the above copyright notice,
     this list of conditions and the following disclaimer.
 
     2. Redistributions in binary form must reproduce the above copyright notice,
     this list of conditions and the following disclaimer in the documentation
     and/or other materials provided with the distribution.
 
     3. Neither the name of the copyright holder nor the names of its
     contributors may be used to endorse or promote products derived from this
     software without specific prior written permission.
 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
 ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
 ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. *
 
 ============================================================================= */
 
 #ifndef SPECMICP_REACTMICP_UNSATURATED_VARIABLES_HPP
 #define SPECMICP_REACTMICP_UNSATURATED_VARIABLES_HPP
 
 //! \file unsaturated/variables.hpp
 //! \name Variables for the unsaturated system
 
 #include "types_fwd.hpp"
 #include "dfpm/meshes/mesh1d.hpp"
 #include "specmicp_database/database_holder.hpp"
 
 #include "variables_sub.hpp"
 
 #include "reactmicp/solver/staggers_base/variables_base.hpp"
 
 #include <functional>
 
 namespace specmicp {
 namespace reactmicp {
 namespace systems {
 
 //! \namespace unsaturated
 //! \brief The unsaturated system
 //!
 //! The system solve the transport reactive in unsaturated porous medium
 namespace unsaturated {
 
 
 // forward declarations
 // declared in variables_box.hpp
 struct SaturationVariableBox;
 struct SaturationPressureVariableBox;
 struct PressureVariableBox;
 struct LiquidAqueousComponentVariableBox;
 struct LiquidGasAqueousVariableBox;
 
 
 //! \brief The unsaturated variables
 //!
 //! This class contains all the variables for the unsaturated system
 //!
 //! Note : it should be initialized using the ::VariablesInterface class.
 class SPECMICP_DLL_PUBLIC UnsaturatedVariables:
         public database::DatabaseHolder,
         public solver::VariablesBase
 {
 public:
     //! \brief Constructor
     //! \internal
     UnsaturatedVariables(
             mesh::Mesh1DPtr the_mesh,
             database::RawDatabasePtr the_database,
             std::vector<bool> has_gas);
-    //! \brief Constructor
     //! \internal
     UnsaturatedVariables(
             mesh::Mesh1DPtr the_mesh,
             database::RawDatabasePtr the_database,
             std::vector<bool> has_gas,
-            units::LengthUnit length_unit);
-
-
+            const units::UnitsSet& units_set);
 
     //! \brief Reset the main variables in case of failure
     void reset_main_variables() override;
 
     //! \brief Return the number of governing equations
     //!
     //! The governing equations are the transport equations
     index_t nb_governing_equations() const {
         return get_database()->nb_aqueous_components()+1;
     }
 
     //! \brief Return the number of nodes
     index_t nb_nodes() const {
         return m_mesh->nb_nodes();
     }
 
     //! \brief Return the mesh
     mesh::Mesh1DPtr get_mesh() {
         return m_mesh;
     }
 
     //! \brief Return true if a component has a gas equation
     //! \param component water or aqueous component
     bool component_has_gas(index_t component) const {
         specmicp_assert_component_bounds(component, get_database());
         return m_has_gas[component];
     }
 
     //! \brief Return the number of gas in the system
     index_t nb_gas();
 
     // Variables Box
     // =============
 
     //! \name VariableBox
     //! \brief This set of variables are passed to the equations
     //!
     //! They contain all the variables a governing equation should need
     //! @{
 
     //! \brief Return the variables for the saturation equation
     SaturationVariableBox get_saturation_variables();
 
     //! \brief Return the variables for the saturation-pressure equation
     SaturationPressureVariableBox get_saturation_pressure_variables();
 
     //! \brief Return the variables for the vapor pressure equation
     PressureVariableBox get_vapor_pressure_variables();
 
     //! \brief Return the variables for the liquid transport of an
     //!  aqueous component
     LiquidAqueousComponentVariableBox get_liquid_aqueous_component_variables(
             index_t component);
 
     //! \brief Return the variables for the liquid and gas transport of an
     //!  aqueous component
     LiquidGasAqueousVariableBox get_liquid_gas_aqueous_variables(
             index_t component);
 
     //! \brief Return the variables for the pressure diffusion equation
     //!
     //! Valid for  water vapor or the gas corresponding to an aqueous component
     PressureVariableBox get_pressure_variables(index_t component);
 
     //! @}
 
     // User model
     // ==========
 
     //! \name User models
     //! \brief Models given by the user to compute saturated-dependant variables
     //!
     //! A model is a function taking a node and the saturation as the argument
     //! @{
 
 
     //! \brief Set the capillary pressure model
     void set_capillary_pressure_model(user_model_saturation_f func) {
         m_capillary_pressure_f = func;
     }
 
     //! \brief Return the capillary pressure model
     user_model_saturation_f get_capillary_pressure_model() {
         return m_capillary_pressure_f;
     }
 
     //! \brief Set the vapor pressure model
     void set_vapor_pressure_model(user_model_saturation_f func) {
         m_vapor_pressure_f = func;
     }
 
     //! \brief Return the vapor pressure model
     user_model_saturation_f get_vapor_pressure_model() {
         return m_vapor_pressure_f;
     }
 
     //! \brief Set the relative liquid permeability model
     void set_relative_liquid_permeability_model(user_model_saturation_f func) {
         m_relative_liquid_permeability_f = func;
     }
 
     //! \brief Set the relative liquid diffusivity model
     void set_relative_liquid_diffusivity_model(user_model_saturation_f func) {
         m_relative_liquid_diffusivity_f = func;
     }
 
     //! \brief Set the relative gas diffusivity model
     void set_relative_gas_diffusivity_model(user_model_saturation_f func) {
         m_relative_gas_diffusivity_f = func;
     }
 
     //! \brief Set the relative variables at 'node'
     void set_relative_variables(index_t node);
     //! \brief Set the relative variables
     void set_relative_variables();
 
     //! \brief Set the partial pressure at 'node'
     void set_vapor_pressure(index_t node);
 
     //! @}
 
 
 
     // Variables Getter
     // ================
 
     //! \name Variables Getter
     //! \brief Return the variables
 
     //! @{
 
     //! \brief Return a reference to the porosity variables
     SecondaryTransientVariable& get_porosity() {
         return m_porosity;}
     //! \brief Return a reference to the water aq. concentration variables
     SecondaryTransientVariable& get_water_aqueous_concentration() {
         return m_water_aq_concentration;}
     //! \brief Return a reference to the liquid permeability
     SecondaryVariable& get_liquid_permeability() {
         return m_liquid_permeability;
     }
     //! \brief Return a reference to the capillary pressure
     SecondaryVariable& get_capillary_pressure() {
         return m_capillary_pressure;
     }
 
     //! \brief Return a reference to the relative liquid permeability
     SecondaryVariable& get_relative_liquid_permeability() {
         return m_relative_liquid_permeability;
     }
     //! \brief Return a reference to the liquid diffusivity
     SecondaryVariable& get_liquid_diffusivity() {
         return m_liquid_diffusivity;
     }
     //! \brief Return a reference to the relative liquid diffusivity
     SecondaryVariable& get_relative_liquid_diffusivity() {
         return m_relative_liquid_diffusivity;
     }
 
     scalar_t & get_binary_gas_diffusivity(index_t component) {
         return m_binary_gas_diffusivity(component);
     }
 
     //! \brief Return a reference to the gas diffusivity resistance factor
     SecondaryVariable& get_resistance_gas_diffusivity() {
         return m_resistance_gas_diffusivity;
     }
     //! \brief Return a reference to the relative gas diffusivity
     SecondaryVariable& get_relative_gas_diffusivity() {
         return m_relative_gas_diffusivity;
     }
     SecondaryVariable& get_advection_flux() {
         return m_advection_flux;
     }
 
     //! \brief Return the liquid saturation
     MainVariable& get_liquid_saturation() {
         return m_liquid_var.water();
     }
 
     //! \brief Return the variables corresponding to the pressure
     //! \param component water or aqueous component
     MainVariable& get_pressure_main_variables(index_t component);
 
     //! \brief Return the aqueous concentration of an aqueous component
     MainVariable& get_aqueous_concentration(index_t aq_component) {
         specmicp_assert(aq_component >= 2 and aq_component < get_database()->nb_component());
         return m_liquid_var.aqueous_component(aq_component);
     }
 
     //! \brief Return the solid variables
     MainVariable& get_solid_concentration(index_t component) {
         specmicp_assert_component_bounds(component, get_database());
         return m_solid_var.component(component);
     }
 
     //! \brief Return the R*T product
     //!
     //! R : ideal gas constant
     //! T : temperature
     scalar_t get_rt() {return m_constants.rt;}
 
     //! \brief return the total pressure
     scalar_t get_total_pressure() {return m_constants.total_pressure;}
 
     //! \brief return the id of a gas
     scalar_t get_id_gas(index_t component) {
         specmicp_assert_component_bounds(component, get_database());
         return m_id_gas[component];
     }
 
     AdimensionalSystemSolution& get_adim_solution(index_t node) {
         return m_chem_sol.solution(node);
     }
 
     void set_adim_solution(index_t node,
                            const AdimensionalSystemSolution& new_solution
                            ) {
         m_chem_sol.update_solution(node, new_solution);
     }
 
     //! @}
 
     //! \name Scaling
     //! \brief Scaling of equations
     //!
     //! @{
 
 
     //! \brief Return the scaling factor for aqueous equation of 'component'
     scalar_t get_aqueous_scaling(index_t component) {
         return m_aqueous_scaling(component);
     }
     //! \brief set the scaling factor for aqueous equation of 'component'
     void set_aqueous_scaling(index_t component, scalar_t value) {
         m_aqueous_scaling(component) = value;
     }
     //! \brief Return the scaling factor for aqueous equation of 'component'
     scalar_t get_gaseous_scaling(index_t component) {
         return m_gaseous_scaling(component);
     }
     //! \brief set the scaling factor for aqueous equation of 'component'
     void set_gaseous_scaling(index_t component, scalar_t value) {
         m_gaseous_scaling(component) = value;
     }
 
 
     //! @}
 
 private:
     mesh::Mesh1DPtr m_mesh;
 
     std::vector<bool> m_has_gas;
     std::vector<index_t> m_id_gas;
 
     // main variables
     // --------------
     ListMainVariable m_liquid_var;
     ListMainVariable m_gas_var;
     ListMainVariable m_solid_var;
 
     // chemistry
     // ---------
     ChemistrySolutions m_chem_sol;
 
     // secondary variables
     // -------------------
     SecondaryTransientVariable m_porosity;
     SecondaryTransientVariable m_water_aq_concentration;
 
     SecondaryVariable m_liquid_permeability;
     SecondaryVariable m_relative_liquid_permeability;
     SecondaryVariable m_liquid_diffusivity;
     SecondaryVariable m_relative_liquid_diffusivity;
 
     Vector m_binary_gas_diffusivity;
     SecondaryVariable m_resistance_gas_diffusivity;
     SecondaryVariable m_relative_gas_diffusivity;
 
     SecondaryVariable m_capillary_pressure;
 
     SecondaryVariable m_advection_flux;
 
     // scaled constants
     // ----------------
     ConstantBox m_constants;
 
     // user models
     // -----------
     user_model_saturation_f m_capillary_pressure_f {nullptr};
     user_model_saturation_f m_vapor_pressure_f {nullptr};
 
     user_model_saturation_f m_relative_liquid_permeability_f {nullptr};
     user_model_saturation_f m_relative_liquid_diffusivity_f {nullptr};
     user_model_saturation_f m_relative_gas_diffusivity_f {nullptr};
 
     // vector of scaling
     Vector m_aqueous_scaling;
     Vector m_gaseous_scaling;
 
     // These struct are declared as friends,
     // To make their initialisation easier
     // (initialisation of referenece, and const reference)
     friend struct SaturationVariableBox;
     friend struct SaturationPressureVariableBox;
     friend struct PressureVariableBox;
     friend struct LiquidAqueousComponentVariableBox;
     friend struct LiquidGasAqueousVariableBox;
 };
 
 //! \brief Shared pointer the unsaturated variables
 using UnsaturatedVariablesPtr = std::shared_ptr<UnsaturatedVariables>;
 
 //! \brief Cast from base variables
 inline UnsaturatedVariablesPtr cast_from_base(
         std::shared_ptr<solver::VariablesBase> ptr) {
     return std::static_pointer_cast<UnsaturatedVariables>(ptr);
 }
 
 } //end namespace unsaturated
 } //end namespace systems
 } //end namespace reactmicp
 } //end namespace specmicp
 
 #endif // SPECMICP_REACTMICP_UNSATURATED_VARIABLES_HPP
diff --git a/src/reactmicp/systems/unsaturated/variables_interface.cpp b/src/reactmicp/systems/unsaturated/variables_interface.cpp
index cb05f0d..9ab6096 100644
--- a/src/reactmicp/systems/unsaturated/variables_interface.cpp
+++ b/src/reactmicp/systems/unsaturated/variables_interface.cpp
@@ -1,820 +1,820 @@
 /* =============================================================================
 
  Copyright (c) 2014 - 2016
  F. Georget <fabieng@princeton.edu> Princeton University
  All rights reserved.
 
 
 Redistribution and use in source and binary forms, with or without modification,
 are permitted provided that the following conditions are met:
 
     1. Redistributions of source code must retain the above copyright notice,
     this list of conditions and the following disclaimer.
 
     2. Redistributions in binary form must reproduce the above copyright notice,
     this list of conditions and the following disclaimer in the documentation
     and/or other materials provided with the distribution.
 
     3. Neither the name of the copyright holder nor the names of its
     contributors may be used to endorse or promote products derived from this
     software without specific prior written permission.
 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
 ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
 ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. *
 
 ============================================================================= */
 
 #include "variables_interface.hpp"
 
 #include "variables.hpp"
 #include "variables_box.hpp"
 
 #include "reactmicp/solver/staggers_base/variables_base.hpp"
 
 #include "specmicp/adimensional/adimensional_system_solution_extractor.hpp"
 
 #include "specmicp_common/compat.hpp"
 #include "specmicp_common/log.hpp"
 #include "specmicp_common/eigen/value_checker.hpp"
 #include "specmicp_common/eigen/vector_checker.hpp"
 
 #include <stdexcept>
 
 // This is a boring class, with only small differences between
 // all the functions and variables. But, they make all the
 // difference at the end. The main purpose of the class is to
 // hide all the variables structure and only present a unified
 // interface to the user.
 
 // None of the structure defined in "variables_sub.hpp" should
 // be leaked outside this class
 
 namespace specmicp {
 namespace reactmicp {
 namespace systems {
 namespace unsaturated {
 
 using namespace specmicp::vector_checker;
 
 
 // typedef
 // -------
 
 using VariablesBase = solver::VariablesBase;
 using VariablesBasePtr = std::shared_ptr<VariablesBase>;
 
 // =================================
 //
 // Declaration of the implementation
 //
 // =================================
 
 // This class just store and retrieve variables group
 
 struct VariablesInterface::VariablesInterfaceImpl
 {
     UnsaturatedVariablesPtr m_vars;
     mesh::Mesh1DPtr m_mesh;
     database::RawDatabasePtr m_database;
 
 
     VariablesInterfaceImpl(
             mesh::Mesh1DPtr the_mesh,
             database::RawDatabasePtr the_database,
             const std::vector<index_t>& component_with_gas,
             units::UnitsSet units_set = {}
             );
 
     VariablesInterfaceImpl(
             std::shared_ptr<UnsaturatedVariables> vars
             );
 
     MainVariable& liquid_saturation() {
         return m_vars->get_liquid_saturation();
     }
     const MainVariable& liquid_saturation() const {
         return m_vars->get_liquid_saturation();
     }
 
     MainVariable& aqueous_concentration(index_t component) {
         return m_vars->get_aqueous_concentration(component);
     }
     const MainVariable& aqueous_concentration(index_t component) const {
         return m_vars->get_aqueous_concentration(component);
     }
 
     MainVariable& partial_pressure(index_t component) {
         if (not m_vars->component_has_gas(component)) {
             throw std::logic_error("Component has no associated gas");
         }
         return m_vars->get_pressure_main_variables(component);
     }
     const MainVariable& partial_pressure(index_t component) const {
         return m_vars->get_pressure_main_variables(component);
     }
 
 
     MainVariable& solid_concentration(index_t component) {
         return m_vars->get_solid_concentration(component);
     }
     const MainVariable& solid_concentration(index_t component) const {
         return m_vars->get_solid_concentration(component);
     }
 
     SecondaryTransientVariable& porosity() {
         return m_vars->get_porosity();
     }
     const SecondaryTransientVariable& porosity() const {
         return m_vars->get_porosity();
     }
 
 
     SecondaryTransientVariable& water_aqueous_concentration() {
         return m_vars->get_water_aqueous_concentration();
     }
     const SecondaryTransientVariable& water_aqueous_concentration() const {
         return m_vars->get_water_aqueous_concentration();
     }
 
     SecondaryVariable& liquid_diffusivity() {
         return m_vars->get_liquid_diffusivity();
     }
     const SecondaryVariable& liquid_diffusivity() const {
         return m_vars->get_liquid_diffusivity();
     }
 
     SecondaryVariable& liquid_permeability() {
         return m_vars->get_liquid_permeability();
     }
     const SecondaryVariable& liquid_permeability() const {
         return m_vars->get_liquid_permeability();
     }
 
 
     scalar_t& binary_gas_diffusivity(index_t component) {
         return m_vars->get_binary_gas_diffusivity(component);
     }
     scalar_t binary_gas_diffusivity(index_t component) const {
         return m_vars->get_binary_gas_diffusivity(component);
     }
 
     SecondaryVariable& resistance_gas_diffusivity() {
         return m_vars->get_resistance_gas_diffusivity();
     }
     const SecondaryVariable& resistance_gas_diffusivity() const {
         return m_vars->get_resistance_gas_diffusivity();
     }
 
     SecondaryVariable& advection_flux() {
         return m_vars->get_advection_flux();
     }
     const SecondaryVariable& advection_flux() const {
         return m_vars->get_advection_flux();
     }
 
 
     VariablesValidity check_variables();
     bool check_bounds(const Vector& to_check, scalar_t lower, scalar_t upper);
 
     VariablesValidity check_saturation();
     VariablesValidity check_relative();
     VariablesValidity check_porosity();
     VariablesValidity check_transport_param();
 };
 
 
 // =================================
 //
 //         Implementation
 //
 // =================================
 
 VariablesInterface::VariablesInterface(
         mesh::Mesh1DPtr the_mesh,
         database::RawDatabasePtr the_database,
         const std::vector<index_t>& component_with_gas):
     m_impl(make_unique<VariablesInterfaceImpl>(
                the_mesh, the_database, component_with_gas))
 {}
 
 
 VariablesInterface::VariablesInterface(
         mesh::Mesh1DPtr the_mesh,
         database::RawDatabasePtr the_database,
         const std::vector<index_t>& component_with_gas,
         const units::UnitsSet& units_set
         ):
     m_impl(make_unique<VariablesInterfaceImpl>(
                the_mesh, the_database, component_with_gas, units_set))
 {}
 
 
 VariablesInterface::VariablesInterface(
         std::shared_ptr<UnsaturatedVariables> vars
         ):
     m_impl(make_unique<VariablesInterfaceImpl>(vars))
 {}
 
 VariablesInterface::VariablesInterfaceImpl::VariablesInterfaceImpl(
         mesh::Mesh1DPtr the_mesh,
         database::RawDatabasePtr the_database,
         const std::vector<index_t>& component_with_gas,
         units::UnitsSet units_set
         ):
     m_vars(nullptr),
     m_mesh(the_mesh),
     m_database(the_database)
 {
     std::vector<bool> comp_has_gas(the_database->nb_component(), false);
     for (auto component: component_with_gas) {
         comp_has_gas[component] = true;
     }
     m_vars = std::make_shared<UnsaturatedVariables>(
-                the_mesh, the_database, comp_has_gas, units_set.length);
+                the_mesh, the_database, comp_has_gas, units_set);
 
 }
 
 
 VariablesInterface::VariablesInterfaceImpl::VariablesInterfaceImpl(
         std::shared_ptr<UnsaturatedVariables> vars
         ):
     m_vars(vars),
     m_mesh(m_vars->get_mesh()),
     m_database(m_vars->get_database())
 {}
 
 VariablesInterface::~VariablesInterface() = default;
 
 std::shared_ptr<VariablesBase> VariablesInterface::get_base_variables()
 {
     return std::static_pointer_cast<VariablesBase>(m_impl->m_vars);
 }
 
 std::shared_ptr<UnsaturatedVariables> VariablesInterface::get_variables()
 {
     return m_impl->m_vars;
 }
 
 
 UnsaturatedVariables* VariablesInterface::get_raw_variables()
 {
     return m_impl->m_vars.get();
 }
 
 // Saturation
 
 const Vector& VariablesInterface::get_liquid_saturation() const
 {
     return m_impl->liquid_saturation().variable;
 }
 
 void VariablesInterface::set_liquid_saturation(index_t node, scalar_t value)
 {
     m_impl->liquid_saturation().variable(node) = value;
     m_impl->liquid_saturation().predictor(node) = value;
 }
 
 void VariablesInterface::set_liquid_saturation(const Vector& value)
 {
     m_impl->liquid_saturation().variable = value;
     m_impl->liquid_saturation().predictor = value;
 }
 
 void VariablesInterface::set_liquid_saturation(scalar_t value)
 {
     m_impl->liquid_saturation().variable.setConstant(value);
     m_impl->liquid_saturation().predictor.setConstant(value);
 }
 
 // Aqueous component concentration
 
 
 const Vector& VariablesInterface::get_aqueous_concentration(index_t component) const
 {
     return m_impl->aqueous_concentration(component).variable;
 }
 
 void VariablesInterface::set_aqueous_concentration(index_t component, index_t node, scalar_t value)
 {
     MainVariable& aq_conc = m_impl->aqueous_concentration(component);
     aq_conc.variable(node) = value;
     aq_conc.predictor(node) = value;
 }
 
 void VariablesInterface::set_aqueous_concentration(index_t component, const Vector& value)
 {
     MainVariable& aq_conc = m_impl->aqueous_concentration(component);
     aq_conc.variable = value;
     aq_conc.predictor = value;
 }
 
 void VariablesInterface::set_aqueous_concentration(index_t component, scalar_t value)
 {
     MainVariable& aq_conc = m_impl->aqueous_concentration(component);
     aq_conc.variable.setConstant(value);
     aq_conc.predictor.setConstant(value);
 }
 
 
 // Partial pressure
 
 
 const Vector& VariablesInterface::get_partial_pressure(index_t component) const
 {
     return m_impl->partial_pressure(component).variable;
 }
 
 void VariablesInterface::set_partial_pressure(index_t component, index_t node, scalar_t value)
 {
     MainVariable& comp_pressure = m_impl->partial_pressure(component);
     comp_pressure.variable(node) = value;
     comp_pressure.predictor(node) = value;
 }
 
 void VariablesInterface::set_partial_pressure(index_t component, const Vector& value)
 {
     MainVariable& comp_pressure = m_impl->partial_pressure(component);
     comp_pressure.variable = value;
     comp_pressure.predictor = value;
 }
 
 void VariablesInterface::set_partial_pressure(index_t component, scalar_t value)
 {
     MainVariable& comp_pressure = m_impl->partial_pressure(component);
     comp_pressure.variable.setConstant(value);
     comp_pressure.predictor.setConstant(value);
 }
 
 
 // Solid concentration
 
 
 const Vector& VariablesInterface::get_solid_concentration(index_t component) const
 {
     return m_impl->solid_concentration(component).variable;
 }
 
 void VariablesInterface::set_solid_concentration(index_t component, index_t node, scalar_t value)
 {
     MainVariable& solid_conc = m_impl->solid_concentration(component);
     solid_conc.variable(node) = value;
     solid_conc.predictor(node) = value;
 }
 
 void VariablesInterface::set_solid_concentration(index_t component, const Vector& value)
 {
     MainVariable& solid_conc = m_impl->solid_concentration(component);
     solid_conc.variable = value;
     solid_conc.predictor = value;
 }
 
 void VariablesInterface::set_solid_concentration(index_t component, scalar_t value)
 {
     MainVariable& solid_conc = m_impl->solid_concentration(component);
     solid_conc.variable.setConstant(value);
     solid_conc.predictor.setConstant(value);
 }
 
 
 // Porosity
 
 const Vector& VariablesInterface::get_porosity() const
 {
     return m_impl->porosity().variable;
 }
 
 void VariablesInterface::set_porosity(index_t node, scalar_t value)
 {
     m_impl->porosity().variable(node) = value;
     m_impl->porosity().predictor(node) = value;
 }
 
 void VariablesInterface::set_porosity(const Vector& value)
 {
     m_impl->porosity().variable = value;
     m_impl->porosity().predictor = value;
 }
 
 void VariablesInterface::set_porosity(scalar_t value)
 {
     m_impl->porosity().variable.setConstant(value);
     m_impl->porosity().predictor.setConstant(value);
 }
 
 // Water aqueous concentration
 
 
 const Vector& VariablesInterface::get_water_aqueous_concentration() const
 {
     return m_impl->water_aqueous_concentration().variable;
 }
 
 void VariablesInterface::set_water_aqueous_concentration(index_t node, scalar_t value)
 {
     m_impl->water_aqueous_concentration().variable(node) = value;
     m_impl->water_aqueous_concentration().predictor(node) = value;
 }
 
 void VariablesInterface::set_water_aqueous_concentration(const Vector& value)
 {
     m_impl->water_aqueous_concentration().variable = value;
     m_impl->water_aqueous_concentration().predictor = value;
 }
 
 void VariablesInterface::set_water_aqueous_concentration(scalar_t value)
 {
     m_impl->water_aqueous_concentration().variable.setConstant(value);
     m_impl->water_aqueous_concentration().predictor.setConstant(value);
 }
 
 
 // Liquid diffusion coefficient
 
 const Vector& VariablesInterface::get_liquid_diffusivity() const
 {
     return m_impl->liquid_diffusivity().variable;
 }
 
 void VariablesInterface::set_liquid_diffusivity(index_t node, scalar_t value)
 {
     m_impl->liquid_diffusivity().variable(node) = value;
 }
 
 void VariablesInterface::set_liquid_diffusivity(const Vector& value)
 {
     m_impl->liquid_diffusivity().variable = value;
 }
 
 void VariablesInterface::set_liquid_diffusivity(scalar_t value)
 {
     m_impl->liquid_diffusivity().variable.setConstant(value);
 }
 
 // Liquid permeability
 
 const Vector& VariablesInterface::get_liquid_permeability() const
 {
     return m_impl->liquid_permeability().variable;
 }
 
 void VariablesInterface::set_liquid_permeability(index_t node, scalar_t value)
 {
     m_impl->liquid_permeability().variable(node) = value;
 }
 
 void VariablesInterface::set_liquid_permeability(const Vector& value)
 {
     m_impl->liquid_permeability().variable = value;
 }
 
 void VariablesInterface::set_liquid_permeability(scalar_t value)
 {
     m_impl->liquid_permeability().variable.setConstant(value);
 }
 
 // Gas diffusivity
 
 scalar_t
 VariablesInterface::get_binary_gas_diffusivity(
         index_t component
         ) const
 {
     return m_impl->binary_gas_diffusivity(component);
 }
 
 void
 VariablesInterface::set_binary_gas_diffusivity(
         index_t component,
         scalar_t value
         )
 {
     m_impl->binary_gas_diffusivity(component) = value;
 }
 
 
 const Vector&
 VariablesInterface::get_resistance_gas_diffusivity() const
 {
     return m_impl->resistance_gas_diffusivity().variable;
 }
 
 void
 VariablesInterface::set_resistance_gas_diffusivity(
         index_t node,
         scalar_t value
         )
 {
     m_impl->resistance_gas_diffusivity().variable(node) = value;
 }
 void
 VariablesInterface::set_resistance_gas_diffusivity(
         const Vector& value)
 {
     m_impl->resistance_gas_diffusivity().variable = value;
 }
 void
 VariablesInterface::set_resistance_gas_diffusivity(
         scalar_t value)
 {
     m_impl->resistance_gas_diffusivity().variable.setConstant(value);
 }
 
 const Vector& VariablesInterface::get_advection_flux() const
 {
     return m_impl->advection_flux().variable;
 }
 
 void VariablesInterface::set_advection_flux(index_t node, scalar_t value)
 {
     m_impl->advection_flux().variable(node) = value;
 }
 
 void VariablesInterface::set_advection_flux(const Vector& value)
 {
     m_impl->advection_flux().variable = value;
 }
 
 void VariablesInterface::set_advection_flux(scalar_t value)
 {
     m_impl->advection_flux().variable.setConstant(value);
 }
 
 // User Models
 
 const Vector& VariablesInterface::get_capillary_pressure() const
 {
     return m_impl->m_vars->get_capillary_pressure().variable;
 }
 
 void VariablesInterface::set_capillary_pressure_model(
         user_model_saturation_f capillary_pressure_model
         )
 {
     m_impl->m_vars->set_capillary_pressure_model(capillary_pressure_model);
     const Vector& sat = get_liquid_saturation();
     Vector& cap_pressure = m_impl->m_vars->get_capillary_pressure().variable;
     for (auto node: m_impl->m_mesh->range_nodes())
     {
         cap_pressure(node) = capillary_pressure_model(node, sat(node));
     }
 }
 
 const Vector& VariablesInterface::get_vapor_pressure() const
 {
     return m_impl->m_vars->get_pressure_main_variables(0).variable;
 }
 void VariablesInterface::set_vapor_pressure_model(
         user_model_saturation_f vapor_pressure_model
         )
 {
     m_impl->m_vars->set_vapor_pressure_model(vapor_pressure_model);
 }
 
 
 const Vector& VariablesInterface::get_relative_liquid_diffusivity() const
 {
     return m_impl->m_vars->get_relative_liquid_diffusivity().variable;
 }
 
 void VariablesInterface::set_relative_liquid_diffusivity_model(
         user_model_saturation_f relative_liquid_diffusivity_model
         )
 {
     m_impl->m_vars->set_relative_liquid_diffusivity_model(relative_liquid_diffusivity_model);
     const Vector& sat = get_liquid_saturation();
     Vector& rel_diff = m_impl->m_vars->get_relative_liquid_diffusivity().variable;
     for (auto node: m_impl->m_mesh->range_nodes())
     {
         rel_diff(node) = relative_liquid_diffusivity_model(node, sat(node));
     }
 }
 
 const Vector& VariablesInterface::get_relative_liquid_permeability() const
 {
     return m_impl->m_vars->get_relative_liquid_permeability().variable;
 }
 
 void VariablesInterface::set_relative_liquid_permeability_model(
         user_model_saturation_f relative_liquid_permeability_model
         )
 {
     m_impl->m_vars->set_relative_liquid_permeability_model(relative_liquid_permeability_model);
     const Vector& sat = get_liquid_saturation();
     Vector& rel_diff = m_impl->m_vars->get_relative_liquid_permeability().variable;
     for (auto node: m_impl->m_mesh->range_nodes())
     {
         rel_diff(node) = relative_liquid_permeability_model(node, sat(node));
     }
 }
 
 const Vector& VariablesInterface::get_relative_gas_diffusivity() const
 {
     return m_impl->m_vars->get_relative_gas_diffusivity().variable;
 }
 
 void VariablesInterface::set_relative_gas_diffusivity_model(
         user_model_saturation_f relative_gas_diffusivity_model
         )
 {
     m_impl->m_vars->set_relative_gas_diffusivity_model(relative_gas_diffusivity_model);
     const Vector& sat = get_liquid_saturation();
     Vector& rel_diff = m_impl->m_vars->get_relative_gas_diffusivity().variable;
     for (auto node: m_impl->m_mesh->range_nodes())
     {
         rel_diff(node) = relative_gas_diffusivity_model(node, sat(node));
     }
 }
 
 
 // chemistry
 void VariablesInterface::initialize_variables(index_t node,
                           const AdimensionalSystemSolutionExtractor& extractor
                           )
 {
     set_porosity(node, extractor.porosity());
     scalar_t saturation = extractor.saturation_water();
     set_liquid_saturation(node, extractor.saturation_water());
     scalar_t rho_l = extractor.density_water();
 
     set_water_aqueous_concentration(node,
             rho_l*extractor.total_aqueous_concentration(0));
     set_solid_concentration(0, node,
             extractor.total_solid_concentration(0));
     if (m_impl->m_vars->component_has_gas(0))
     {
         set_partial_pressure(0, node,
                              m_impl->m_vars->get_vapor_pressure_model()(node, saturation)
                              );
     }
 
     for (index_t component: m_impl->m_database->range_aqueous_component())
     {
         set_aqueous_concentration(component, node,
                 rho_l*extractor.total_aqueous_concentration(component));
         set_solid_concentration(component, node,
                 extractor.total_solid_concentration(component));
         if (m_impl->m_vars->component_has_gas(component))
         {
             set_partial_pressure(component, node,
                                  extractor.fugacity_gas(m_impl->m_vars->get_id_gas(component))
                                  *m_impl->m_vars->get_total_pressure());
         }
     }
 
     m_impl->m_vars->set_adim_solution(node,
                                       std::move(extractor.get_solution()));
 }
 
 VariablesValidity VariablesInterface::check_variables()
 {
     return m_impl->check_variables();
 }
 
 
 bool VariablesInterface::VariablesInterfaceImpl::check_bounds(
         const Vector& to_check,
         scalar_t lower,
         scalar_t upper
         )
 {
     bool flag = true;
     for (index_t ind=0; ind< to_check.rows(); ++ind)
     {
         const scalar_t& value = to_check(ind);
         if (value < lower or value > upper) {
             flag = false;
             break;
         }
     }
     return flag;
 }
 
 VariablesValidity max_error_level(VariablesValidity err1, VariablesValidity err2)
 {
     return (err1>err2)?err1:err2;
 }
 
 VariablesValidity VariablesInterface::VariablesInterfaceImpl::check_variables()
 {
 
     VariablesValidity return_code = VariablesValidity::good;
 
     auto tmpret = check_saturation();
     return_code = max_error_level(return_code, tmpret);
 
     tmpret = check_porosity();
     return_code = max_error_level(return_code, tmpret);
 
     tmpret = check_relative();
     return_code = max_error_level(return_code, tmpret);
 
     return return_code;
 }
 
 static bool SPECMICP_DLL_LOCAL saturation_check(const Vector& x)
 {
     return (is_bounded(x, 0.0, 1.0) && (any(x) > 0.0));
 }
 
 VariablesValidity
 VariablesInterface::VariablesInterfaceImpl::check_saturation()
 {
     using namespace vector_checker;
     const auto& x = m_vars->get_liquid_saturation().variable;
 
     if (not saturation_check(x))
     {
         ERROR << "Saturation vector is incorrect !";
         return VariablesValidity::error;
     }
     return VariablesValidity::good;
 }
 
 static bool SPECMICP_DLL_LOCAL porosity_check(const Vector& x)
 {
     return (is_bounded(x, 0.0, 1.0) && (any(x) > 0.0));
 }
 
 VariablesValidity
 VariablesInterface::VariablesInterfaceImpl::check_porosity()
 {
     using namespace vector_checker;
     const auto& x = m_vars->get_porosity().variable;
 
     if (not porosity_check(x)) {
         ERROR << "Porosity is incorrect";
         return VariablesValidity::error;
     }
     return VariablesValidity::good;
 }
 
 static bool SPECMICP_DLL_LOCAL relative_variable_check(const Vector& x)
 {
     return (is_bounded(x, 0.0, 1.0) && (any(x) > 0.0));
 }
 
 VariablesValidity
 VariablesInterface::VariablesInterfaceImpl::check_relative()
 {
     auto is_valid = VariablesValidity::good;
 
     if (not relative_variable_check(
             m_vars->get_relative_liquid_diffusivity().variable)) {
         ERROR << "Relative liquid diffusivity is incorrect";
         is_valid = VariablesValidity::error;
     }
     if (not relative_variable_check(
             m_vars->get_relative_liquid_permeability().variable)) {
         ERROR << "Relative liquid permeability is incorrect";
         is_valid = VariablesValidity::error;
     }
     if (not relative_variable_check(
             m_vars->get_relative_gas_diffusivity().variable)) {
         ERROR << "Relative gas diffusivity is incorrect";
         is_valid = VariablesValidity::error;
     }
     if (not relative_variable_check(
             m_vars->get_resistance_gas_diffusivity().variable)) {
         ERROR << "Resistance gas diffusivity is incorrect";
         is_valid = VariablesValidity::error;
     }
 
     return is_valid;
 }
 
 static bool SPECMICP_DLL_LOCAL transport_parameters_check(const Vector& x)
 {
     return (is_lower_bounded(x, 0.0) && (any(x) > 0.0));
 }
 
 VariablesValidity
 VariablesInterface::VariablesInterfaceImpl::check_transport_param()
 {
     auto is_valid = VariablesValidity::good;
 
     if (not transport_parameters_check(
             m_vars->get_liquid_diffusivity().variable)) {
         ERROR << "Liquid diffusivity is incorrect";
         is_valid = VariablesValidity::error;
     }
     if (not transport_parameters_check(
             m_vars->get_liquid_permeability().variable)) {
         ERROR << "Liquid permeability is incorrect";
         is_valid = VariablesValidity::error;
     }
 
     return is_valid;
 }
 
 } //end namespace unsaturated
 } //end namespace systems
 } //end namespace reactmicp
 } //end namespace specmicp
diff --git a/src/reactmicp/systems/unsaturated/variables_sub.hpp b/src/reactmicp/systems/unsaturated/variables_sub.hpp
index 7558438..abafd8d 100644
--- a/src/reactmicp/systems/unsaturated/variables_sub.hpp
+++ b/src/reactmicp/systems/unsaturated/variables_sub.hpp
@@ -1,269 +1,268 @@
 /* =============================================================================
 
  Copyright (c) 2014 - 2016
  F. Georget <fabieng@princeton.edu> Princeton University
  All rights reserved.
 
 
 Redistribution and use in source and binary forms, with or without modification,
 are permitted provided that the following conditions are met:
 
     1. Redistributions of source code must retain the above copyright notice,
     this list of conditions and the following disclaimer.
 
     2. Redistributions in binary form must reproduce the above copyright notice,
     this list of conditions and the following disclaimer in the documentation
     and/or other materials provided with the distribution.
 
     3. Neither the name of the copyright holder nor the names of its
     contributors may be used to endorse or promote products derived from this
     software without specific prior written permission.
 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
 ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
 ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. *
 
 ============================================================================= */
 
 //! \file unsaturated/variables_sub.hpp
 //! \brief base structures for the variables
 
 #ifndef SPECMICP_REACTMICP_UNSATURATED_VARIABLESSUB_HPP
 #define SPECMICP_REACTMICP_UNSATURATED_VARIABLESSUB_HPP
 
 #include "specmicp_common/types.hpp"
 #include "specmicp/adimensional/adimensional_system_solution.hpp"
 
 #include "specmicp_common/compat.hpp"
 
 #include "specmicp_common/physics/constants.hpp"
 #include "specmicp_common/physics/units.hpp"
 
 #include <vector>
 #include <memory>
 
 namespace specmicp {
 namespace reactmicp {
 namespace systems {
 namespace unsaturated {
 
 //! \brief A simple variable
 struct BaseVariable
 {
 
     Vector variable;      //!< The variable
 
     //! \brief Return value of the variable
     scalar_t& operator() (index_t node) {
         return variable(node);
     }
 
     //! \brief Return value of the variable
     scalar_t operator() (index_t node) const {
         return variable(node);
     }
 
     //! \brief Return the size of the vector
     index_t size() const {
         return variable.rows();
     }
 
     //! \brief The variables are all set to 'value'
     void set_constant(scalar_t value) {
         variable.setConstant(value);
     }
 
     //! \brief Set the variables to zero
     void set_zero() {
         variable.setZero();
     }
 
     //! \brief Build a variable of size 'size'
     //!
     //! The varaibles are initialized to zero
     BaseVariable(index_t size):
         variable(Vector::Zero(size))
     {}
 };
 
 //! \brief A transient variable
 //!
 //! These type of variable store the rate of change of the variable (velocity)
 //! and the value at the beginning of the timestep (predictor)
 struct BaseTransientVariable:
         public BaseVariable
 {
     Vector velocity;      //!< Rate of change of the variable
     Vector predictor;     //!< Predictor, value before first iteration
 
     BaseTransientVariable(index_t size):
         BaseVariable(size),
         velocity(Vector::Zero(size)),
         predictor(Vector::Zero(size))
     {}
 
     //! \brief Update the variable
     void update(scalar_t dt) {
         variable = predictor + dt*velocity;
     }
 };
 
 //! \brief A main variable, used to solve the governing equations
 //!
 //! These variables also store the chemistry and transport rates for the
 //!  coupling
 struct MainVariable:
         public BaseTransientVariable
 {
     //! \brief Transport fluxes of the corresponding governing equations
     Vector transport_fluxes;
     //! \brief Chemical exchange rate
     Vector chemistry_rate;
 
     MainVariable(index_t size):
         BaseTransientVariable(size),
         transport_fluxes(Vector::Zero(size)),
         chemistry_rate(Vector::Zero(size))
     {}
 
     //! \brief Reset the variables
     //!
     //! This function is called if the solver failed and must be restarted
     void reset() {
         transport_fluxes.setZero();
         chemistry_rate.setZero();
         velocity.setZero();
         variable = predictor;
     }
 };
 
 //! \brief List of main variables
 struct ListMainVariable
 {
     // Variables are stored as pointer so we can skip components
     // that do not exist
     // But still access their variables through their index
     using value_type = std::unique_ptr<MainVariable>;
     using vector_type = std::vector<value_type>;
     using iterator = vector_type::iterator;
     using const_iterator = vector_type::const_iterator;
 
     vector_type variables;
 
     MainVariable& water() {
         return *variables[0];
     }
 
     MainVariable& aqueous_component(index_t aq_component) {
         specmicp_assert(variables[aq_component] != nullptr);
         return *variables[aq_component];
     }
 
     MainVariable& component(index_t component) {
         specmicp_assert(variables[component] != nullptr);
         return *variables[component];
     }
 
     MainVariable* operator() (index_t component) {
         return variables[component].get();
     }
 
     // Implementations of the following methods
     // is in variables.cpp
 
     //! \brief Initialize the variables
     ListMainVariable(index_t nb_vars, index_t nb_nodes);
     //! \brief Initialize some variables
     //!
     //! Only initialize values given by 'to_init'
     ListMainVariable(index_t nb_vars, index_t nb_nodes, std::vector<bool> to_init);
 
     //! \brief Reset the variables
     void reset();
 
     //! \brief Hard reset, erase all info for a component;
     void hard_reset(index_t component, index_t nb_nodes);
 };
 
 //! \brief A secondary variable
 struct SecondaryVariable:
         public BaseVariable
 {
     SecondaryVariable(index_t size):
         BaseVariable(size)
     {}
 
 };
 
 //! \brief A secondary transient variable
 //!
 //! e.g. the porosity
 struct SecondaryTransientVariable:
         public BaseTransientVariable
 {
     SecondaryTransientVariable(index_t size):
         BaseTransientVariable(size)
     {}
 };
 
 //! \brief The chemistry solutions
 struct ChemistrySolutions
 {
     std::vector<AdimensionalSystemSolution> solutions;
 
     ChemistrySolutions(index_t size):
         solutions(size)
     {}
 
     //! \brief Return a solution
     AdimensionalSystemSolution& solution(index_t node) {
         return solutions[node];
     }
 
     //! \brief Return a const solution
     const AdimensionalSystemSolution& solution(index_t node) const {
         return solutions[node];
     }
 
     //! \brief Set a solution
     void update_solution(
             index_t node,
             const AdimensionalSystemSolution& new_solution
             ) {
         solutions[node] = new_solution;
     }
 };
 
 //! \brief List of constants used in the problems
 struct ConstantBox
 {
     //! \brief R*T
     //!
     //! R : ideal gas constant
     //! T : temperature
     scalar_t rt {constants::gas_constant*units::celsius(25.0)};
     //! \brief  Viscosity of liquid water
     scalar_t viscosity_liquid_water {constants::water_viscosity};
     //! \brief Total pressure
     scalar_t total_pressure {1.01325e5};
 
-    //! \brief Scale the constants
-    void scale(units::LengthUnit length_unit);
+    void scale(const units::UnitsSet& units_set);
 };
 
 
 
 
 } //end namespace unsaturated
 } //end namespace systems
 } //end namespace reactmicp
 } //end namespace specmicp
 
 #endif // SPECMICP_REACTMICP_UNSATURATED_VARIABLESSUB_HPP