Repository: fast-aircraft-design/FAST-OAD (language: python)

Issue:
Sub models
This PR implements the idea in #255, that is to provide a way of replacing only a part of a FAST-OAD module.

Usage is described [here](https://fast-oad.readthedocs.io/en/sub-models-complete/documentation/custom_modules/add_submodels.html).

I aimed at having the `RegisterSubModel` class sharing the same base as the existing `RegisterOpenMDAOSystem` and `RegisterPropulsion` classes, so all these classes have been reworked somehow.

Also, this PR has been the occasion to update some Sphinx settings.


Discussion / hints:
# [Codecov](https://codecov.io/gh/fast-aircraft-design/FAST-OAD/pull/379?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=fast-aircraft-design) Report
> Merging [#379](https://codecov.io/gh/fast-aircraft-design/FAST-OAD/pull/379?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=fast-aircraft-design) (f634bda) into [master](https://codecov.io/gh/fast-aircraft-design/FAST-OAD/commit/6ccd0ae973555f2ed846250fef6271d65bf21e41?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=fast-aircraft-design) (6ccd0ae) will **increase** coverage by `0.07%`.
> The diff coverage is `100.00%`.

[![Impacted file tree graph](https://codecov.io/gh/fast-aircraft-design/FAST-OAD/pull/379/graphs/tree.svg?width=650&height=150&src=pr&token=1J1R5DT61A&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=fast-aircraft-design)](https://codecov.io/gh/fast-aircraft-design/FAST-OAD/pull/379?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=fast-aircraft-design)

```diff
@@            Coverage Diff             @@
##           master     #379      +/-   ##
==========================================
+ Coverage   91.44%   91.52%   +0.07%     
==========================================
  Files         192      192              
  Lines        7761     7831      +70     
  Branches      734      740       +6     
==========================================
+ Hits         7097     7167      +70     
  Misses        523      523              
  Partials      141      141              
```


| [Impacted Files](https://codecov.io/gh/fast-aircraft-design/FAST-OAD/pull/379?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=fast-aircraft-design) | Coverage Δ | |
|---|---|---|
| [src/fastoad/api.py](https://codecov.io/gh/fast-aircraft-design/FAST-OAD/pull/379/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=fast-aircraft-design#diff-c3JjL2Zhc3RvYWQvYXBpLnB5) | `0.00% <ø> (ø)` | |
| [src/fastoad/io/configuration/configuration.py](https://codecov.io/gh/fast-aircraft-design/FAST-OAD/pull/379/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=fast-aircraft-design#diff-c3JjL2Zhc3RvYWQvaW8vY29uZmlndXJhdGlvbi9jb25maWd1cmF0aW9uLnB5) | `92.19% <100.00%> (+0.33%)` | :arrow_up: |
| [src/fastoad/module\_management/exceptions.py](https://codecov.io/gh/fast-aircraft-design/FAST-OAD/pull/379/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=fast-aircraft-design#diff-c3JjL2Zhc3RvYWQvbW9kdWxlX21hbmFnZW1lbnQvZXhjZXB0aW9ucy5weQ==) | `100.00% <100.00%> (ø)` | |
| [src/fastoad/module\_management/service\_registry.py](https://codecov.io/gh/fast-aircraft-design/FAST-OAD/pull/379/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=fast-aircraft-design#diff-c3JjL2Zhc3RvYWQvbW9kdWxlX21hbmFnZW1lbnQvc2VydmljZV9yZWdpc3RyeS5weQ==) | `100.00% <100.00%> (ø)` | |

------

[Continue to review full report at Codecov](https://codecov.io/gh/fast-aircraft-design/FAST-OAD/pull/379?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=fast-aircraft-design).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=fast-aircraft-design)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/fast-aircraft-design/FAST-OAD/pull/379?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=fast-aircraft-design). Last update [6ccd0ae...f634bda](https://codecov.io/gh/fast-aircraft-design/FAST-OAD/pull/379?src=pr&el=lastupdated&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=fast-aircraft-design). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=fast-aircraft-design).


Fix patch:
diff --git a/src/fastoad/models/aerodynamics/__init__.py b/src/fastoad/models/aerodynamics/__init__.py
index b75be4c43..385ae0611 100644
--- a/src/fastoad/models/aerodynamics/__init__.py
+++ b/src/fastoad/models/aerodynamics/__init__.py
@@ -1,5 +1,5 @@
 #  This file is part of FAST-OAD : A framework for rapid Overall Aircraft Design
-#  Copyright (C) 2020  ONERA & ISAE-SUPAERO
+#  Copyright (C) 2021 ONERA & ISAE-SUPAERO
 #  FAST is free software: you can redistribute it and/or modify
 #  it under the terms of the GNU General Public License as published by
 #  the Free Software Foundation, either version 3 of the License, or
@@ -10,4 +10,3 @@
 #  GNU General Public License for more details.
 #  You should have received a copy of the GNU General Public License
 #  along with this program.  If not, see <https://www.gnu.org/licenses/>.
-from .aerodynamics import Aerodynamics
diff --git a/src/fastoad/models/aerodynamics/aerodynamics.py b/src/fastoad/models/aerodynamics/aerodynamics.py
deleted file mode 100644
index ce49b9f05..000000000
--- a/src/fastoad/models/aerodynamics/aerodynamics.py
+++ /dev/null
@@ -1,34 +0,0 @@
-"""
-    FAST - Copyright (c) 2016 ONERA ISAE
-"""
-
-#  This file is part of FAST-OAD : A framework for rapid Overall Aircraft Design
-#  Copyright (C) 2020  ONERA & ISAE-SUPAERO
-#  FAST is free software: you can redistribute it and/or modify
-#  it under the terms of the GNU General Public License as published by
-#  the Free Software Foundation, either version 3 of the License, or
-#  (at your option) any later version.
-#  This program is distributed in the hope that it will be useful,
-#  but WITHOUT ANY WARRANTY; without even the implied warranty of
-#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-#  GNU General Public License for more details.
-#  You should have received a copy of the GNU General Public License
-#  along with this program.  If not, see <https://www.gnu.org/licenses/>.
-
-from openmdao.api import Group
-
-from fastoad.models.aerodynamics.aerodynamics_high_speed import AerodynamicsHighSpeed
-from fastoad.models.aerodynamics.aerodynamics_landing import AerodynamicsLanding
-from fastoad.models.aerodynamics.aerodynamics_low_speed import AerodynamicsLowSpeed
-
-
-class Aerodynamics(Group):
-    def setup(self):
-        # Compute the low speed aero (Cl alpha at takeoff and Cl0)
-        self.add_subsystem("aero_low", AerodynamicsLowSpeed(), promotes=["*"])
-
-        # Compute CL max landing
-        self.add_subsystem("aero_landing", AerodynamicsLanding(), promotes=["*"])
-
-        # Compute cruise characteristics
-        self.add_subsystem("aero_high", AerodynamicsHighSpeed(), promotes=["*"])
diff --git a/src/fastoad/models/aerodynamics/aerodynamics_high_speed.py b/src/fastoad/models/aerodynamics/aerodynamics_high_speed.py
index d36fe19cb..e9f30f27e 100644
--- a/src/fastoad/models/aerodynamics/aerodynamics_high_speed.py
+++ b/src/fastoad/models/aerodynamics/aerodynamics_high_speed.py
@@ -1,6 +1,4 @@
-"""
-    FAST - Copyright (c) 2016 ONERA ISAE
-"""
+"""Computation of aerodynamic polar in cruise conditions."""
 #  This file is part of FAST-OAD : A framework for rapid Overall Aircraft Design
 #  Copyright (C) 2021 ONERA & ISAE-SUPAERO
 #  FAST is free software: you can redistribute it and/or modify
@@ -17,19 +15,17 @@
 import openmdao.api as om
 
 from fastoad.module_management.constants import ModelDomain
-from fastoad.module_management.service_registry import RegisterOpenMDAOSystem
-from .components.cd0_fuselage import Cd0Fuselage
-from .components.cd0_ht import Cd0HorizontalTail
-from .components.cd0_nacelle_pylons import Cd0NacelleAndPylons
-from .components.cd0_total import Cd0Total
-from .components.cd0_vt import Cd0VerticalTail
-from .components.cd0_wing import Cd0Wing
-from .components.cd_compressibility import CdCompressibility
-from .components.cd_trim import CdTrim
-from .components.compute_polar import ComputePolar
-from .components.compute_reynolds import ComputeReynolds
-from .components.initialize_cl import InitializeClPolar
-from .components.oswald import OswaldCoefficient
+from fastoad.module_management.service_registry import RegisterOpenMDAOSystem, RegisterSubmodel
+from .constants import (
+    SERVICE_CD0,
+    SERVICE_CD_COMPRESSIBILITY,
+    SERVICE_CD_TRIM,
+    SERVICE_INDUCED_DRAG_COEFFICIENT,
+    SERVICE_INITIALIZE_CL,
+    SERVICE_OSWALD_COEFFICIENT,
+    SERVICE_POLAR,
+    SERVICE_REYNOLDS_COEFFICIENT,
+)
 
 
 @RegisterOpenMDAOSystem("fastoad.aerodynamics.highspeed.legacy", domain=ModelDomain.AERODYNAMICS)
@@ -42,15 +38,29 @@ class AerodynamicsHighSpeed(om.Group):
     """
 
     def setup(self):
-        self.add_subsystem("compute_oswald_coeff", OswaldCoefficient(), promotes=["*"])
-        self.add_subsystem("comp_re", ComputeReynolds(), promotes=["*"])
-        self.add_subsystem("initialize_cl", InitializeClPolar(), promotes=["*"])
-        self.add_subsystem("cd0_wing", Cd0Wing(), promotes=["*"])
-        self.add_subsystem("cd0_fuselage", Cd0Fuselage(), promotes=["*"])
-        self.add_subsystem("cd0_ht", Cd0HorizontalTail(), promotes=["*"])
-        self.add_subsystem("cd0_vt", Cd0VerticalTail(), promotes=["*"])
-        self.add_subsystem("cd0_nac_pylons", Cd0NacelleAndPylons(), promotes=["*"])
-        self.add_subsystem("cd0_total", Cd0Total(), promotes=["*"])
-        self.add_subsystem("cd_comp", CdCompressibility(), promotes=["*"])
-        self.add_subsystem("cd_trim", CdTrim(), promotes=["*"])
-        self.add_subsystem("get_polar", ComputePolar(), promotes=["*"])
+        self.add_subsystem(
+            "compute_oswald_coeff",
+            RegisterSubmodel.get_submodel(SERVICE_OSWALD_COEFFICIENT),
+            promotes=["*"],
+        )
+        self.add_subsystem(
+            "compute_induced_drag_coeff",
+            RegisterSubmodel.get_submodel(SERVICE_INDUCED_DRAG_COEFFICIENT),
+            promotes=["*"],
+        )
+        self.add_subsystem(
+            "comp_re", RegisterSubmodel.get_submodel(SERVICE_REYNOLDS_COEFFICIENT), promotes=["*"]
+        )
+        self.add_subsystem(
+            "initialize_cl", RegisterSubmodel.get_submodel(SERVICE_INITIALIZE_CL), promotes=["*"]
+        )
+        self.add_subsystem("cd0_wing", RegisterSubmodel.get_submodel(SERVICE_CD0), promotes=["*"])
+        self.add_subsystem(
+            "cd_comp", RegisterSubmodel.get_submodel(SERVICE_CD_COMPRESSIBILITY), promotes=["*"]
+        )
+        self.add_subsystem(
+            "cd_trim", RegisterSubmodel.get_submodel(SERVICE_CD_TRIM), promotes=["*"]
+        )
+        self.add_subsystem(
+            "get_polar", RegisterSubmodel.get_submodel(SERVICE_POLAR), promotes=["*"]
+        )
diff --git a/src/fastoad/models/aerodynamics/aerodynamics_landing.py b/src/fastoad/models/aerodynamics/aerodynamics_landing.py
index 688091a2e..c22c87683 100644
--- a/src/fastoad/models/aerodynamics/aerodynamics_landing.py
+++ b/src/fastoad/models/aerodynamics/aerodynamics_landing.py
@@ -19,10 +19,14 @@
 
 from fastoad.model_base import Atmosphere
 from fastoad.module_management.constants import ModelDomain
-from fastoad.module_management.service_registry import RegisterOpenMDAOSystem
-from .components.compute_max_cl_landing import ComputeMaxClLanding
-from .components.high_lift_aero import ComputeDeltaHighLift
-from .external.xfoil import XfoilPolar
+from fastoad.module_management.service_registry import RegisterOpenMDAOSystem, RegisterSubmodel
+from .constants import (
+    SERVICE_HIGH_LIFT,
+    SERVICE_LANDING_MACH_REYNOLDS,
+    SERVICE_LANDING_MAX_CL,
+    SERVICE_LANDING_MAX_CL_CLEAN,
+    SERVICE_XFOIL,
+)
 from .external.xfoil.xfoil_polar import (
     OPTION_ALPHA_END,
     OPTION_ALPHA_START,
@@ -85,25 +89,46 @@ def initialize(self):
         )
 
     def setup(self):
-        self.add_subsystem("mach_reynolds", ComputeMachReynolds(), promotes=["*"])
+        self.add_subsystem(
+            "mach_reynolds",
+            RegisterSubmodel.get_submodel(SERVICE_LANDING_MACH_REYNOLDS),
+            promotes=["*"],
+        )
+
         if self.options["use_xfoil"]:
             start = self.options["xfoil_alpha_min"]
             end = self.options["xfoil_alpha_max"]
             iter_limit = self.options["xfoil_iter_limit"]
-            kwargs = {
+            xfoil_options = {
                 OPTION_ALPHA_START: start,
                 OPTION_ALPHA_END: end,
                 OPTION_ITER_LIMIT: iter_limit,
                 OPTION_XFOIL_EXE_PATH: self.options[OPTION_XFOIL_EXE_PATH],
             }
             self.add_subsystem(
-                "xfoil_run", XfoilPolar(**kwargs), promotes=["data:geometry:wing:thickness_ratio"],
+                "xfoil_run",
+                RegisterSubmodel.get_submodel(SERVICE_XFOIL, xfoil_options),
+                promotes=["data:geometry:wing:thickness_ratio"],
             )
-        self.add_subsystem("CL_2D_to_3D", Compute3DMaxCL(), promotes=["*"])
+
         self.add_subsystem(
-            "delta_cl_landing", ComputeDeltaHighLift(landing_flag=True), promotes=["*"]
+            "CL_2D_to_3D",
+            RegisterSubmodel.get_submodel(SERVICE_LANDING_MAX_CL_CLEAN),
+            promotes=["*"],
+        )
+
+        landing_flag_option = {"landing_flag": True}
+        self.add_subsystem(
+            "delta_cl_landing",
+            RegisterSubmodel.get_submodel(SERVICE_HIGH_LIFT, landing_flag_option),
+            promotes=["*"],
+        )
+
+        self.add_subsystem(
+            "compute_max_cl_landing",
+            RegisterSubmodel.get_submodel(SERVICE_LANDING_MAX_CL),
+            promotes=["*"],
         )
-        self.add_subsystem("compute_max_cl_landing", ComputeMaxClLanding(), promotes=["*"])
 
         if self.options["use_xfoil"]:
             self.connect("data:aerodynamics:aircraft:landing:mach", "xfoil_run.xfoil:mach")
@@ -113,6 +138,9 @@ def setup(self):
             )
 
 
+@RegisterSubmodel(
+    SERVICE_LANDING_MACH_REYNOLDS, "fastoad.submodel.aerodynamics.landing.mach_reynolds.legacy"
+)
 class ComputeMachReynolds(om.ExplicitComponent):
     """
     Mach and Reynolds computation
@@ -139,6 +167,9 @@ def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
         outputs["data:aerodynamics:wing:landing:reynolds"] = reynolds
 
 
+@RegisterSubmodel(
+    SERVICE_LANDING_MAX_CL_CLEAN, "fastoad.submodel.aerodynamics.landing.max_CL_clean.legacy"
+)
 class Compute3DMaxCL(om.ExplicitComponent):
     """
     Computes 3D max CL from 2D CL (XFOIL-computed) and sweep angle
diff --git a/src/fastoad/models/aerodynamics/aerodynamics_low_speed.py b/src/fastoad/models/aerodynamics/aerodynamics_low_speed.py
index cffa65988..d6c05bf22 100644
--- a/src/fastoad/models/aerodynamics/aerodynamics_low_speed.py
+++ b/src/fastoad/models/aerodynamics/aerodynamics_low_speed.py
@@ -1,7 +1,4 @@
-"""
-    FAST - Copyright (c) 2016 ONERA ISAE
-"""
-
+"""Computation of aerodynamic polar in low speed conditions."""
 #  This file is part of FAST-OAD : A framework for rapid Overall Aircraft Design
 #  Copyright (C) 2021 ONERA & ISAE-SUPAERO
 #  FAST is free software: you can redistribute it and/or modify
@@ -18,19 +15,18 @@
 import openmdao.api as om
 
 from fastoad.module_management.constants import ModelDomain
-from fastoad.module_management.service_registry import RegisterOpenMDAOSystem
-from .components.cd0_fuselage import Cd0Fuselage
-from .components.cd0_ht import Cd0HorizontalTail
-from .components.cd0_nacelle_pylons import Cd0NacelleAndPylons
-from .components.cd0_total import Cd0Total
-from .components.cd0_vt import Cd0VerticalTail
-from .components.cd0_wing import Cd0Wing
-from .components.cd_trim import CdTrim
-from .components.compute_low_speed_aero import ComputeAerodynamicsLowSpeed
-from .components.compute_polar import ComputePolar, PolarType
-from .components.compute_reynolds import ComputeReynolds
-from .components.initialize_cl import InitializeClPolar
-from .components.oswald import OswaldCoefficient
+from fastoad.module_management.service_registry import RegisterOpenMDAOSystem, RegisterSubmodel
+from .constants import (
+    PolarType,
+    SERVICE_CD0,
+    SERVICE_CD_TRIM,
+    SERVICE_INDUCED_DRAG_COEFFICIENT,
+    SERVICE_INITIALIZE_CL,
+    SERVICE_OSWALD_COEFFICIENT,
+    SERVICE_POLAR,
+    SERVICE_REYNOLDS_COEFFICIENT,
+    SERVICE_LOW_SPEED_CL_AOA,
+)
 
 
 @RegisterOpenMDAOSystem("fastoad.aerodynamics.lowspeed.legacy", domain=ModelDomain.AERODYNAMICS)
@@ -40,21 +36,49 @@ class AerodynamicsLowSpeed(om.Group):
     """
 
     def setup(self):
-        self.add_subsystem("compute_low_speed_aero", ComputeAerodynamicsLowSpeed(), promotes=["*"])
+        low_speed_option = {"low_speed_aero": True}
+
+        self.add_subsystem(
+            "compute_low_speed_aero",
+            RegisterSubmodel.get_submodel(SERVICE_LOW_SPEED_CL_AOA),
+            promotes=["*"],
+        )
         ivc = om.IndepVarComp("data:aerodynamics:aircraft:takeoff:mach", val=0.2)
         self.add_subsystem("mach_low_speed", ivc, promotes=["*"])
+
+        self.add_subsystem(
+            "compute_oswald_coeff",
+            RegisterSubmodel.get_submodel(SERVICE_OSWALD_COEFFICIENT, low_speed_option),
+            promotes=["*"],
+        )
+        self.add_subsystem(
+            "compute_induced_drag_coeff",
+            RegisterSubmodel.get_submodel(SERVICE_INDUCED_DRAG_COEFFICIENT, low_speed_option),
+            promotes=["*"],
+        )
+        self.add_subsystem(
+            "comp_re",
+            RegisterSubmodel.get_submodel(SERVICE_REYNOLDS_COEFFICIENT, low_speed_option),
+            promotes=["*"],
+        )
+        self.add_subsystem(
+            "initialize_cl",
+            RegisterSubmodel.get_submodel(SERVICE_INITIALIZE_CL, low_speed_option),
+            promotes=["*"],
+        )
+        self.add_subsystem(
+            "cd0_wing",
+            RegisterSubmodel.get_submodel(SERVICE_CD0, low_speed_option),
+            promotes=["*"],
+        )
         self.add_subsystem(
-            "compute_oswald_coeff", OswaldCoefficient(low_speed_aero=True), promotes=["*"]
+            "cd_trim",
+            RegisterSubmodel.get_submodel(SERVICE_CD_TRIM, low_speed_option),
+            promotes=["*"],
         )
-        self.add_subsystem("comp_re", ComputeReynolds(low_speed_aero=True), promotes=["*"])
-        self.add_subsystem("initialize_cl", InitializeClPolar(low_speed_aero=True), promotes=["*"])
-        self.add_subsystem("cd0_wing", Cd0Wing(low_speed_aero=True), promotes=["*"])
-        self.add_subsystem("cd0_fuselage", Cd0Fuselage(low_speed_aero=True), promotes=["*"])
-        self.add_subsystem("cd0_ht", Cd0HorizontalTail(low_speed_aero=True), promotes=["*"])
-        self.add_subsystem("cd0_vt", Cd0VerticalTail(low_speed_aero=True), promotes=["*"])
+        polar_type_option = {"polar_type": PolarType.LOW_SPEED}
         self.add_subsystem(
-            "cd0_nac_pylons", Cd0NacelleAndPylons(low_speed_aero=True), promotes=["*"]
+            "get_polar",
+            RegisterSubmodel.get_submodel(SERVICE_POLAR, polar_type_option),
+            promotes=["*"],
         )
-        self.add_subsystem("cd0_total", Cd0Total(low_speed_aero=True), promotes=["*"])
-        self.add_subsystem("cd_trim", CdTrim(low_speed_aero=True), promotes=["*"])
-        self.add_subsystem("get_polar", ComputePolar(type=PolarType.LOW_SPEED), promotes=["*"])
diff --git a/src/fastoad/models/aerodynamics/aerodynamics_takeoff.py b/src/fastoad/models/aerodynamics/aerodynamics_takeoff.py
index e3a805e46..59841eb9d 100644
--- a/src/fastoad/models/aerodynamics/aerodynamics_takeoff.py
+++ b/src/fastoad/models/aerodynamics/aerodynamics_takeoff.py
@@ -1,3 +1,4 @@
+"""Computation of aerodynamic characteristics at takeoff."""
 #  This file is part of FAST-OAD : A framework for rapid Overall Aircraft Design
 #  Copyright (C) 2021 ONERA & ISAE-SUPAERO
 #  FAST is free software: you can redistribute it and/or modify
@@ -14,9 +15,8 @@
 import openmdao.api as om
 
 from fastoad.module_management.constants import ModelDomain
-from fastoad.module_management.service_registry import RegisterOpenMDAOSystem
-from .components.compute_polar import ComputePolar, PolarType
-from .components.high_lift_aero import ComputeDeltaHighLift
+from fastoad.module_management.service_registry import RegisterOpenMDAOSystem, RegisterSubmodel
+from .constants import PolarType, SERVICE_POLAR, SERVICE_HIGH_LIFT
 
 
 @RegisterOpenMDAOSystem("fastoad.aerodynamics.takeoff.legacy", domain=ModelDomain.AERODYNAMICS)
@@ -28,5 +28,16 @@ class AerodynamicsTakeoff(om.Group):
     """
 
     def setup(self):
-        self.add_subsystem("delta_cl_cd", ComputeDeltaHighLift(landing_flag=False), promotes=["*"])
-        self.add_subsystem("polar", ComputePolar(type=PolarType.TAKEOFF), promotes=["*"])
+        landing_flag_option = {"landing_flag": False}
+        self.add_subsystem(
+            "delta_cl_cd",
+            RegisterSubmodel.get_submodel(SERVICE_HIGH_LIFT, landing_flag_option),
+            promotes=["*"],
+        )
+
+        polar_type_option = {"polar_type": PolarType.TAKEOFF}
+        self.add_subsystem(
+            "polar",
+            RegisterSubmodel.get_submodel(SERVICE_POLAR, polar_type_option),
+            promotes=["*"],
+        )
diff --git a/src/fastoad/models/aerodynamics/components/cd0.py b/src/fastoad/models/aerodynamics/components/cd0.py
index 9e20582e4..ad2e8dd26 100644
--- a/src/fastoad/models/aerodynamics/components/cd0.py
+++ b/src/fastoad/models/aerodynamics/components/cd0.py
@@ -1,5 +1,6 @@
+"""Computation of form drag for whole aircraft."""
 #  This file is part of FAST-OAD : A framework for rapid Overall Aircraft Design
-#  Copyright (C) 2020  ONERA & ISAE-SUPAERO
+#  Copyright (C) 2021 ONERA & ISAE-SUPAERO
 #  FAST is free software: you can redistribute it and/or modify
 #  it under the terms of the GNU General Public License as published by
 #  the Free Software Foundation, either version 3 of the License, or
@@ -13,40 +14,59 @@
 
 import openmdao.api as om
 
-from .cd0_fuselage import Cd0Fuselage
-from .cd0_ht import Cd0HorizontalTail
-from .cd0_nacelle_pylons import Cd0NacelleAndPylons
-from .cd0_total import Cd0Total
-from .cd0_vt import Cd0VerticalTail
-from .cd0_wing import Cd0Wing
+from fastoad.module_management.service_registry import RegisterSubmodel
+from ..constants import (
+    SERVICE_CD0,
+    SERVICE_CD0_FUSELAGE,
+    SERVICE_CD0_HORIZONTAL_TAIL,
+    SERVICE_CD0_NACELLES_PYLONS,
+    SERVICE_CD0_SUM,
+    SERVICE_CD0_VERTICAL_TAIL,
+    SERVICE_CD0_WING,
+)
 
 
+@RegisterSubmodel(SERVICE_CD0, "fastoad.submodel.aerodynamics.CD0.legacy")
 class CD0(om.Group):
+    """
+    Computation of form drag for whole aircraft.
+
+    Computes and sums the drag coefficients of all components.
+    Interaction drag is assumed to be taken into account at component level.
+    """
+
     def initialize(self):
         self.options.declare("low_speed_aero", default=False, types=bool)
 
     def setup(self):
+        low_speed_option = {"low_speed_aero": self.options["low_speed_aero"]}
         self.add_subsystem(
-            "cd0_wing", Cd0Wing(low_speed_aero=self.options["low_speed_aero"]), promotes=["*"]
+            "cd0_wing",
+            RegisterSubmodel.get_submodel(SERVICE_CD0_WING, low_speed_option),
+            promotes=["*"],
         )
         self.add_subsystem(
             "cd0_fuselage",
-            Cd0Fuselage(low_speed_aero=self.options["low_speed_aero"]),
+            RegisterSubmodel.get_submodel(SERVICE_CD0_FUSELAGE, low_speed_option),
             promotes=["*"],
         )
         self.add_subsystem(
             "cd0_ht",
-            Cd0HorizontalTail(low_speed_aero=self.options["low_speed_aero"]),
+            RegisterSubmodel.get_submodel(SERVICE_CD0_HORIZONTAL_TAIL, low_speed_option),
             promotes=["*"],
         )
         self.add_subsystem(
-            "cd0_vt", Cd0VerticalTail(low_speed_aero=self.options["low_speed_aero"]), promotes=["*"]
+            "cd0_vt",
+            RegisterSubmodel.get_submodel(SERVICE_CD0_VERTICAL_TAIL, low_speed_option),
+            promotes=["*"],
         )
         self.add_subsystem(
             "cd0_nac_pylons",
-            Cd0NacelleAndPylons(low_speed_aero=self.options["low_speed_aero"]),
+            RegisterSubmodel.get_submodel(SERVICE_CD0_NACELLES_PYLONS, low_speed_option),
             promotes=["*"],
         )
         self.add_subsystem(
-            "cd0_total", Cd0Total(low_speed_aero=self.options["low_speed_aero"]), promotes=["*"]
+            "cd0_total",
+            RegisterSubmodel.get_submodel(SERVICE_CD0_SUM, low_speed_option),
+            promotes=["*"],
         )
diff --git a/src/fastoad/models/aerodynamics/components/cd0_fuselage.py b/src/fastoad/models/aerodynamics/components/cd0_fuselage.py
index 6399d9800..17d219703 100644
--- a/src/fastoad/models/aerodynamics/components/cd0_fuselage.py
+++ b/src/fastoad/models/aerodynamics/components/cd0_fuselage.py
@@ -1,6 +1,4 @@
-"""
-    FAST - Copyright (c) 2016 ONERA ISAE
-"""
+"""Computation of form drag for fuselage."""
 #  This file is part of FAST-OAD : A framework for rapid Overall Aircraft Design
 #  Copyright (C) 2021 ONERA & ISAE-SUPAERO
 #  FAST is free software: you can redistribute it and/or modify
@@ -14,19 +12,23 @@
 #  You should have received a copy of the GNU General Public License
 #  along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
-import math
-
 import numpy as np
-from openmdao.core.explicitcomponent import ExplicitComponent
+import openmdao.api as om
+
+from fastoad.module_management.service_registry import RegisterSubmodel
+from .utils.friction_drag import get_flat_plate_friction_drag_coefficient
+from ..constants import SERVICE_CD0_FUSELAGE
 
 
-class Cd0Fuselage(ExplicitComponent):
+@RegisterSubmodel(SERVICE_CD0_FUSELAGE, "fastoad.submodel.aerodynamics.CD0.fuselage.legacy")
+class Cd0Fuselage(om.ExplicitComponent):
+    """Computation of form drag for fuselage."""
+
     def initialize(self):
         self.options.declare("low_speed_aero", default=False, types=bool)
 
     def setup(self):
-        self.low_speed_aero = self.options["low_speed_aero"]
-        if self.low_speed_aero:
+        if self.options["low_speed_aero"]:
             self.add_input("data:aerodynamics:wing:low_speed:reynolds", val=np.nan)
             self.add_input(
                 "data:aerodynamics:aircraft:low_speed:CL", shape_by_conn=True, val=np.nan
@@ -54,13 +56,13 @@ def setup(self):
     def setup_partials(self):
         self.declare_partials("*", "*", method="fd")
 
-    def compute(self, inputs, outputs):
+    def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
         height_max = inputs["data:geometry:fuselage:maximum_height"]
         width_max = inputs["data:geometry:fuselage:maximum_width"]
         wet_area_fus = inputs["data:geometry:fuselage:wetted_area"]
         wing_area = inputs["data:geometry:wing:area"]
         fus_length = inputs["data:geometry:fuselage:length"]
-        if self.low_speed_aero:
+        if self.options["low_speed_aero"]:
             cl = inputs["data:aerodynamics:aircraft:low_speed:CL"]
             mach = inputs["data:aerodynamics:aircraft:takeoff:mach"]
             reynolds = inputs["data:aerodynamics:wing:low_speed:reynolds"]
@@ -69,12 +71,10 @@ def compute(self, inputs, outputs):
             mach = inputs["data:TLAR:cruise_mach"]
             reynolds = inputs["data:aerodynamics:wing:cruise:reynolds"]
 
-        cf_fus = 0.455 / (
-            (1 + 0.144 * mach ** 2) ** 0.65 * (math.log10(reynolds * fus_length)) ** 2.58
-        )
+        cf_fus = get_flat_plate_friction_drag_coefficient(fus_length, mach, reynolds)
 
         cd0_friction_fus = (
-            (0.98 + 0.745 * math.sqrt(height_max * width_max) / fus_length)
+            (0.98 + 0.745 * np.sqrt(height_max * width_max) / fus_length)
             * cf_fus
             * wet_area_fus
             / wing_area
@@ -86,7 +86,7 @@ def compute(self, inputs, outputs):
         )
         cd0_fus = cd0_friction_fus + cd0_upsweep_fus
 
-        if self.low_speed_aero:
+        if self.options["low_speed_aero"]:
             outputs["data:aerodynamics:fuselage:low_speed:CD0"] = cd0_fus
         else:
             outputs["data:aerodynamics:fuselage:cruise:CD0"] = cd0_fus
diff --git a/src/fastoad/models/aerodynamics/components/cd0_ht.py b/src/fastoad/models/aerodynamics/components/cd0_ht.py
index 86eea6ba3..3b96b72d7 100644
--- a/src/fastoad/models/aerodynamics/components/cd0_ht.py
+++ b/src/fastoad/models/aerodynamics/components/cd0_ht.py
@@ -1,7 +1,4 @@
-"""
-    FAST - Copyright (c) 2016 ONERA ISAE
-"""
-
+"""Computation of form drag for Horizontal Tail Plane."""
 #  This file is part of FAST-OAD : A framework for rapid Overall Aircraft Design
 #  Copyright (C) 2021 ONERA & ISAE-SUPAERO
 #  FAST is free software: you can redistribute it and/or modify
@@ -15,25 +12,35 @@
 #  You should have received a copy of the GNU General Public License
 #  along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
-import math
-
 import numpy as np
-from openmdao.core.explicitcomponent import ExplicitComponent
+import openmdao.api as om
+
+from fastoad.module_management.service_registry import RegisterSubmodel
+from .utils.cd0_lifting_surface import LiftingSurfaceGeometry, compute_cd0_lifting_surface
+from ..constants import SERVICE_CD0_HORIZONTAL_TAIL
+
+
+@RegisterSubmodel(
+    SERVICE_CD0_HORIZONTAL_TAIL, "fastoad.submodel.aerodynamics.CD0.horizontal_tail.legacy"
+)
+class Cd0HorizontalTail(om.ExplicitComponent):
+    """
+    Computation of form drag for Horizontal Tail Plane.
 
+    See :meth:`~fastoad.models.aerodynamics.components.utils.cd0_lifting_surface` for used method.
+    """
 
-class Cd0HorizontalTail(ExplicitComponent):
     def initialize(self):
         self.options.declare("low_speed_aero", default=False, types=bool)
 
     def setup(self):
-        self.low_speed_aero = self.options["low_speed_aero"]
-
         self.add_input("data:geometry:horizontal_tail:MAC:length", val=np.nan, units="m")
         self.add_input("data:geometry:horizontal_tail:thickness_ratio", val=np.nan)
         self.add_input("data:geometry:horizontal_tail:sweep_25", val=np.nan, units="deg")
         self.add_input("data:geometry:horizontal_tail:wetted_area", val=np.nan, units="m**2")
         self.add_input("data:geometry:wing:area", val=np.nan, units="m**2")
-        if self.low_speed_aero:
+        self.add_input("settings:aerodynamics:wing:CD:fuselage_interaction", val=0.04)
+        if self.options["low_speed_aero"]:
             self.add_input("data:aerodynamics:wing:low_speed:reynolds", val=np.nan)
             self.add_input("data:aerodynamics:aircraft:takeoff:mach", val=np.nan)
             self.add_output("data:aerodynamics:horizontal_tail:low_speed:CD0")
@@ -45,29 +52,26 @@ def setup(self):
     def setup_partials(self):
         self.declare_partials("*", "*", method="fd")
 
-    def compute(self, inputs, outputs):
-        el_ht = inputs["data:geometry:horizontal_tail:thickness_ratio"]
-        ht_length = inputs["data:geometry:horizontal_tail:MAC:length"]
-        sweep_25_ht = inputs["data:geometry:horizontal_tail:sweep_25"]
-        wet_area_ht = inputs["data:geometry:horizontal_tail:wetted_area"]
+    def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
+        ht_geometry = LiftingSurfaceGeometry(
+            thickness_ratio=inputs["data:geometry:horizontal_tail:thickness_ratio"],
+            MAC_length=inputs["data:geometry:horizontal_tail:MAC:length"],
+            sweep_angle_25=inputs["data:geometry:horizontal_tail:sweep_25"],
+            wet_area=inputs["data:geometry:horizontal_tail:wetted_area"],
+            cambered=False,
+            interaction_coeff=0.01,
+        )
         wing_area = inputs["data:geometry:wing:area"]
-        if self.low_speed_aero:
+        if self.options["low_speed_aero"]:
             mach = inputs["data:aerodynamics:aircraft:takeoff:mach"]
             reynolds = inputs["data:aerodynamics:wing:low_speed:reynolds"]
         else:
             mach = inputs["data:TLAR:cruise_mach"]
             reynolds = inputs["data:aerodynamics:wing:cruise:reynolds"]
 
-        ki_arrow_cd0 = 0.04
-
-        cf_ht = 0.455 / (
-            (1 + 0.144 * mach ** 2) ** 0.65 * (math.log10(reynolds * ht_length)) ** 2.58
-        )
-        ke_cd0_ht = 4.688 * el_ht ** 2 + 3.146 * el_ht
-        k_phi_cd0_ht = 1 - 0.000178 * sweep_25_ht ** 2 - 0.0065 * sweep_25_ht
-        cd0_ht = (ke_cd0_ht * k_phi_cd0_ht + ki_arrow_cd0 / 4 + 1) * cf_ht * wet_area_ht / wing_area
+        cd0_ht = compute_cd0_lifting_surface(ht_geometry, mach, reynolds, wing_area)
 
-        if self.low_speed_aero:
+        if self.options["low_speed_aero"]:
             outputs["data:aerodynamics:horizontal_tail:low_speed:CD0"] = cd0_ht
         else:
             outputs["data:aerodynamics:horizontal_tail:cruise:CD0"] = cd0_ht
diff --git a/src/fastoad/models/aerodynamics/components/cd0_nacelle_pylons.py b/src/fastoad/models/aerodynamics/components/cd0_nacelles_pylons.py
similarity index 73%
rename from src/fastoad/models/aerodynamics/components/cd0_nacelle_pylons.py
rename to src/fastoad/models/aerodynamics/components/cd0_nacelles_pylons.py
index b053bcd27..62ea895e3 100644
--- a/src/fastoad/models/aerodynamics/components/cd0_nacelle_pylons.py
+++ b/src/fastoad/models/aerodynamics/components/cd0_nacelles_pylons.py
@@ -1,7 +1,4 @@
-"""
-    FAST - Copyright (c) 2016 ONERA ISAE
-"""
-
+"""Computation of form drag for nacelles and pylons."""
 #  This file is part of FAST-OAD : A framework for rapid Overall Aircraft Design
 #  Copyright (C) 2021 ONERA & ISAE-SUPAERO
 #  FAST is free software: you can redistribute it and/or modify
@@ -15,20 +12,25 @@
 #  You should have received a copy of the GNU General Public License
 #  along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
-import math
-
 import numpy as np
-from openmdao.core.explicitcomponent import ExplicitComponent
+import openmdao.api as om
+
+from fastoad.module_management.service_registry import RegisterSubmodel
+from .utils.friction_drag import get_flat_plate_friction_drag_coefficient
+from ..constants import SERVICE_CD0_NACELLES_PYLONS
+
 
+@RegisterSubmodel(
+    SERVICE_CD0_NACELLES_PYLONS, "fastoad.submodel.aerodynamics.CD0.nacelles_pylons.legacy"
+)
+class Cd0NacellesAndPylons(om.ExplicitComponent):
+    """Computation of form drag for nacelles and pylons."""
 
-class Cd0NacelleAndPylons(ExplicitComponent):
     def initialize(self):
         self.options.declare("low_speed_aero", default=False, types=bool)
 
     def setup(self):
-        self.low_speed_aero = self.options["low_speed_aero"]
-
-        if self.low_speed_aero:
+        if self.options["low_speed_aero"]:
             self.add_input("data:aerodynamics:wing:low_speed:reynolds", val=np.nan)
             self.add_input("data:aerodynamics:aircraft:takeoff:mach", val=np.nan)
             self.add_output("data:aerodynamics:nacelles:low_speed:CD0")
@@ -50,42 +52,46 @@ def setup(self):
     def setup_partials(self):
         self.declare_partials("*", "*", method="fd")
 
-    def compute(self, inputs, outputs):
-        pylon_length = inputs["data:geometry:propulsion:pylon:length"]
-        nac_length = inputs["data:geometry:propulsion:nacelle:length"]
-        wet_area_pylon = inputs["data:geometry:propulsion:pylon:wetted_area"]
-        wet_area_nac = inputs["data:geometry:propulsion:nacelle:wetted_area"]
+    def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
         n_engines = inputs["data:geometry:propulsion:engine:count"]
-        fan_length = inputs["data:geometry:propulsion:fan:length"]
         wing_area = inputs["data:geometry:wing:area"]
-        if self.low_speed_aero:
+        if self.options["low_speed_aero"]:
             mach = inputs["data:aerodynamics:aircraft:takeoff:mach"]
             reynolds = inputs["data:aerodynamics:wing:low_speed:reynolds"]
         else:
             mach = inputs["data:TLAR:cruise_mach"]
             reynolds = inputs["data:aerodynamics:wing:cruise:reynolds"]
 
-        cf_pylon = 0.455 / (
-            (1 + 0.144 * mach ** 2) ** 0.65 * (math.log10(reynolds * pylon_length)) ** 2.58
-        )
-        cf_nac = 0.455 / (
-            (1 + 0.144 * mach ** 2) ** 0.65 * (math.log10(reynolds * nac_length)) ** 2.58
-        )
+        cd0_pylon = self._compute_cd0_for_pylons(inputs, n_engines, wing_area, mach, reynolds)
+        cd0_nac = self._compute_cd0_for_nacelles(inputs, n_engines, wing_area, mach, reynolds)
 
-        # cd0 Pylon
+        if self.options["low_speed_aero"]:
+            outputs["data:aerodynamics:pylons:low_speed:CD0"] = cd0_pylon
+            outputs["data:aerodynamics:nacelles:low_speed:CD0"] = cd0_nac
+        else:
+            outputs["data:aerodynamics:pylons:cruise:CD0"] = cd0_pylon
+            outputs["data:aerodynamics:nacelles:cruise:CD0"] = cd0_nac
+
+    @staticmethod
+    def _compute_cd0_for_pylons(inputs, n_engines, wing_area, mach, reynolds):
+        pylon_length = inputs["data:geometry:propulsion:pylon:length"]
+        wet_area_pylon = inputs["data:geometry:propulsion:pylon:wetted_area"]
+
+        cf_pylon = get_flat_plate_friction_drag_coefficient(pylon_length, mach, reynolds)
         el_pylon = 0.06
         ke_cd0_pylon = 4.688 * el_pylon ** 2 + 3.146 * el_pylon
         cd0_pylon = n_engines * (1 + ke_cd0_pylon) * cf_pylon * wet_area_pylon / wing_area
+        return cd0_pylon
+
+    @staticmethod
+    def _compute_cd0_for_nacelles(inputs, n_engines, wing_area, mach, reynolds):
+        nac_length = inputs["data:geometry:propulsion:nacelle:length"]
+        wet_area_nac = inputs["data:geometry:propulsion:nacelle:wetted_area"]
+        fan_length = inputs["data:geometry:propulsion:fan:length"]
 
-        # cd0 Nacelles
+        cf_nac = get_flat_plate_friction_drag_coefficient(nac_length, mach, reynolds)
         e_fan = 0.22
         kn_cd0_nac = 1 + 0.05 + 5.8 * e_fan / fan_length
         cd0_int_nac = 0.0002
         cd0_nac = n_engines * (kn_cd0_nac * cf_nac * wet_area_nac / wing_area + cd0_int_nac)
-
-        if self.low_speed_aero:
-            outputs["data:aerodynamics:pylons:low_speed:CD0"] = cd0_pylon
-            outputs["data:aerodynamics:nacelles:low_speed:CD0"] = cd0_nac
-        else:
-            outputs["data:aerodynamics:pylons:cruise:CD0"] = cd0_pylon
-            outputs["data:aerodynamics:nacelles:cruise:CD0"] = cd0_nac
+        return cd0_nac
diff --git a/src/fastoad/models/aerodynamics/components/cd0_total.py b/src/fastoad/models/aerodynamics/components/cd0_total.py
index 19aa9fc8c..0ed2c6629 100644
--- a/src/fastoad/models/aerodynamics/components/cd0_total.py
+++ b/src/fastoad/models/aerodynamics/components/cd0_total.py
@@ -1,6 +1,4 @@
-"""
-    FAST - Copyright (c) 2016 ONERA ISAE
-"""
+"""Sum of form drags from aircraft components."""
 #  This file is part of FAST-OAD : A framework for rapid Overall Aircraft Design
 #  Copyright (C) 2021 ONERA & ISAE-SUPAERO
 #  FAST is free software: you can redistribute it and/or modify
@@ -15,19 +13,23 @@
 #  along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
 import numpy as np
-from openmdao.core.explicitcomponent import ExplicitComponent
+import openmdao.api as om
 
+from fastoad.module_management.service_registry import RegisterSubmodel
+from ..constants import SERVICE_CD0_SUM
+
+
+@RegisterSubmodel(SERVICE_CD0_SUM, "fastoad.submodel.aerodynamics.CD0.sum.legacy")
+class Cd0Total(om.ExplicitComponent):
+    """Computes the sum of form drags from aircraft components."""
 
-class Cd0Total(ExplicitComponent):
     def initialize(self):
         self.options.declare("low_speed_aero", default=False, types=bool)
 
     def setup(self):
-        self.low_speed_aero = self.options["low_speed_aero"]
-
         self.add_input("data:geometry:aircraft:wetted_area", val=np.nan, units="m**2")
 
-        if self.low_speed_aero:
+        if self.options["low_speed_aero"]:
             self.add_input("data:aerodynamics:wing:low_speed:CD0", shape_by_conn=True, val=np.nan)
             self.add_input(
                 "data:aerodynamics:fuselage:low_speed:CD0", shape_by_conn=True, val=np.nan
@@ -40,6 +42,14 @@ def setup(self):
                 "data:aerodynamics:aircraft:low_speed:CD0",
                 copy_shape="data:aerodynamics:wing:low_speed:CD0",
             )
+            self.add_output(
+                "data:aerodynamics:aircraft:low_speed:CD0:clean",
+                copy_shape="data:aerodynamics:wing:low_speed:CD0",
+            )
+            self.add_output(
+                "data:aerodynamics:aircraft:low_speed:CD0:parasitic",
+                copy_shape="data:aerodynamics:wing:low_speed:CD0",
+            )
         else:
             self.add_input("data:aerodynamics:wing:cruise:CD0", shape_by_conn=True, val=np.nan)
             self.add_input("data:aerodynamics:fuselage:cruise:CD0", shape_by_conn=True, val=np.nan)
@@ -51,13 +61,21 @@ def setup(self):
                 "data:aerodynamics:aircraft:cruise:CD0",
                 copy_shape="data:aerodynamics:wing:cruise:CD0",
             )
+            self.add_output(
+                "data:aerodynamics:aircraft:cruise:CD0:clean",
+                copy_shape="data:aerodynamics:wing:cruise:CD0",
+            )
+            self.add_output(
+                "data:aerodynamics:aircraft:cruise:CD0:parasitic",
+                copy_shape="data:aerodynamics:wing:cruise:CD0",
+            )
 
     def setup_partials(self):
         self.declare_partials("*", "*", method="fd")
 
-    def compute(self, inputs, outputs):
+    def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
         wet_area_total = inputs["data:geometry:aircraft:wetted_area"]
-        if self.low_speed_aero:
+        if self.options["low_speed_aero"]:
             cd0_wing = inputs["data:aerodynamics:wing:low_speed:CD0"]
             cd0_fus = inputs["data:aerodynamics:fuselage:low_speed:CD0"]
             cd0_ht = inputs["data:aerodynamics:horizontal_tail:low_speed:CD0"]
@@ -72,7 +90,6 @@ def compute(self, inputs, outputs):
             cd0_nac = inputs["data:aerodynamics:nacelles:cruise:CD0"]
             cd0_pylon = inputs["data:aerodynamics:pylons:cruise:CD0"]
 
-        k_techno = 1
         k_parasite = (
             -2.39 * pow(10, -12) * wet_area_total ** 3
             + 2.58 * pow(10, -8) * wet_area_total ** 2
@@ -80,10 +97,16 @@ def compute(self, inputs, outputs):
             + 0.163
         )
 
-        cd0_total_hs = cd0_wing + cd0_fus + cd0_ht + cd0_vt + cd0_nac + cd0_pylon
-        cd0_total = cd0_total_hs * (1.0 + k_parasite * k_techno)
+        cd0_total_clean = cd0_wing + cd0_fus + cd0_ht + cd0_vt + cd0_nac + cd0_pylon
+        cd0_total = cd0_total_clean * (1.0 + k_parasite)
 
-        if self.low_speed_aero:
+        if self.options["low_speed_aero"]:
             outputs["data:aerodynamics:aircraft:low_speed:CD0"] = cd0_total
+            outputs["data:aerodynamics:aircraft:low_speed:CD0:clean"] = cd0_total_clean
+            outputs["data:aerodynamics:aircraft:low_speed:CD0:parasitic"] = (
+                cd0_total - cd0_total_clean
+            )
         else:
             outputs["data:aerodynamics:aircraft:cruise:CD0"] = cd0_total
+            outputs["data:aerodynamics:aircraft:cruise:CD0:clean"] = cd0_total_clean
+            outputs["data:aerodynamics:aircraft:cruise:CD0:parasitic"] = cd0_total - cd0_total_clean
diff --git a/src/fastoad/models/aerodynamics/components/cd0_vt.py b/src/fastoad/models/aerodynamics/components/cd0_vt.py
index c1a26979c..ef6715074 100644
--- a/src/fastoad/models/aerodynamics/components/cd0_vt.py
+++ b/src/fastoad/models/aerodynamics/components/cd0_vt.py
@@ -1,7 +1,4 @@
-"""
-    FAST - Copyrigvt (c) 2016 ONERA ISAE
-"""
-
+"""Computation of form drag for Vertical Tail Plane."""
 #  This file is part of FAST-OAD : A framework for rapid Overall Aircraft Design
 #  Copyright (C) 2021 ONERA & ISAE-SUPAERO
 #  FAST is free software: you can redistribute it and/or modify
@@ -15,25 +12,34 @@
 #  You should have received a copy of the GNU General Public License
 #  along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
-import math
-
 import numpy as np
-from openmdao.core.explicitcomponent import ExplicitComponent
+import openmdao.api as om
+
+from fastoad.module_management.service_registry import RegisterSubmodel
+from .utils.cd0_lifting_surface import LiftingSurfaceGeometry, compute_cd0_lifting_surface
+from ..constants import SERVICE_CD0_VERTICAL_TAIL
+
+
+@RegisterSubmodel(
+    SERVICE_CD0_VERTICAL_TAIL, "fastoad.submodel.aerodynamics.CD0.vertical_tail.legacy"
+)
+class Cd0VerticalTail(om.ExplicitComponent):
+    """
+    Computation of form drag for Vertical Tail Plane.
 
+    See :meth:`~fastoad.models.aerodynamics.components.utils.cd0_lifting_surface` for used method.
+    """
 
-class Cd0VerticalTail(ExplicitComponent):
     def initialize(self):
         self.options.declare("low_speed_aero", default=False, types=bool)
 
     def setup(self):
-        self.low_speed_aero = self.options["low_speed_aero"]
-
         self.add_input("data:geometry:vertical_tail:MAC:length", val=np.nan, units="m")
         self.add_input("data:geometry:vertical_tail:thickness_ratio", val=np.nan)
         self.add_input("data:geometry:vertical_tail:sweep_25", val=np.nan, units="deg")
         self.add_input("data:geometry:vertical_tail:wetted_area", val=np.nan, units="m**2")
         self.add_input("data:geometry:wing:area", val=np.nan, units="m**2")
-        if self.low_speed_aero:
+        if self.options["low_speed_aero"]:
             self.add_input("data:aerodynamics:wing:low_speed:reynolds", val=np.nan)
             self.add_input("data:aerodynamics:aircraft:takeoff:mach", val=np.nan)
             self.add_output("data:aerodynamics:vertical_tail:low_speed:CD0")
@@ -45,29 +51,26 @@ def setup(self):
     def setup_partials(self):
         self.declare_partials("*", "*", method="fd")
 
-    def compute(self, inputs, outputs):
-        el_vt = inputs["data:geometry:vertical_tail:thickness_ratio"]
-        vt_length = inputs["data:geometry:vertical_tail:MAC:length"]
-        sweep_25_vt = inputs["data:geometry:vertical_tail:sweep_25"]
-        wet_area_vt = inputs["data:geometry:vertical_tail:wetted_area"]
+    def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
+        vt_geometry = LiftingSurfaceGeometry(
+            thickness_ratio=inputs["data:geometry:vertical_tail:thickness_ratio"],
+            MAC_length=inputs["data:geometry:vertical_tail:MAC:length"],
+            sweep_angle_25=inputs["data:geometry:vertical_tail:sweep_25"],
+            wet_area=inputs["data:geometry:vertical_tail:wetted_area"],
+            cambered=False,
+            interaction_coeff=0.005,
+        )
         wing_area = inputs["data:geometry:wing:area"]
-        if self.low_speed_aero:
+        if self.options["low_speed_aero"]:
             mach = inputs["data:aerodynamics:aircraft:takeoff:mach"]
             reynolds = inputs["data:aerodynamics:wing:low_speed:reynolds"]
         else:
             mach = inputs["data:TLAR:cruise_mach"]
             reynolds = inputs["data:aerodynamics:wing:cruise:reynolds"]
 
-        ki_arrow_cd0 = 0.04
-
-        cf_vt = 0.455 / (
-            (1 + 0.144 * mach ** 2) ** 0.65 * (math.log10(reynolds * vt_length)) ** 2.58
-        )
-        ke_cd0_vt = 4.688 * el_vt ** 2 + 3.146 * el_vt
-        k_phi_cd0_vt = 1 - 0.000178 * sweep_25_vt ** 2 - 0.0065 * sweep_25_vt
-        cd0_vt = (ke_cd0_vt * k_phi_cd0_vt + ki_arrow_cd0 / 8 + 1) * cf_vt * wet_area_vt / wing_area
+        cd0_vt = compute_cd0_lifting_surface(vt_geometry, mach, reynolds, wing_area)
 
-        if self.low_speed_aero:
+        if self.options["low_speed_aero"]:
             outputs["data:aerodynamics:vertical_tail:low_speed:CD0"] = cd0_vt
         else:
             outputs["data:aerodynamics:vertical_tail:cruise:CD0"] = cd0_vt
diff --git a/src/fastoad/models/aerodynamics/components/cd0_wing.py b/src/fastoad/models/aerodynamics/components/cd0_wing.py
index 778039c6f..1955e20f1 100644
--- a/src/fastoad/models/aerodynamics/components/cd0_wing.py
+++ b/src/fastoad/models/aerodynamics/components/cd0_wing.py
@@ -1,6 +1,4 @@
-"""
-    FAST - Copyright (c) 2016 ONERA ISAE
-"""
+"""Computation of form drag for wing."""
 #  This file is part of FAST-OAD : A framework for rapid Overall Aircraft Design
 #  Copyright (C) 2021 ONERA & ISAE-SUPAERO
 #  FAST is free software: you can redistribute it and/or modify
@@ -15,17 +13,26 @@
 #  along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
 import numpy as np
-from openmdao.core.explicitcomponent import ExplicitComponent
+import openmdao.api as om
 
+from fastoad.module_management.service_registry import RegisterSubmodel
+from .utils.cd0_lifting_surface import LiftingSurfaceGeometry, compute_cd0_lifting_surface
+from ..constants import SERVICE_CD0_WING
+
+
+@RegisterSubmodel(SERVICE_CD0_WING, "fastoad.submodel.aerodynamics.CD0.wing.legacy")
+class Cd0Wing(om.ExplicitComponent):
+    """
+    Computation of form drag for wing.
+
+    See :meth:`~fastoad.models.aerodynamics.components.utils.cd0_lifting_surface` for used method.
+    """
 
-class Cd0Wing(ExplicitComponent):
     def initialize(self):
         self.options.declare("low_speed_aero", default=False, types=bool)
 
     def setup(self):
-        self.low_speed_aero = self.options["low_speed_aero"]
-
-        if self.low_speed_aero:
+        if self.options["low_speed_aero"]:
             self.add_input("data:aerodynamics:wing:low_speed:reynolds", val=np.nan)
             self.add_input(
                 "data:aerodynamics:aircraft:low_speed:CL", shape_by_conn=True, val=np.nan
@@ -53,13 +60,17 @@ def setup(self):
     def setup_partials(self):
         self.declare_partials("*", "*", method="fd")
 
-    def compute(self, inputs, outputs):
+    def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
+        wing_geometry = LiftingSurfaceGeometry(
+            thickness_ratio=inputs["data:geometry:wing:thickness_ratio"],
+            MAC_length=inputs["data:geometry:wing:MAC:length"],
+            sweep_angle_25=inputs["data:geometry:wing:sweep_25"],
+            wet_area=inputs["data:geometry:wing:wetted_area"],
+            cambered=True,
+            interaction_coeff=0.04,
+        )
         wing_area = inputs["data:geometry:wing:area"]
-        wet_area_wing = inputs["data:geometry:wing:wetted_area"]
-        el_aero = inputs["data:geometry:wing:thickness_ratio"]
-        sweep_25 = inputs["data:geometry:wing:sweep_25"]
-        l0_wing = inputs["data:geometry:wing:MAC:length"]
-        if self.low_speed_aero:
+        if self.options["low_speed_aero"]:
             cl = inputs["data:aerodynamics:aircraft:low_speed:CL"]
             mach = inputs["data:aerodynamics:aircraft:takeoff:mach"]
             reynolds = inputs["data:aerodynamics:wing:low_speed:reynolds"]
@@ -68,30 +79,9 @@ def compute(self, inputs, outputs):
             mach = inputs["data:TLAR:cruise_mach"]
             reynolds = inputs["data:aerodynamics:wing:cruise:reynolds"]
 
-        ki_arrow_cd0 = 0.04
-        # Friction coefficients
-        cf_wing = 0.455 / ((1 + 0.144 * mach ** 2) ** 0.65 * (np.log10(reynolds * l0_wing)) ** 2.58)
-
-        # cd0 wing
-        # factor of relative thickness
-        ke_cd0_wing = 4.688 * el_aero ** 2 + 3.146 * el_aero
-        k_phi_cd0_wing = 1 - 0.000178 * sweep_25 ** 2 - 0.0065 * sweep_25
-
-        kc_cd0_wing = (
-            2.859 * (cl / np.cos(np.radians(sweep_25)) ** 2) ** 3
-            - 1.849 * (cl / np.cos(np.radians(sweep_25)) ** 2) ** 2
-            + 0.382 * (cl / np.cos(np.radians(sweep_25)) ** 2)
-            + 0.06
-        )  # sweep factor
-
-        cd0_wing = (
-            ((ke_cd0_wing + kc_cd0_wing) * k_phi_cd0_wing + ki_arrow_cd0 + 1)
-            * cf_wing
-            * wet_area_wing
-            / wing_area
-        )
+        cd0_wing = compute_cd0_lifting_surface(wing_geometry, mach, reynolds, wing_area, cl)
 
-        if self.low_speed_aero:
+        if self.options["low_speed_aero"]:
             outputs["data:aerodynamics:wing:low_speed:CD0"] = cd0_wing
         else:
             outputs["data:aerodynamics:wing:cruise:CD0"] = cd0_wing
diff --git a/src/fastoad/models/aerodynamics/components/cd_compressibility.py b/src/fastoad/models/aerodynamics/components/cd_compressibility.py
index ffc4a9ab1..b85f7c3d0 100644
--- a/src/fastoad/models/aerodynamics/components/cd_compressibility.py
+++ b/src/fastoad/models/aerodynamics/components/cd_compressibility.py
@@ -1,7 +1,4 @@
-"""
-Compressibility drag computation.
-"""
-
+"""Compressibility drag computation."""
 #  This file is part of FAST-OAD : A framework for rapid Overall Aircraft Design
 #  Copyright (C) 2021 ONERA & ISAE-SUPAERO
 #  FAST is free software: you can redistribute it and/or modify
@@ -16,15 +13,22 @@
 #  along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
 import numpy as np
-from openmdao.core.explicitcomponent import ExplicitComponent
+import openmdao.api as om
+
+from fastoad.module_management.service_registry import RegisterSubmodel
+from ..constants import SERVICE_CD_COMPRESSIBILITY
 
 
-class CdCompressibility(ExplicitComponent):
+@RegisterSubmodel(
+    SERVICE_CD_COMPRESSIBILITY, "fastoad.submodel.aerodynamics.CD.compressibility.legacy"
+)
+class CdCompressibility(om.ExplicitComponent):
     """
     Computation of drag increment due to compressibility effects.
 
     Formula from §4.2.4 of :cite:`supaero:2014`. This formula can be used for aircraft
     before year 2000.
+
     Earlier aircraft have more optimized wing profiles that are expected to limit the
     compressibility drag below 2 drag counts. Until a better model can be provided, the
     variable `tuning:aerodynamics:aircraft:cruise:CD:compressibility:characteristic_mach_increment`
@@ -51,7 +55,7 @@ def setup(self):
     def setup_partials(self):
         self.declare_partials("*", "*", method="fd")
 
-    def compute(self, inputs, outputs):
+    def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
         cl = inputs["data:aerodynamics:aircraft:cruise:CL"]
         m = inputs["data:TLAR:cruise_mach"]
         max_cd_comp = inputs["tuning:aerodynamics:aircraft:cruise:CD:compressibility:max_value"]
diff --git a/src/fastoad/models/aerodynamics/components/cd_trim.py b/src/fastoad/models/aerodynamics/components/cd_trim.py
index 165c391af..93624a084 100644
--- a/src/fastoad/models/aerodynamics/components/cd_trim.py
+++ b/src/fastoad/models/aerodynamics/components/cd_trim.py
@@ -1,7 +1,4 @@
-"""
-    FAST - Copyright (c) 2016 ONERA ISAE
-"""
-
+"""Computation of trim drag."""
 #  This file is part of FAST-OAD : A framework for rapid Overall Aircraft Design
 #  Copyright (C) 2021 ONERA & ISAE-SUPAERO
 #  FAST is free software: you can redistribute it and/or modify
@@ -16,17 +13,21 @@
 #  along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
 import numpy as np
-from openmdao.core.explicitcomponent import ExplicitComponent
+import openmdao.api as om
+
+from fastoad.module_management.service_registry import RegisterSubmodel
+from ..constants import SERVICE_CD_TRIM
 
 
-class CdTrim(ExplicitComponent):
+@RegisterSubmodel(SERVICE_CD_TRIM, "fastoad.submodel.aerodynamics.CD.trim.legacy")
+class CdTrim(om.ExplicitComponent):
+    """Computation of trim drag."""
+
     def initialize(self):
         self.options.declare("low_speed_aero", default=False, types=bool)
 
     def setup(self):
-        self.low_speed_aero = self.options["low_speed_aero"]
-
-        if self.low_speed_aero:
+        if self.options["low_speed_aero"]:
             self.add_input(
                 "data:aerodynamics:aircraft:low_speed:CL", shape_by_conn=True, val=np.nan
             )
@@ -44,8 +45,8 @@ def setup(self):
     def setup_partials(self):
         self.declare_partials("*", "*", method="fd")
 
-    def compute(self, inputs, outputs):
-        if self.low_speed_aero:
+    def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
+        if self.options["low_speed_aero"]:
             cl = inputs["data:aerodynamics:aircraft:low_speed:CL"]
         else:
             cl = inputs["data:aerodynamics:aircraft:cruise:CL"]
@@ -55,7 +56,7 @@ def compute(self, inputs, outputs):
         for cl_val in cl:
             cd_trim.append(5.89 * pow(10, -4) * cl_val)
 
-        if self.low_speed_aero:
+        if self.options["low_speed_aero"]:
             outputs["data:aerodynamics:aircraft:low_speed:CD:trim"] = cd_trim
         else:
             outputs["data:aerodynamics:aircraft:cruise:CD:trim"] = cd_trim
diff --git a/src/fastoad/models/aerodynamics/components/compute_low_speed_aero.py b/src/fastoad/models/aerodynamics/components/compute_low_speed_aero.py
index fe22a507c..b310cc21d 100644
--- a/src/fastoad/models/aerodynamics/components/compute_low_speed_aero.py
+++ b/src/fastoad/models/aerodynamics/components/compute_low_speed_aero.py
@@ -1,6 +1,4 @@
-"""
-    FAST - Copyright (c) 2016 ONERA ISAE
-"""
+"""Computation of CL characteristics at low speed."""
 #  This file is part of FAST-OAD : A framework for rapid Overall Aircraft Design
 #  Copyright (C) 2021 ONERA & ISAE-SUPAERO
 #  FAST is free software: you can redistribute it and/or modify
@@ -14,18 +12,23 @@
 #  You should have received a copy of the GNU General Public License
 #  along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
-from math import pi, sqrt, tan
-
 import numpy as np
-from openmdao.core.explicitcomponent import ExplicitComponent
+import openmdao.api as om
+
+from fastoad.module_management.service_registry import RegisterSubmodel
+from ..constants import SERVICE_LOW_SPEED_CL_AOA
 
 
-class ComputeAerodynamicsLowSpeed(ExplicitComponent):
+@RegisterSubmodel(SERVICE_LOW_SPEED_CL_AOA, "fastoad.submodel.aerodynamics.low_speed.AoA.legacy")
+class ComputeAerodynamicsLowSpeed(om.ExplicitComponent):
     """
-    Czalpha from Raymer Eq 12.6
-    TODO: complete source
+    Computes CL gradient and CL at low speed.
+
+    CL gradient from :cite:`raymer:1999` Eq 12.6
     """
 
+    # TODO: complete source
+
     def setup(self):
         self.add_input("data:geometry:fuselage:maximum_width", val=np.nan, units="m")
         self.add_input("data:geometry:fuselage:maximum_height", val=np.nan, units="m")
@@ -43,7 +46,7 @@ def setup(self):
     def setup_partials(self):
         self.declare_partials("*", "*", method="fd")
 
-    def compute(self, inputs, outputs):
+    def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
         width_max = inputs["data:geometry:fuselage:maximum_width"]
         height_max = inputs["data:geometry:fuselage:maximum_height"]
         span = inputs["data:geometry:wing:span"]
@@ -56,27 +59,27 @@ def compute(self, inputs, outputs):
 
         mach = 0.2
 
-        beta = sqrt(1 - mach ** 2)
-        d_f = sqrt(width_max * height_max)
-        fact_F = 1.07 * (1 + d_f / span) ** 2
+        beta = np.sqrt(1 - mach ** 2)
+        d_f = np.sqrt(width_max * height_max)
+        fuselage_lift_factor = 1.07 * (1 + d_f / span) ** 2
         lambda_wing_eff = lambda_wing * (1 + 1.9 * l4_wing * el_ext / span)
         cl_alpha_wing_low = (
             2
-            * pi
+            * np.pi
             * lambda_wing_eff
             / (
                 2
-                + sqrt(
+                + np.sqrt(
                     4
                     + lambda_wing_eff ** 2
                     * beta ** 2
                     / 0.95 ** 2
-                    * (1 + (tan(sweep_25 / 180.0 * pi)) ** 2 / beta ** 2)
+                    * (1 + (np.tan(sweep_25 / 180.0 * np.pi)) ** 2 / beta ** 2)
                 )
             )
             * (wing_area - l2_wing * width_max)
             / wing_area
-            * fact_F
+            * fuselage_lift_factor
         )
 
         outputs["data:aerodynamics:aircraft:takeoff:CL_alpha"] = cl_alpha_wing_low
diff --git a/src/fastoad/models/aerodynamics/components/compute_max_cl_landing.py b/src/fastoad/models/aerodynamics/components/compute_max_cl_landing.py
index 4443e8ce9..11775c6fd 100644
--- a/src/fastoad/models/aerodynamics/components/compute_max_cl_landing.py
+++ b/src/fastoad/models/aerodynamics/components/compute_max_cl_landing.py
@@ -1,7 +1,4 @@
-"""
-    FAST - Copyright (c) 2016 ONERA ISAE
-"""
-
+"""Computation of max CL in landing conditions."""
 #  This file is part of FAST-OAD : A framework for rapid Overall Aircraft Design
 #  Copyright (C) 2021 ONERA & ISAE-SUPAERO
 #  FAST is free software: you can redistribute it and/or modify
@@ -18,8 +15,14 @@
 import numpy as np
 from openmdao.core.explicitcomponent import ExplicitComponent
 
+from fastoad.module_management.service_registry import RegisterSubmodel
+from ..constants import SERVICE_LANDING_MAX_CL
+
 
+@RegisterSubmodel(SERVICE_LANDING_MAX_CL, "fastoad.submodel.aerodynamics.landing.max_CL.legacy")
 class ComputeMaxClLanding(ExplicitComponent):
+    """Computation of max CL in landing conditions."""
+
     def setup(self):
         self.add_input("data:aerodynamics:aircraft:landing:CL_max_clean", val=np.nan)
         self.add_input("data:aerodynamics:high_lift_devices:landing:CL", val=np.nan)
diff --git a/src/fastoad/models/aerodynamics/components/compute_polar.py b/src/fastoad/models/aerodynamics/components/compute_polar.py
index 007e7952f..60919c7fa 100644
--- a/src/fastoad/models/aerodynamics/components/compute_polar.py
+++ b/src/fastoad/models/aerodynamics/components/compute_polar.py
@@ -1,6 +1,4 @@
-"""
-    FAST - Copyright (c) 2016 ONERA ISAE
-"""
+"""Computation of CL and CD for whole aircraft."""
 #  This file is part of FAST-OAD : A framework for rapid Overall Aircraft Design
 #  Copyright (C) 2021 ONERA & ISAE-SUPAERO
 #  FAST is free software: you can redistribute it and/or modify
@@ -15,20 +13,18 @@
 #  along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
 import numpy as np
-from aenum import Enum
-from openmdao.core.explicitcomponent import ExplicitComponent
+import openmdao.api as om
 
+from fastoad.module_management.service_registry import RegisterSubmodel
+from ..constants import PolarType, SERVICE_POLAR
 
-class PolarType(Enum):
-    HIGH_SPEED = "high_speed"
-    LOW_SPEED = "low_speed"
-    TAKEOFF = "takeoff"
-    LANDING = "landing"
 
+@RegisterSubmodel(SERVICE_POLAR, "fastoad.submodel.aerodynamics.polar.legacy")
+class ComputePolar(om.ExplicitComponent):
+    """Computation of CL and CD for whole aircraft."""
 
-class ComputePolar(ExplicitComponent):
     def initialize(self):
-        self.options.declare("type", default=PolarType.HIGH_SPEED, types=PolarType)
+        self.options.declare("polar_type", default=PolarType.HIGH_SPEED, types=PolarType)
 
     def setup(self):
         self.add_input("tuning:aerodynamics:aircraft:cruise:CD:k", val=np.nan)
@@ -36,7 +32,7 @@ def setup(self):
         self.add_input("tuning:aerodynamics:aircraft:cruise:CD:winglet_effect:k", val=np.nan)
         self.add_input("tuning:aerodynamics:aircraft:cruise:CD:winglet_effect:offset", val=np.nan)
 
-        if self.options["type"] != PolarType.HIGH_SPEED:
+        if self.options["polar_type"] != PolarType.HIGH_SPEED:
             self.add_input(
                 "data:aerodynamics:aircraft:low_speed:CL", shape_by_conn=True, val=np.nan
             )
@@ -50,7 +46,7 @@ def setup(self):
                 "data:aerodynamics:aircraft:low_speed:induced_drag_coefficient", val=np.nan
             )
 
-            if self.options["type"] == PolarType.TAKEOFF:
+            if self.options["polar_type"] == PolarType.TAKEOFF:
                 self.add_input("data:aerodynamics:high_lift_devices:takeoff:CL", val=np.nan)
                 self.add_input("data:aerodynamics:high_lift_devices:takeoff:CD", val=np.nan)
                 self.add_output(
@@ -62,7 +58,7 @@ def setup(self):
                     copy_shape="data:aerodynamics:aircraft:low_speed:CL",
                 )
 
-            elif self.options["type"] == PolarType.LANDING:
+            elif self.options["polar_type"] == PolarType.LANDING:
                 self.add_input("data:aerodynamics:high_lift_devices:landing:CL", val=np.nan)
                 self.add_input("data:aerodynamics:high_lift_devices:landing:CD", val=np.nan)
                 self.add_output(
@@ -77,13 +73,15 @@ def setup(self):
                     "data:aerodynamics:aircraft:landing:CD",
                     copy_shape="data:aerodynamics:aircraft:low_speed:CL",
                 )
-            else:
+            elif self.options["polar_type"] == PolarType.LOW_SPEED:
                 self.add_output(
                     "data:aerodynamics:aircraft:low_speed:CD",
                     copy_shape="data:aerodynamics:aircraft:low_speed:CL",
                 )
+            else:
+                raise AttributeError(f'Unknown polar type: {self.options["polar_type"]}')
 
-        else:
+        elif self.options["polar_type"] == PolarType.HIGH_SPEED:
             self.add_input("data:aerodynamics:aircraft:cruise:CL", shape_by_conn=True, val=np.nan)
             self.add_input("data:aerodynamics:aircraft:cruise:CD0", shape_by_conn=True, val=np.nan)
             self.add_input(
@@ -104,31 +102,36 @@ def setup(self):
             self.add_output("data:aerodynamics:aircraft:cruise:optimal_CL")
             self.add_output("data:aerodynamics:aircraft:cruise:optimal_CD")
 
+        else:
+            raise AttributeError(f'Unknown polar type: {self.options["polar_type"]}')
+
     def setup_partials(self):
         self.declare_partials("*", "*", method="fd")
 
-    def compute(self, inputs, outputs):
+    def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
         k_cd = inputs["tuning:aerodynamics:aircraft:cruise:CD:k"]
         offset_cd = inputs["tuning:aerodynamics:aircraft:cruise:CD:offset"]
         k_winglet_cd = inputs["tuning:aerodynamics:aircraft:cruise:CD:winglet_effect:k"]
         offset_winglet_cd = inputs["tuning:aerodynamics:aircraft:cruise:CD:winglet_effect:offset"]
-        if self.options["type"] != PolarType.HIGH_SPEED:
+        if self.options["polar_type"] != PolarType.HIGH_SPEED:
             cl = inputs["data:aerodynamics:aircraft:low_speed:CL"]
             cd0 = inputs["data:aerodynamics:aircraft:low_speed:CD0"]
             cd_trim = inputs["data:aerodynamics:aircraft:low_speed:CD:trim"]
             cd_c = 0.0
             coef_k = inputs["data:aerodynamics:aircraft:low_speed:induced_drag_coefficient"]
-            if self.options["type"] == PolarType.TAKEOFF:
+            if self.options["polar_type"] == PolarType.TAKEOFF:
                 delta_cl_hl = inputs["data:aerodynamics:high_lift_devices:takeoff:CL"]
                 delta_cd_hl = inputs["data:aerodynamics:high_lift_devices:takeoff:CD"]
-            elif self.options["type"] == PolarType.LANDING:
+            elif self.options["polar_type"] == PolarType.LANDING:
                 delta_cl_hl = inputs["data:aerodynamics:high_lift_devices:landing:CL"]
                 delta_cd_hl = inputs["data:aerodynamics:high_lift_devices:landing:CD"]
-            else:
+            elif self.options["polar_type"] == PolarType.LOW_SPEED:
                 delta_cl_hl = 0.0
                 delta_cd_hl = 0.0
+            else:
+                raise AttributeError(f'Unknown polar type: {self.options["polar_type"]}')
 
-        else:
+        elif self.options["polar_type"] == PolarType.HIGH_SPEED:
             cl = inputs["data:aerodynamics:aircraft:cruise:CL"]
             cd0 = inputs["data:aerodynamics:aircraft:cruise:CD0"]
             cd_trim = inputs["data:aerodynamics:aircraft:cruise:CD:trim"]
@@ -136,18 +139,20 @@ def compute(self, inputs, outputs):
             coef_k = inputs["data:aerodynamics:aircraft:cruise:induced_drag_coefficient"]
             delta_cl_hl = 0.0
             delta_cd_hl = 0.0
+        else:
+            raise AttributeError(f'Unknown polar type: {self.options["polar_type"]}')
 
         cl = cl + delta_cl_hl
         cd = (
             cd0 + cd_c + cd_trim + coef_k * cl ** 2 * k_winglet_cd + offset_winglet_cd + delta_cd_hl
         ) * k_cd + offset_cd
 
-        if self.options["type"] == PolarType.LOW_SPEED:
+        if self.options["polar_type"] == PolarType.LOW_SPEED:
             outputs["data:aerodynamics:aircraft:low_speed:CD"] = cd
-        elif self.options["type"] == PolarType.TAKEOFF:
+        elif self.options["polar_type"] == PolarType.TAKEOFF:
             outputs["data:aerodynamics:aircraft:takeoff:CL"] = cl
             outputs["data:aerodynamics:aircraft:takeoff:CD"] = cd
-        elif self.options["type"] == PolarType.LANDING:
+        elif self.options["polar_type"] == PolarType.LANDING:
             outputs["data:aerodynamics:aircraft:landing:CL"] = cl
             outputs["data:aerodynamics:aircraft:landing:CD"] = cd
         else:
diff --git a/src/fastoad/models/aerodynamics/components/compute_reynolds.py b/src/fastoad/models/aerodynamics/components/compute_reynolds.py
index 37e988d98..87b0f02ed 100644
--- a/src/fastoad/models/aerodynamics/components/compute_reynolds.py
+++ b/src/fastoad/models/aerodynamics/components/compute_reynolds.py
@@ -1,7 +1,4 @@
-"""
-    FAST - Copyright (c) 2016 ONERA ISAE
-"""
-
+"""Computation of Reynolds number"""
 #  This file is part of FAST-OAD : A framework for rapid Overall Aircraft Design
 #  Copyright (C) 2021 ONERA & ISAE-SUPAERO
 #  FAST is free software: you can redistribute it and/or modify
@@ -20,16 +17,21 @@
 from openmdao.core.explicitcomponent import ExplicitComponent
 
 from fastoad.model_base import AtmosphereSI
+from fastoad.module_management.service_registry import RegisterSubmodel
+from ..constants import SERVICE_REYNOLDS_COEFFICIENT
 
 
+@RegisterSubmodel(
+    SERVICE_REYNOLDS_COEFFICIENT, "fastoad.submodel.aerodynamics.reynolds_coefficient.legacy"
+)
 class ComputeReynolds(ExplicitComponent):
+    """Computation of Reynolds number"""
+
     def initialize(self):
         self.options.declare("low_speed_aero", default=False, types=bool)
 
     def setup(self):
-        self.low_speed_aero = self.options["low_speed_aero"]
-
-        if self.low_speed_aero:
+        if self.options["low_speed_aero"]:
             self.add_input("data:aerodynamics:aircraft:takeoff:mach", val=np.nan)
             self.add_output("data:aerodynamics:wing:low_speed:reynolds")
         else:
@@ -40,8 +42,8 @@ def setup(self):
     def setup_partials(self):
         self.declare_partials("*", "*", method="fd")
 
-    def compute(self, inputs, outputs):
-        if self.low_speed_aero:
+    def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
+        if self.options["low_speed_aero"]:
             mach = inputs["data:aerodynamics:aircraft:takeoff:mach"]
             altitude = 0.0
         else:
@@ -52,7 +54,7 @@ def compute(self, inputs, outputs):
         atm.mach = mach
         reynolds = atm.unitary_reynolds
 
-        if self.low_speed_aero:
+        if self.options["low_speed_aero"]:
             outputs["data:aerodynamics:wing:low_speed:reynolds"] = reynolds
         else:
             outputs["data:aerodynamics:wing:cruise:reynolds"] = reynolds
diff --git a/src/fastoad/models/aerodynamics/components/high_lift_aero.py b/src/fastoad/models/aerodynamics/components/high_lift_aero.py
index 589100019..f2383494a 100644
--- a/src/fastoad/models/aerodynamics/components/high_lift_aero.py
+++ b/src/fastoad/models/aerodynamics/components/high_lift_aero.py
@@ -20,11 +20,14 @@
 import openmdao.api as om
 from scipy import interpolate
 
+from fastoad.module_management.service_registry import RegisterSubmodel
 from . import resources
+from ..constants import SERVICE_HIGH_LIFT
 
 LIFT_EFFECTIVENESS_FILENAME = "interpolation of lift effectiveness.txt"
 
 
+@RegisterSubmodel(SERVICE_HIGH_LIFT, "fastoad.submodel.aerodynamics.high_lift.legacy")
 class ComputeDeltaHighLift(om.ExplicitComponent):
     """
     Provides lift and drag increments due to high-lift devices
diff --git a/src/fastoad/models/aerodynamics/components/initialize_cl.py b/src/fastoad/models/aerodynamics/components/initialize_cl.py
index b846fdd4a..a122b3f13 100644
--- a/src/fastoad/models/aerodynamics/components/initialize_cl.py
+++ b/src/fastoad/models/aerodynamics/components/initialize_cl.py
@@ -1,6 +1,4 @@
-"""
-    FAST - Copyright (c) 2016 ONERA ISAE
-"""
+"""Initialization of CL vector."""
 #  This file is part of FAST-OAD : A framework for rapid Overall Aircraft Design
 #  Copyright (C) 2021 ONERA & ISAE-SUPAERO
 #  FAST is free software: you can redistribute it and/or modify
@@ -15,24 +13,26 @@
 #  along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
 import numpy as np
-from openmdao.core.explicitcomponent import ExplicitComponent
+import openmdao.api as om
 
-from fastoad.models.aerodynamics.constants import POLAR_POINT_COUNT
+from fastoad.module_management.service_registry import RegisterSubmodel
+from ..constants import POLAR_POINT_COUNT, SERVICE_INITIALIZE_CL
 
 
-class InitializeClPolar(ExplicitComponent):
+@RegisterSubmodel(SERVICE_INITIALIZE_CL, "fastoad.submodel.aerodynamics.initialize_CL.legacy")
+class InitializeClPolar(om.ExplicitComponent):
+    """Initialization of CL vector."""
+
     def initialize(self):
         self.options.declare("low_speed_aero", default=False, types=bool)
 
     def setup(self):
-        self.low_speed_aero = self.options["low_speed_aero"]
-
         self.add_input("tuning:aerodynamics:aircraft:cruise:CL:k", val=np.nan)
         self.add_input("tuning:aerodynamics:aircraft:cruise:CL:offset", val=np.nan)
         self.add_input("tuning:aerodynamics:aircraft:cruise:CL:winglet_effect:k", val=np.nan)
         self.add_input("tuning:aerodynamics:aircraft:cruise:CL:winglet_effect:offset", val=np.nan)
 
-        if self.low_speed_aero:
+        if self.options["low_speed_aero"]:
             self.add_output("data:aerodynamics:aircraft:low_speed:CL", shape=POLAR_POINT_COUNT)
         else:
             self.add_output("data:aerodynamics:aircraft:cruise:CL", shape=POLAR_POINT_COUNT)
@@ -40,7 +40,7 @@ def setup(self):
     def setup_partials(self):
         self.declare_partials("*", "*", method="fd")
 
-    def compute(self, inputs, outputs):
+    def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
         k_cl = inputs["tuning:aerodynamics:aircraft:cruise:CL:k"]
         offset_cl = inputs["tuning:aerodynamics:aircraft:cruise:CL:offset"]
         k_winglet_cl = inputs["tuning:aerodynamics:aircraft:cruise:CL:winglet_effect:k"]
@@ -49,7 +49,7 @@ def compute(self, inputs, outputs):
         # FIXME: initialization of CL range should be done more directly, without these coefficients
         cl = np.arange(0.0, 1.5, 0.01) * k_cl * k_winglet_cl + offset_cl + offset_winglet_cl
 
-        if self.low_speed_aero:
+        if self.options["low_speed_aero"]:
             outputs["data:aerodynamics:aircraft:low_speed:CL"] = cl
         else:
             outputs["data:aerodynamics:aircraft:cruise:CL"] = cl
diff --git a/src/fastoad/models/aerodynamics/components/oswald.py b/src/fastoad/models/aerodynamics/components/oswald.py
index d507c37c6..ccafa8a6d 100644
--- a/src/fastoad/models/aerodynamics/components/oswald.py
+++ b/src/fastoad/models/aerodynamics/components/oswald.py
@@ -1,7 +1,6 @@
 """
 Computation of Oswald coefficient
 """
-
 #  This file is part of FAST-OAD : A framework for rapid Overall Aircraft Design
 #  Copyright (C) 2021 ONERA & ISAE-SUPAERO
 #  FAST is free software: you can redistribute it and/or modify
@@ -15,17 +14,63 @@
 #  You should have received a copy of the GNU General Public License
 #  along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
-import math
-
 import numpy as np
-from openmdao.core.explicitcomponent import ExplicitComponent
+import openmdao.api as om
 
+from fastoad.module_management.service_registry import RegisterSubmodel
+from ..constants import SERVICE_INDUCED_DRAG_COEFFICIENT, SERVICE_OSWALD_COEFFICIENT
 
-class OswaldCoefficient(ExplicitComponent):
-    # TODO: Document equations. Cite sources (M. Nita and D. Scholz)
-    # FIXME: output the real Oswald coefficient (coef_e instead of coef_k)
+
+@RegisterSubmodel(
+    SERVICE_INDUCED_DRAG_COEFFICIENT,
+    "fastoad.submodel.aerodynamics.induced_drag_coefficient.legacy",
+)
+class InducedDragCoefficient(om.ExplicitComponent):
+    """Computes the coefficient that should be multiplied by CL**2 to get induced drag."""
+
+    def initialize(self):
+        self.options.declare("low_speed_aero", default=False, types=bool)
+
+    def setup(self):
+        self.add_input("data:geometry:wing:area", val=np.nan, units="m**2")
+        self.add_input("data:geometry:wing:span", val=np.nan, units="m")
+
+        if self.options["low_speed_aero"]:
+            self.add_input("data:aerodynamics:aircraft:low_speed:oswald_coefficient", val=np.nan)
+            self.add_output("data:aerodynamics:aircraft:low_speed:induced_drag_coefficient")
+        else:
+            self.add_input("data:aerodynamics:aircraft:cruise:oswald_coefficient", val=np.nan)
+            self.add_output("data:aerodynamics:aircraft:cruise:induced_drag_coefficient")
+
+    def setup_partials(self):
+        self.declare_partials("*", "*", method="fd")
+
+    def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
+        wing_area = inputs["data:geometry:wing:area"]
+        span = inputs["data:geometry:wing:span"] / np.cos(5.0 / 180 * np.pi)
+        aspect_ratio = span ** 2 / wing_area
+
+        if self.options["low_speed_aero"]:
+            coef_e = inputs["data:aerodynamics:aircraft:low_speed:oswald_coefficient"]
+        else:
+            coef_e = inputs["data:aerodynamics:aircraft:cruise:oswald_coefficient"]
+
+        coef_k = 1.0 / (np.pi * aspect_ratio * coef_e)
+
+        if self.options["low_speed_aero"]:
+            outputs["data:aerodynamics:aircraft:low_speed:induced_drag_coefficient"] = coef_k
+        else:
+            outputs["data:aerodynamics:aircraft:cruise:induced_drag_coefficient"] = coef_k
+
+
+@RegisterSubmodel(
+    SERVICE_OSWALD_COEFFICIENT, "fastoad.submodel.aerodynamics.oswald_coefficient.legacy"
+)
+class OswaldCoefficient(om.ExplicitComponent):
     """ Computes Oswald efficiency number """
 
+    # TODO: Document equations. Cite sources (M. Nita and D. Scholz)
+
     def initialize(self):
         self.options.declare("low_speed_aero", default=False, types=bool)
 
@@ -40,17 +85,17 @@ def setup(self):
 
         if self.options["low_speed_aero"]:
             self.add_input("data:aerodynamics:aircraft:takeoff:mach", val=np.nan)
-            self.add_output("data:aerodynamics:aircraft:low_speed:induced_drag_coefficient")
+            self.add_output("data:aerodynamics:aircraft:low_speed:oswald_coefficient")
         else:
             self.add_input("data:TLAR:cruise_mach", val=np.nan)
-            self.add_output("data:aerodynamics:aircraft:cruise:induced_drag_coefficient")
+            self.add_output("data:aerodynamics:aircraft:cruise:oswald_coefficient")
 
     def setup_partials(self):
         self.declare_partials("*", "*", method="fd")
 
-    def compute(self, inputs, outputs):
+    def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
         wing_area = inputs["data:geometry:wing:area"]
-        span = inputs["data:geometry:wing:span"] / math.cos(5.0 / 180 * math.pi)
+        span = inputs["data:geometry:wing:span"] / np.cos(5.0 / 180 * np.pi)
         height_fus = inputs["data:geometry:fuselage:maximum_height"]
         width_fus = inputs["data:geometry:fuselage:maximum_width"]
         l2_wing = inputs["data:geometry:wing:root:chord"]
@@ -62,9 +107,9 @@ def compute(self, inputs, outputs):
             mach = inputs["data:TLAR:cruise_mach"]
 
         aspect_ratio = span ** 2 / wing_area
-        df = math.sqrt(width_fus * height_fus)
+        df = np.sqrt(width_fus * height_fus)
         lamda = l4_wing / l2_wing
-        delta_lamda = -0.357 + 0.45 * math.exp(0.0375 * sweep_25 / 180.0 * math.pi)
+        delta_lamda = -0.357 + 0.45 * np.exp(0.0375 * sweep_25 / 180.0 * np.pi)
         lamda = lamda - delta_lamda
         f_lamda = (
             0.0524 * lamda ** 4 - 0.15 * lamda ** 3 + 0.1659 * lamda ** 2 - 0.0706 * lamda + 0.0119
@@ -78,9 +123,8 @@ def compute(self, inputs, outputs):
 
         ke_f = 1 - 2 * (df / span) ** 2
         coef_e = e_theory * ke_f * ke_m * 0.9
-        coef_k = 1.0 / (math.pi * aspect_ratio * coef_e)
 
         if self.options["low_speed_aero"]:
-            outputs["data:aerodynamics:aircraft:low_speed:induced_drag_coefficient"] = coef_k
+            outputs["data:aerodynamics:aircraft:low_speed:oswald_coefficient"] = coef_e
         else:
-            outputs["data:aerodynamics:aircraft:cruise:induced_drag_coefficient"] = coef_k
+            outputs["data:aerodynamics:aircraft:cruise:oswald_coefficient"] = coef_e
diff --git a/src/fastoad/models/aerodynamics/components/utils/__init__.py b/src/fastoad/models/aerodynamics/components/utils/__init__.py
new file mode 100644
index 000000000..385ae0611
--- /dev/null
+++ b/src/fastoad/models/aerodynamics/components/utils/__init__.py
@@ -0,0 +1,12 @@
+#  This file is part of FAST-OAD : A framework for rapid Overall Aircraft Design
+#  Copyright (C) 2021 ONERA & ISAE-SUPAERO
+#  FAST is free software: you can redistribute it and/or modify
+#  it under the terms of the GNU General Public License as published by
+#  the Free Software Foundation, either version 3 of the License, or
+#  (at your option) any later version.
+#  This program is distributed in the hope that it will be useful,
+#  but WITHOUT ANY WARRANTY; without even the implied warranty of
+#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#  GNU General Public License for more details.
+#  You should have received a copy of the GNU General Public License
+#  along with this program.  If not, see <https://www.gnu.org/licenses/>.
diff --git a/src/fastoad/models/aerodynamics/components/utils/cd0_lifting_surface.py b/src/fastoad/models/aerodynamics/components/utils/cd0_lifting_surface.py
new file mode 100644
index 000000000..1c09d89f1
--- /dev/null
+++ b/src/fastoad/models/aerodynamics/components/utils/cd0_lifting_surface.py
@@ -0,0 +1,91 @@
+"""Computation of CD0 for a lifting surface."""
+#  This file is part of FAST-OAD : A framework for rapid Overall Aircraft Design
+#  Copyright (C) 2021 ONERA & ISAE-SUPAERO
+#  FAST is free software: you can redistribute it and/or modify
+#  it under the terms of the GNU General Public License as published by
+#  the Free Software Foundation, either version 3 of the License, or
+#  (at your option) any later version.
+#  This program is distributed in the hope that it will be useful,
+#  but WITHOUT ANY WARRANTY; without even the implied warranty of
+#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#  GNU General Public License for more details.
+#  You should have received a copy of the GNU General Public License
+#  along with this program.  If not, see <https://www.gnu.org/licenses/>.
+
+from dataclasses import dataclass
+
+import numpy as np
+
+from fastoad.models.aerodynamics.components.utils.friction_drag import (
+    get_flat_plate_friction_drag_coefficient,
+)
+
+
+@dataclass
+class LiftingSurfaceGeometry:
+    """Minimum geometry data for computation of CD0 of lifting surfaces."""
+
+    thickness_ratio: float  #: average thickness ratio
+    MAC_length: float  #: length of Mean Aerodynamic Chord
+    sweep_angle_25: float  #: sweep angle at 25% chord, in degrees
+    cambered: bool  #: True if airfoil is cambered
+    wet_area: float  #: wet surface area of the lifting surface
+    interaction_coeff: float  #: ratio of additional drag due to interaction effects
+
+
+def compute_cd0_lifting_surface(
+    geometry: LiftingSurfaceGeometry,
+    mach: float,
+    reynolds: float,
+    wing_area: float,
+    lift_coefficient: float = 0.0,
+):
+    """
+    Computes CD0 for a lifting surface.
+
+    Friction coefficient is assessed from :cite:`raymer:1999` (Eq 12.27).
+    Corrections for lifting surfaces are from :cite:`supaero:2014`.
+
+    :param geometry: definition of lifting surface geometry
+    :param mach: Mach number
+    :param reynolds: Reynolds number
+    :param wing_area: wing area (will be used for getting CD specific to wing area
+    :param lift_coefficient: needed if wing is cambered
+    :return: CD0 value
+    """
+    # Drag coefficient for flat plate (
+    cf = get_flat_plate_friction_drag_coefficient(geometry.MAC_length, mach, reynolds)
+
+    # Contribution of relative thickness
+    thickness_contribution = (
+        4.688 * geometry.thickness_ratio ** 2 + 3.146 * geometry.thickness_ratio
+    )
+
+    # Contribution of camber
+    sweep_25 = np.radians(geometry.sweep_angle_25)
+    if geometry.cambered:
+        camber_contribution = (
+            2.859 * (lift_coefficient / np.cos(sweep_25) ** 2) ** 3
+            - 1.849 * (lift_coefficient / np.cos(sweep_25) ** 2) ** 2
+            + 0.382 * (lift_coefficient / np.cos(sweep_25) ** 2)
+            + 0.06
+        )
+    else:
+        camber_contribution = 0.0
+
+    # Correction from sweep angle
+    sweep_correction = (
+        1 - 0.000178 * geometry.sweep_angle_25 ** 2 - 0.0065 * geometry.sweep_angle_25
+    )
+
+    cd0_ht = (
+        (
+            (thickness_contribution + camber_contribution) * sweep_correction
+            + geometry.interaction_coeff
+            + 1
+        )
+        * cf
+        * geometry.wet_area
+        / wing_area
+    )
+    return cd0_ht
diff --git a/src/fastoad/models/aerodynamics/components/utils/friction_drag.py b/src/fastoad/models/aerodynamics/components/utils/friction_drag.py
new file mode 100644
index 000000000..0b04beded
--- /dev/null
+++ b/src/fastoad/models/aerodynamics/components/utils/friction_drag.py
@@ -0,0 +1,27 @@
+"""Computation of friction drag."""
+#  This file is part of FAST-OAD : A framework for rapid Overall Aircraft Design
+#  Copyright (C) 2021 ONERA & ISAE-SUPAERO
+#  FAST is free software: you can redistribute it and/or modify
+#  it under the terms of the GNU General Public License as published by
+#  the Free Software Foundation, either version 3 of the License, or
+#  (at your option) any later version.
+#  This program is distributed in the hope that it will be useful,
+#  but WITHOUT ANY WARRANTY; without even the implied warranty of
+#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#  GNU General Public License for more details.
+#  You should have received a copy of the GNU General Public License
+#  along with this program.  If not, see <https://www.gnu.org/licenses/>.
+
+import numpy as np
+
+
+def get_flat_plate_friction_drag_coefficient(length, mach, reynolds):
+    """
+
+    :param length: flat plate length in meters
+    :param mach: Mach number
+    :param reynolds: Reynolds number
+    :return: Drag coefficient w.r.t. a surface of area length*1 m**2
+    """
+    c_f = 0.455 / ((1 + 0.144 * mach ** 2) ** 0.65 * (np.log10(reynolds * length)) ** 2.58)
+    return c_f
diff --git a/src/fastoad/models/aerodynamics/constants.py b/src/fastoad/models/aerodynamics/constants.py
index 55c44f4fc..fce802572 100644
--- a/src/fastoad/models/aerodynamics/constants.py
+++ b/src/fastoad/models/aerodynamics/constants.py
@@ -1,5 +1,6 @@
+"""Constants for aerodynamics models."""
 #  This file is part of FAST-OAD : A framework for rapid Overall Aircraft Design
-#  Copyright (C) 2020  ONERA & ISAE-SUPAERO
+#  Copyright (C) 2021 ONERA & ISAE-SUPAERO
 #  FAST is free software: you can redistribute it and/or modify
 #  it under the terms of the GNU General Public License as published by
 #  the Free Software Foundation, either version 3 of the License, or
@@ -11,4 +12,36 @@
 #  You should have received a copy of the GNU General Public License
 #  along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
+from aenum import Enum
+
 POLAR_POINT_COUNT = 150
+
+SERVICE_OSWALD_COEFFICIENT = "service.aerodynamics.oswald_coefficient"
+SERVICE_INDUCED_DRAG_COEFFICIENT = "service.aerodynamics.induced_drag_coefficient"
+SERVICE_REYNOLDS_COEFFICIENT = "service.aerodynamics.reynolds_coefficient"
+SERVICE_INITIALIZE_CL = "service.aerodynamics.initialize_CL"
+SERVICE_CD0 = "service.aerodynamics.CD0"
+SERVICE_CD0_WING = "service.aerodynamics.CD0.wing"
+SERVICE_CD0_FUSELAGE = "service.aerodynamics.CD0.fuselage"
+SERVICE_CD0_HORIZONTAL_TAIL = "service.aerodynamics.CD0.horizontal_tail"
+SERVICE_CD0_VERTICAL_TAIL = "service.aerodynamics.CD0.vertical_tail"
+SERVICE_CD0_NACELLES_PYLONS = "service.aerodynamics.CD0.nacelles_pylons"
+SERVICE_CD0_SUM = "service.aerodynamics.CD0.sum"
+SERVICE_CD_COMPRESSIBILITY = "service.aerodynamics.CD.compressibility"
+SERVICE_CD_TRIM = "service.aerodynamics.CD.trim"
+SERVICE_POLAR = "service.aerodynamics.polar"
+SERVICE_HIGH_LIFT = "service.aerodynamics.high_lift"
+SERVICE_XFOIL = "service.aerodynamics.xfoil"
+SERVICE_LANDING_MAX_CL_CLEAN = "service.aerodynamics.landing.max_CL_clean"
+SERVICE_LANDING_MAX_CL = "service.aerodynamics.landing.max_CL"
+SERVICE_LANDING_MACH_REYNOLDS = "service.aerodynamics.landing.mach_reynolds"
+SERVICE_LOW_SPEED_CL_AOA = "service.aerodynamics.low_speed.CL_AoA"
+
+
+class PolarType(Enum):
+    """Enumeration of polar types to be computed."""
+
+    HIGH_SPEED = "high_speed"
+    LOW_SPEED = "low_speed"
+    TAKEOFF = "takeoff"
+    LANDING = "landing"
diff --git a/src/fastoad/models/aerodynamics/external/xfoil/__init__.py b/src/fastoad/models/aerodynamics/external/xfoil/__init__.py
index 932e2e327..e15826c5d 100644
--- a/src/fastoad/models/aerodynamics/external/xfoil/__init__.py
+++ b/src/fastoad/models/aerodynamics/external/xfoil/__init__.py
@@ -1,6 +1,6 @@
 """ Module for OpenMDAO-embedded XFOIL """
 #  This file is part of FAST-OAD : A framework for rapid Overall Aircraft Design
-#  Copyright (C) 2020  ONERA & ISAE-SUPAERO
+#  Copyright (C) 2021 ONERA & ISAE-SUPAERO
 #  FAST is free software: you can redistribute it and/or modify
 #  it under the terms of the GNU General Public License as published by
 #  the Free Software Foundation, either version 3 of the License, or
@@ -11,5 +11,3 @@
 #  GNU General Public License for more details.
 #  You should have received a copy of the GNU General Public License
 #  along with this program.  If not, see <https://www.gnu.org/licenses/>.
-
-from .xfoil_polar import XfoilPolar
diff --git a/src/fastoad/models/aerodynamics/external/xfoil/xfoil_polar.py b/src/fastoad/models/aerodynamics/external/xfoil/xfoil_polar.py
index 266466287..1742d005e 100644
--- a/src/fastoad/models/aerodynamics/external/xfoil/xfoil_polar.py
+++ b/src/fastoad/models/aerodynamics/external/xfoil/xfoil_polar.py
@@ -23,13 +23,15 @@
 from tempfile import TemporaryDirectory
 
 import numpy as np
-from openmdao.components.external_code_comp import ExternalCodeComp
+import openmdao.api as om
 from openmdao.utils.file_wrap import InputFileGenerator
 
 from fastoad._utils.resource_management.copy import copy_resource
 from fastoad.models.aerodynamics.external.xfoil import xfoil699
-from fastoad.models.geometry.profiles import get_profile
+from fastoad.models.geometry.profiles.profile_getter import get_profile
+from fastoad.module_management.service_registry import RegisterSubmodel
 from . import resources
+from ...constants import SERVICE_XFOIL
 
 OPTION_RESULT_POLAR_FILENAME = "result_polar_filename"
 OPTION_RESULT_FOLDER_PATH = "result_folder_path"
@@ -55,7 +57,8 @@
 _XFOIL_PATH_LIMIT = 64
 
 
-class XfoilPolar(ExternalCodeComp):
+@RegisterSubmodel(SERVICE_XFOIL, "fastoad.submodel.aerodynamics.xfoil")
+class XfoilPolar(om.ExternalCodeComp):
     """
     Runs a polar computation with XFOIL and returns the 2D max lift coefficient
     """
diff --git a/src/fastoad/models/geometry/profiles/__init__.py b/src/fastoad/models/geometry/profiles/__init__.py
index ddc128117..0c52195cc 100644
--- a/src/fastoad/models/geometry/profiles/__init__.py
+++ b/src/fastoad/models/geometry/profiles/__init__.py
@@ -2,7 +2,7 @@
 Management of wing profiles
 """
 #  This file is part of FAST-OAD : A framework for rapid Overall Aircraft Design
-#  Copyright (C) 2020  ONERA & ISAE-SUPAERO
+#  Copyright (C) 2021 ONERA & ISAE-SUPAERO
 #  FAST is free software: you can redistribute it and/or modify
 #  it under the terms of the GNU General Public License as published by
 #  the Free Software Foundation, either version 3 of the License, or
@@ -13,5 +13,3 @@
 #  GNU General Public License for more details.
 #  You should have received a copy of the GNU General Public License
 #  along with this program.  If not, see <https://www.gnu.org/licenses/>.
-
-from .get_profile import get_profile
diff --git a/src/fastoad/models/geometry/profiles/get_profile.py b/src/fastoad/models/geometry/profiles/profile_getter.py
similarity index 97%
rename from src/fastoad/models/geometry/profiles/get_profile.py
rename to src/fastoad/models/geometry/profiles/profile_getter.py
index 27a19b9f3..fd671c75b 100644
--- a/src/fastoad/models/geometry/profiles/get_profile.py
+++ b/src/fastoad/models/geometry/profiles/profile_getter.py
@@ -2,7 +2,7 @@
 Airfoil reshape function
 """
 #  This file is part of FAST-OAD : A framework for rapid Overall Aircraft Design
-#  Copyright (C) 2020  ONERA & ISAE-SUPAERO
+#  Copyright (C) 2021 ONERA & ISAE-SUPAERO
 #  FAST is free software: you can redistribute it and/or modify
 #  it under the terms of the GNU General Public License as published by
 #  the Free Software Foundation, either version 3 of the License, or
diff --git a/src/fastoad/models/variable_descriptions.txt b/src/fastoad/models/variable_descriptions.txt
index 6ad084d75..49dacb7bc 100644
--- a/src/fastoad/models/variable_descriptions.txt
+++ b/src/fastoad/models/variable_descriptions.txt
@@ -8,15 +8,19 @@ data:TLAR:cruise_mach || top-level requirement: cruise Mach number
 data:TLAR:max_TOFL || top-level requirement: maximum takeoff field length
 data:TLAR:range || top-level requirement: design range
 data:aerodynamics:aircraft:cruise:CD || drag coefficient in cruise conditions w.r.t. data:aerodynamics:aircraft:cruise:CL
-data:aerodynamics:aircraft:cruise:CD0 || profile drag coefficient for whole aircraft w.r.t. data:aerodynamics:aircraft:cruise:CL
+data:aerodynamics:aircraft:cruise:CD0 || profile drag coefficient for whole aircraft in cruise conditions w.r.t. data:aerodynamics:aircraft:cruise:CL
+data:aerodynamics:aircraft:cruise:CD0:clean || like data:aerodynamics:aircraft:cruise:CD0 but without parasitic drag
+data:aerodynamics:aircraft:cruise:CD0:parasitic || estimated parasitic drag for whole aircraft in cruise conditions (no high-lift) w.r.t. data:aerodynamics:aircraft:low_speed:CL
 data:aerodynamics:aircraft:cruise:CD:compressibility || increment of drag coefficient due to compressibility effects (shock waves) w.r.t. data:aerodynamics:aircraft:cruise:CL
 data:aerodynamics:aircraft:cruise:CD:trim || increment of drag coefficient due to aircraft trim w.r.t. data:aerodynamics:aircraft:cruise:CL
 data:aerodynamics:aircraft:cruise:CL || scale of lift coefficient values for drag computations in cruise conditions
+data:aerodynamics:aircraft:cruise:CL || scale of lift coefficient values for drag computations in cruise conditions
 data:aerodynamics:aircraft:cruise:CL_alpha || derivative of lift coefficient with respect to angle of attack in cruise conditions
 data:aerodynamics:aircraft:cruise:L_D_max || max lift/drag ratio in cruise conditions
 data:aerodynamics:aircraft:cruise:induced_drag_coefficient || multiply squared lift coefficient by this coefficient to get induced drag coefficient
 data:aerodynamics:aircraft:cruise:optimal_CD || drag coefficient at maximum lift/drag ratio in cruise conditions
 data:aerodynamics:aircraft:cruise:optimal_CL || lift coefficient at maximum lift/drag ratio in cruise conditions
+data:aerodynamics:aircraft:cruise:oswald_coefficient || Oswald coefficient for cruise conditions
 data:aerodynamics:aircraft:landing:CL_max || maximum lift coefficient in landing conditions
 data:aerodynamics:aircraft:landing:CL_max_clean || maximum lift coefficient in landing conditions without high-lift devices
 data:aerodynamics:aircraft:landing:CL_max_clean_2D || maximum lift coefficient of 2D average profile in landing conditions without high-lift devices
@@ -24,10 +28,13 @@ data:aerodynamics:aircraft:landing:additional_CL_capacity || at landing, is equa
 data:aerodynamics:aircraft:landing:mach || considered Mach number for landing phase
 data:aerodynamics:aircraft:low_speed:CD || drag coefficient in low speed conditions (no high-lift) w.r.t. data:aerodynamics:aircraft:low_speed:CL
 data:aerodynamics:aircraft:low_speed:CD0 || profile drag coefficient for whole aircraft in low speed conditions (no high-lift) w.r.t. data:aerodynamics:aircraft:low_speed:CL
+data:aerodynamics:aircraft:low_speed:CD0:clean || like data:aerodynamics:aircraft:low_speed:CD0 but without parasitic drag
+data:aerodynamics:aircraft:low_speed:CD0:parasitic || estimated parasitic drag for whole aircraft in low speed conditions (no high-lift) w.r.t. data:aerodynamics:aircraft:low_speed:CL
 data:aerodynamics:aircraft:low_speed:CD:trim || increment of drag coefficient due to aircraft trim in low speed conditions (no high-lift) w.r.t. data:aerodynamics:aircraft:low_speed:CL
 data:aerodynamics:aircraft:low_speed:CL || scale of lift coefficient values for drag computations in low speed conditions
 data:aerodynamics:aircraft:low_speed:CL_alpha || derivative of lift coefficient with respect to angle of attack in low speed conditions
 data:aerodynamics:aircraft:low_speed:induced_drag_coefficient || multiply squared lift coefficient by this coefficient to get induced drag coefficient
+data:aerodynamics:aircraft:low_speed:oswald_coefficient || Oswald coefficient for low speed conditions
 data:aerodynamics:aircraft:takeoff:CD || drag coefficient in takeoff conditions w.r.t. data:aerodynamics:aircraft:takeoff:CL
 data:aerodynamics:aircraft:takeoff:CL || scale of lift coefficient values for drag computations in takeoff conditions
 data:aerodynamics:aircraft:takeoff:CL0_clean || value of lift coefficient at angle of attack 0 in clean configuration at takeoff
diff --git a/src/fastoad/models/weight/cg/cg_components/compute_cg_tanks.py b/src/fastoad/models/weight/cg/cg_components/compute_cg_tanks.py
index 2c77f3738..d4ff602e6 100644
--- a/src/fastoad/models/weight/cg/cg_components/compute_cg_tanks.py
+++ b/src/fastoad/models/weight/cg/cg_components/compute_cg_tanks.py
@@ -20,7 +20,7 @@
 import openmdao.api as om
 from scipy import interpolate
 
-from fastoad.models.geometry.profiles import get_profile
+from fastoad.models.geometry.profiles.profile_getter import get_profile
 from fastoad.module_management.service_registry import RegisterSubmodel
 from ..constants import SERVICE_TANKS_CG
 


Test patch:
diff --git a/src/fastoad/models/aerodynamics/components/tests/test_components.py b/src/fastoad/models/aerodynamics/components/tests/test_components.py
index bf6e56361..1e48fa353 100644
--- a/src/fastoad/models/aerodynamics/components/tests/test_components.py
+++ b/src/fastoad/models/aerodynamics/components/tests/test_components.py
@@ -29,10 +29,11 @@
 from ..cd_compressibility import CdCompressibility
 from ..cd_trim import CdTrim
 from ..compute_low_speed_aero import ComputeAerodynamicsLowSpeed
-from ..compute_polar import ComputePolar, PolarType
+from ..compute_polar import ComputePolar
 from ..compute_reynolds import ComputeReynolds
 from ..high_lift_aero import ComputeDeltaHighLift
-from ..oswald import OswaldCoefficient
+from ..oswald import InducedDragCoefficient, OswaldCoefficient
+from ...constants import PolarType
 
 
 def get_indep_var_comp(var_names):
@@ -106,7 +107,7 @@ def get_cl_cd(slat_angle, flap_angle, mach, landing_flag, multi_slot_flag):
     assert cd == approx(0.02230, abs=1e-5)
 
 
-def test_oswald():
+def test_oswald_coefficient():
     """ Tests OswaldCoefficient """
     input_list = [
         "data:geometry:wing:area",
@@ -118,14 +119,46 @@ def test_oswald():
         "data:geometry:wing:sweep_25",
     ]
 
-    def get_coeff(mach):
+    def get_coeff(mach, low_speed_aero=False):
         ivc = get_indep_var_comp(input_list)
-        ivc.add_output("data:TLAR:cruise_mach", mach)
-        problem = run_system(OswaldCoefficient(), ivc)
-        return problem["data:aerodynamics:aircraft:cruise:induced_drag_coefficient"]
+        if low_speed_aero:
+            ivc.add_output("data:aerodynamics:aircraft:takeoff:mach", mach)
+        else:
+            ivc.add_output("data:TLAR:cruise_mach", mach)
+        problem = run_system(OswaldCoefficient(low_speed_aero=low_speed_aero), ivc)
+        if low_speed_aero:
+            return problem["data:aerodynamics:aircraft:low_speed:oswald_coefficient"]
+        else:
+            return problem["data:aerodynamics:aircraft:cruise:oswald_coefficient"]
+
+    assert get_coeff(0.2, True) == approx(0.8490, abs=1e-4)
+    # FIXME: Obtained Oswald coeff is too low, especially for high speed.
+    assert get_coeff(0.8, False) == approx(0.7451, abs=1e-4)
+
+
+def test_induced_drag_coefficient():
+    """ Tests InducedDragCoefficient """
+    input_list = [
+        "data:geometry:wing:area",
+        "data:geometry:wing:span",
+        "data:aerodynamics:aircraft:low_speed:oswald_coefficient",
+        "data:aerodynamics:aircraft:cruise:oswald_coefficient",
+    ]
+
+    def get_coeff(mach, low_speed_aero=False):
+        ivc = get_indep_var_comp(input_list)
+        if low_speed_aero:
+            ivc.add_output("data:aerodynamics:aircraft:low_speed:oswald_coefficient", mach)
+        else:
+            ivc.add_output("data:aerodynamics:aircraft:cruise:oswald_coefficient", mach)
+        problem = run_system(InducedDragCoefficient(low_speed_aero=low_speed_aero), ivc)
+        if low_speed_aero:
+            return problem["data:aerodynamics:aircraft:low_speed:induced_drag_coefficient"]
+        else:
+            return problem["data:aerodynamics:aircraft:cruise:induced_drag_coefficient"]
 
-    assert get_coeff(0.2) == approx(0.0465, abs=1e-4)
-    assert get_coeff(0.8) == approx(0.0530, abs=1e-4)
+    assert get_coeff(0.8490, True) == approx(0.0465, abs=1e-4)
+    assert get_coeff(0.7451, False) == approx(0.0530, abs=1e-4)
 
 
 def test_cd0():
@@ -288,6 +321,7 @@ def test_polar_high_speed():
     group = Group()
     group.add_subsystem("reynolds", ComputeReynolds(), promotes=["*"])
     group.add_subsystem("oswald", OswaldCoefficient(), promotes=["*"])
+    group.add_subsystem("induced_drag_coeff", InducedDragCoefficient(), promotes=["*"])
     group.add_subsystem("cd0", CD0(), promotes=["*"])
     group.add_subsystem("cd_compressibility", CdCompressibility(), promotes=["*"])
     group.add_subsystem("cd_trim", CdTrim(), promotes=["*"])
@@ -352,9 +386,12 @@ def test_polar_low_speed():
     group = Group()
     group.add_subsystem("reynolds", ComputeReynolds(low_speed_aero=True), promotes=["*"])
     group.add_subsystem("oswald", OswaldCoefficient(low_speed_aero=True), promotes=["*"])
+    group.add_subsystem(
+        "induced_drag_coeff", InducedDragCoefficient(low_speed_aero=True), promotes=["*"]
+    )
     group.add_subsystem("cd0", CD0(low_speed_aero=True), promotes=["*"])
     group.add_subsystem("cd_trim", CdTrim(low_speed_aero=True), promotes=["*"])
-    group.add_subsystem("polar", ComputePolar(type=PolarType.LOW_SPEED), promotes=["*"])
+    group.add_subsystem("polar", ComputePolar(polar_type=PolarType.LOW_SPEED), promotes=["*"])
 
     ivc = get_indep_var_comp(input_list)
     ivc.add_output("data:aerodynamics:aircraft:low_speed:CL", np.arange(0.0, 1.5, 0.01))
@@ -412,9 +449,12 @@ def test_polar_high_lift():
     group = Group()
     group.add_subsystem("reynolds", ComputeReynolds(low_speed_aero=True), promotes=["*"])
     group.add_subsystem("oswald", OswaldCoefficient(low_speed_aero=True), promotes=["*"])
+    group.add_subsystem(
+        "induced_drag_coeff", InducedDragCoefficient(low_speed_aero=True), promotes=["*"]
+    )
     group.add_subsystem("cd0", CD0(low_speed_aero=True), promotes=["*"])
     group.add_subsystem("cd_trim", CdTrim(low_speed_aero=True), promotes=["*"])
-    group.add_subsystem("polar", ComputePolar(type=PolarType.TAKEOFF), promotes=["*"])
+    group.add_subsystem("polar", ComputePolar(polar_type=PolarType.TAKEOFF), promotes=["*"])
 
     ivc = get_indep_var_comp(input_list)
     ivc.add_output("data:aerodynamics:aircraft:low_speed:CL", np.arange(0.0, 1.5, 0.01))
diff --git a/src/fastoad/models/aerodynamics/tests/test_aerodynamics.py b/src/fastoad/models/aerodynamics/tests/test_aerodynamics.py
deleted file mode 100644
index b0b70a326..000000000
--- a/src/fastoad/models/aerodynamics/tests/test_aerodynamics.py
+++ /dev/null
@@ -1,224 +0,0 @@
-"""
-Test module for aerodynamics groups
-"""
-#  This file is part of FAST-OAD : A framework for rapid Overall Aircraft Design
-#  Copyright (C) 2020  ONERA & ISAE-SUPAERO
-#  FAST is free software: you can redistribute it and/or modify
-#  it under the terms of the GNU General Public License as published by
-#  the Free Software Foundation, either version 3 of the License, or
-#  (at your option) any later version.
-#  This program is distributed in the hope that it will be useful,
-#  but WITHOUT ANY WARRANTY; without even the implied warranty of
-#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-#  GNU General Public License for more details.
-#  You should have received a copy of the GNU General Public License
-#  along with this program.  If not, see <https://www.gnu.org/licenses/>.
-
-import os.path as pth
-from platform import system
-
-import pytest
-from pytest import approx
-
-from fastoad.io import VariableIO
-from tests.testing_utilities import run_system
-from tests.xfoil_exe.get_xfoil import get_xfoil_path
-from ..aerodynamics_high_speed import AerodynamicsHighSpeed
-from ..aerodynamics_landing import AerodynamicsLanding
-from ..aerodynamics_low_speed import AerodynamicsLowSpeed
-
-xfoil_path = None if system() == "Windows" else get_xfoil_path()
-
-
-def get_indep_var_comp(var_names):
-    """ Reads required input data and returns an IndepVarcomp() instance"""
-    reader = VariableIO(pth.join(pth.dirname(__file__), "data", "aerodynamics_inputs.xml"))
-    reader.path_separator = ":"
-    ivc = reader.read(only=var_names).to_ivc()
-    return ivc
-
-
-@pytest.mark.skipif(
-    system() != "Windows" and xfoil_path is None, reason="No XFOIL executable available"
-)
-def test_aerodynamics_landing_with_xfoil():
-    """ Tests AerodynamicsHighSpeed """
-    input_list = [
-        "data:TLAR:approach_speed",
-        "data:mission:sizing:landing:flap_angle",
-        "data:mission:sizing:landing:slat_angle",
-        "data:geometry:wing:MAC:length",
-        "data:geometry:wing:thickness_ratio",
-        "data:geometry:wing:sweep_25",
-        "data:geometry:wing:sweep_0",
-        "data:geometry:wing:sweep_100_outer",
-        "data:geometry:flap:chord_ratio",
-        "data:geometry:flap:span_ratio",
-        "data:geometry:slat:chord_ratio",
-        "data:geometry:slat:span_ratio",
-        "xfoil:mach",
-        "tuning:aerodynamics:aircraft:landing:CL_max:landing_gear_effect:k",
-    ]
-
-    ivc = get_indep_var_comp(input_list)
-    problem = run_system(AerodynamicsLanding(xfoil_exe_path=xfoil_path), ivc)
-    # Reference values are for Windows XFOIL version, but as results can be slightly different on other
-    # platforms, tolerance is extended to 1e-2
-    assert problem["data:aerodynamics:aircraft:landing:CL_max_clean"] == approx(1.59359, abs=1e-2)
-    assert problem["data:aerodynamics:aircraft:landing:CL_max"] == approx(2.82178, abs=1e-2)
-
-
-def test_aerodynamics_landing_without_xfoil():
-    """ Tests AerodynamicsHighSpeed """
-    input_list = [
-        "data:TLAR:approach_speed",
-        "data:mission:sizing:landing:flap_angle",
-        "data:mission:sizing:landing:slat_angle",
-        "data:geometry:wing:MAC:length",
-        "data:geometry:wing:thickness_ratio",
-        "data:geometry:wing:sweep_25",
-        "data:geometry:wing:sweep_0",
-        "data:geometry:wing:sweep_100_outer",
-        "data:geometry:flap:chord_ratio",
-        "data:geometry:flap:span_ratio",
-        "data:geometry:slat:chord_ratio",
-        "data:geometry:slat:span_ratio",
-        "xfoil:mach",
-        "tuning:aerodynamics:aircraft:landing:CL_max:landing_gear_effect:k",
-        "data:aerodynamics:aircraft:landing:CL_max_clean_2D",
-    ]
-
-    ivc = get_indep_var_comp(input_list)
-    problem = run_system(AerodynamicsLanding(use_xfoil=False), ivc)
-    assert problem["data:aerodynamics:aircraft:landing:CL_max_clean"] == approx(1.54978, abs=1e-5)
-    assert problem["data:aerodynamics:aircraft:landing:CL_max"] == approx(2.77798, abs=1e-5)
-
-
-def test_aerodynamics_high_speed():
-    """ Tests AerodynamicsHighSpeed """
-    input_list = [
-        "data:geometry:propulsion:engine:count",
-        "data:geometry:propulsion:fan:length",
-        "data:geometry:flap:chord_ratio",
-        "data:geometry:flap:span_ratio",
-        "data:geometry:fuselage:maximum_height",
-        "data:geometry:fuselage:length",
-        "data:geometry:fuselage:wetted_area",
-        "data:geometry:fuselage:maximum_width",
-        "data:geometry:horizontal_tail:MAC:length",
-        "data:geometry:horizontal_tail:sweep_25",
-        "data:geometry:horizontal_tail:thickness_ratio",
-        "data:geometry:horizontal_tail:wetted_area",
-        "data:geometry:propulsion:nacelle:length",
-        "data:geometry:propulsion:nacelle:wetted_area",
-        "data:geometry:propulsion:pylon:length",
-        "data:geometry:propulsion:pylon:wetted_area",
-        "data:geometry:aircraft:wetted_area",
-        "data:geometry:slat:chord_ratio",
-        "data:geometry:slat:span_ratio",
-        "data:geometry:vertical_tail:MAC:length",
-        "data:geometry:vertical_tail:sweep_25",
-        "data:geometry:vertical_tail:thickness_ratio",
-        "data:geometry:vertical_tail:wetted_area",
-        "data:geometry:wing:area",
-        "data:geometry:wing:MAC:length",
-        "data:geometry:wing:root:chord",
-        "data:geometry:wing:tip:chord",
-        "data:geometry:wing:span",
-        "data:geometry:wing:sweep_0",
-        "data:geometry:wing:sweep_100_outer",
-        "data:geometry:wing:sweep_25",
-        "data:geometry:wing:thickness_ratio",
-        "data:geometry:wing:wetted_area",
-        "tuning:aerodynamics:aircraft:cruise:CD:k",
-        "tuning:aerodynamics:aircraft:cruise:CL:k",
-        "tuning:aerodynamics:aircraft:landing:CL_max:landing_gear_effect:k",
-        "tuning:aerodynamics:aircraft:cruise:CD:winglet_effect:k",
-        "tuning:aerodynamics:aircraft:cruise:CL:winglet_effect:k",
-        "tuning:aerodynamics:aircraft:cruise:CD:offset",
-        "tuning:aerodynamics:aircraft:cruise:CL:offset",
-        "tuning:aerodynamics:aircraft:cruise:CD:winglet_effect:offset",
-        "tuning:aerodynamics:aircraft:cruise:CL:winglet_effect:offset",
-        "data:mission:sizing:main_route:cruise:altitude",
-        "data:mission:sizing:landing:flap_angle",
-        "data:mission:sizing:landing:slat_angle",
-        "data:TLAR:cruise_mach",
-        "data:TLAR:approach_speed",
-    ]
-
-    ivc = get_indep_var_comp(input_list)
-    problem = run_system(AerodynamicsHighSpeed(), ivc)
-
-    cd = problem["data:aerodynamics:aircraft:cruise:CD"]
-    cl = problem["data:aerodynamics:aircraft:cruise:CL"]
-
-    assert cd[cl == 0.0] == approx(0.02102, abs=1e-5)
-    assert cd[cl == 0.2] == approx(0.02282, abs=1e-5)
-    assert cd[cl == 0.42] == approx(0.02977, abs=1e-5)
-    assert cd[cl == 0.85] == approx(0.24041, abs=1e-5)
-
-    assert problem["data:aerodynamics:aircraft:cruise:optimal_CL"] == approx(0.51, abs=1e-3)
-    assert problem["data:aerodynamics:aircraft:cruise:optimal_CD"] == approx(0.03487, abs=1e-5)
-
-
-def test_aerodynamics_low_speed():
-    """ Tests AerodynamicsLowSpeed """
-    input_list = [
-        "data:geometry:propulsion:engine:count",
-        "data:geometry:propulsion:fan:length",
-        "data:geometry:flap:chord_ratio",
-        "data:geometry:flap:span_ratio",
-        "data:geometry:fuselage:maximum_height",
-        "data:geometry:fuselage:length",
-        "data:geometry:fuselage:wetted_area",
-        "data:geometry:fuselage:maximum_width",
-        "data:geometry:horizontal_tail:MAC:length",
-        "data:geometry:horizontal_tail:sweep_25",
-        "data:geometry:horizontal_tail:thickness_ratio",
-        "data:geometry:horizontal_tail:wetted_area",
-        "data:geometry:propulsion:nacelle:length",
-        "data:geometry:propulsion:nacelle:wetted_area",
-        "data:geometry:propulsion:pylon:length",
-        "data:geometry:propulsion:pylon:wetted_area",
-        "data:geometry:aircraft:wetted_area",
-        "data:geometry:slat:chord_ratio",
-        "data:geometry:slat:span_ratio",
-        "data:geometry:vertical_tail:MAC:length",
-        "data:geometry:vertical_tail:sweep_25",
-        "data:geometry:vertical_tail:thickness_ratio",
-        "data:geometry:vertical_tail:wetted_area",
-        "data:geometry:wing:area",
-        "data:geometry:wing:aspect_ratio",
-        "data:geometry:wing:MAC:length",
-        "data:geometry:wing:root:chord",
-        "data:geometry:wing:tip:chord",
-        "data:geometry:wing:span",
-        "data:geometry:wing:sweep_0",
-        "data:geometry:wing:sweep_100_outer",
-        "data:geometry:wing:sweep_25",
-        "data:geometry:wing:thickness_ratio",
-        "data:geometry:wing:tip:thickness_ratio",
-        "data:geometry:wing:wetted_area",
-        "tuning:aerodynamics:aircraft:cruise:CD:k",
-        "tuning:aerodynamics:aircraft:cruise:CL:k",
-        "tuning:aerodynamics:aircraft:cruise:CD:winglet_effect:k",
-        "tuning:aerodynamics:aircraft:cruise:CL:winglet_effect:k",
-        "tuning:aerodynamics:aircraft:cruise:CD:offset",
-        "tuning:aerodynamics:aircraft:cruise:CL:offset",
-        "tuning:aerodynamics:aircraft:cruise:CD:winglet_effect:offset",
-        "tuning:aerodynamics:aircraft:cruise:CL:winglet_effect:offset",
-        "data:mission:sizing:landing:flap_angle",
-        "data:mission:sizing:landing:slat_angle",
-        "data:TLAR:approach_speed",
-    ]
-
-    ivc = get_indep_var_comp(input_list)
-    problem = run_system(AerodynamicsLowSpeed(), ivc)
-
-    cd = problem["data:aerodynamics:aircraft:low_speed:CD"]
-    cl = problem["data:aerodynamics:aircraft:low_speed:CL"]
-
-    assert cd[cl == 0.0] == approx(0.02173, abs=1e-5)
-    assert cd[cl == 0.2] == approx(0.02339, abs=1e-5)
-    assert cd[cl == 0.42] == approx(0.02974, abs=1e-5)
-    assert cd[cl == 0.85] == approx(0.06010, abs=1e-5)
diff --git a/src/fastoad/models/aerodynamics/tests/test_aerodynamics_landing.py b/src/fastoad/models/aerodynamics/tests/test_aerodynamics_landing.py
new file mode 100644
index 000000000..03b2a04e5
--- /dev/null
+++ b/src/fastoad/models/aerodynamics/tests/test_aerodynamics_landing.py
@@ -0,0 +1,92 @@
+"""
+Test module for aerodynamics groups
+"""
+#  This file is part of FAST-OAD : A framework for rapid Overall Aircraft Design
+#  Copyright (C) 2021 ONERA & ISAE-SUPAERO
+#  FAST is free software: you can redistribute it and/or modify
+#  it under the terms of the GNU General Public License as published by
+#  the Free Software Foundation, either version 3 of the License, or
+#  (at your option) any later version.
+#  This program is distributed in the hope that it will be useful,
+#  but WITHOUT ANY WARRANTY; without even the implied warranty of
+#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#  GNU General Public License for more details.
+#  You should have received a copy of the GNU General Public License
+#  along with this program.  If not, see <https://www.gnu.org/licenses/>.
+
+import os.path as pth
+from platform import system
+
+import pytest
+from pytest import approx
+
+from fastoad.io import VariableIO
+from tests.testing_utilities import run_system
+from tests.xfoil_exe.get_xfoil import get_xfoil_path
+from ..aerodynamics_landing import AerodynamicsLanding
+
+xfoil_path = None if system() == "Windows" else get_xfoil_path()
+
+
+def get_indep_var_comp(var_names):
+    """ Reads required input data and returns an IndepVarcomp() instance"""
+    reader = VariableIO(pth.join(pth.dirname(__file__), "data", "aerodynamics_inputs.xml"))
+    reader.path_separator = ":"
+    ivc = reader.read(only=var_names).to_ivc()
+    return ivc
+
+
+@pytest.mark.skipif(
+    system() != "Windows" and xfoil_path is None, reason="No XFOIL executable available"
+)
+def test_aerodynamics_landing_with_xfoil():
+    """ Tests AerodynamicsHighSpeed """
+    input_list = [
+        "data:TLAR:approach_speed",
+        "data:mission:sizing:landing:flap_angle",
+        "data:mission:sizing:landing:slat_angle",
+        "data:geometry:wing:MAC:length",
+        "data:geometry:wing:thickness_ratio",
+        "data:geometry:wing:sweep_25",
+        "data:geometry:wing:sweep_0",
+        "data:geometry:wing:sweep_100_outer",
+        "data:geometry:flap:chord_ratio",
+        "data:geometry:flap:span_ratio",
+        "data:geometry:slat:chord_ratio",
+        "data:geometry:slat:span_ratio",
+        "xfoil:mach",
+        "tuning:aerodynamics:aircraft:landing:CL_max:landing_gear_effect:k",
+    ]
+
+    ivc = get_indep_var_comp(input_list)
+    problem = run_system(AerodynamicsLanding(xfoil_exe_path=xfoil_path), ivc)
+    # Reference values are for Windows XFOIL version, but as results can be slightly different on other
+    # platforms, tolerance is extended to 1e-2
+    assert problem["data:aerodynamics:aircraft:landing:CL_max_clean"] == approx(1.59359, abs=1e-2)
+    assert problem["data:aerodynamics:aircraft:landing:CL_max"] == approx(2.82178, abs=1e-2)
+
+
+def test_aerodynamics_landing_without_xfoil():
+    """ Tests AerodynamicsHighSpeed """
+    input_list = [
+        "data:TLAR:approach_speed",
+        "data:mission:sizing:landing:flap_angle",
+        "data:mission:sizing:landing:slat_angle",
+        "data:geometry:wing:MAC:length",
+        "data:geometry:wing:thickness_ratio",
+        "data:geometry:wing:sweep_25",
+        "data:geometry:wing:sweep_0",
+        "data:geometry:wing:sweep_100_outer",
+        "data:geometry:flap:chord_ratio",
+        "data:geometry:flap:span_ratio",
+        "data:geometry:slat:chord_ratio",
+        "data:geometry:slat:span_ratio",
+        "xfoil:mach",
+        "tuning:aerodynamics:aircraft:landing:CL_max:landing_gear_effect:k",
+        "data:aerodynamics:aircraft:landing:CL_max_clean_2D",
+    ]
+
+    ivc = get_indep_var_comp(input_list)
+    problem = run_system(AerodynamicsLanding(use_xfoil=False), ivc)
+    assert problem["data:aerodynamics:aircraft:landing:CL_max_clean"] == approx(1.54978, abs=1e-5)
+    assert problem["data:aerodynamics:aircraft:landing:CL_max"] == approx(2.77798, abs=1e-5)
diff --git a/src/fastoad/models/geometry/profiles/tests/test_profile.py b/src/fastoad/models/geometry/profiles/tests/test_profile.py
index f782e290a..605e016a7 100644
--- a/src/fastoad/models/geometry/profiles/tests/test_profile.py
+++ b/src/fastoad/models/geometry/profiles/tests/test_profile.py
@@ -1,6 +1,6 @@
 """ Test module for profile.py """
 #  This file is part of FAST-OAD : A framework for rapid Overall Aircraft Design
-#  Copyright (C) 2020  ONERA & ISAE-SUPAERO
+#  Copyright (C) 2021 ONERA & ISAE-SUPAERO
 #  FAST is free software: you can redistribute it and/or modify
 #  it under the terms of the GNU General Public License as published by
 #  the Free Software Foundation, either version 3 of the License, or
@@ -12,13 +12,14 @@
 #  You should have received a copy of the GNU General Public License
 #  along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
+# pylint: disable=redefined-outer-name  # false positive on pytest fixtures
+
 import numpy as np
 import pytest
 from numpy.testing import assert_allclose
 
-# pylint: disable=redefined-outer-name  # false positive on pytest fixtures
-from ..get_profile import get_profile
 from ..profile import Profile
+from ..profile_getter import get_profile
 
 
 @pytest.fixture()