diff --git a/src/Air_EQ.m b/src/+dryair/assembleCoefficientMatrices.m
similarity index 97%
rename from src/Air_EQ.m
rename to src/+dryair/assembleCoefficientMatrices.m
index 31b02107..0f55e318 100644
--- a/src/Air_EQ.m
+++ b/src/+dryair/assembleCoefficientMatrices.m
@@ -1,78 +1,78 @@
-function [RHS, AirMatrices, SAVE] = Air_EQ(AirMatrices, Delt_t, P_g)
-
-    C1 = AirMatrices.C1
-    C2 = AirMatrices.C2
-    C3 = AirMatrices.C3
-    C4 = AirMatrices.C4
-    C4_a = AirMatrices.C4_a
-    C5 = AirMatrices.C5
-    C5_a = AirMatrices.C5_a
-    C6 = AirMatrices.C6
-    C7 = AirMatrices.C7
-
-    ModelSettings = io.getModelSettings();
-    n = ModelSettings.NN;
-
-    # Alias of SoilVariables
-    SV = SoilVariables;
-
-    if ModelSettings.Thmrlefc
-        RHS(1) = -C7(1) + (C3(1, 1) * P_g(1) + C3(1, 2) * P_g(2)) / Delt_t ...
-            - (C2(1, 1) / Delt_t + C5(1, 1)) * SV.TT(1) - (C2(1, 2) / Delt_t + C5(1, 2)) * SV.TT(2) ...
-            - (C1(1, 1) / Delt_t + C4(1, 1)) * SV.hh(1) - (C1(1, 2) / Delt_t + C4(1, 2)) * SV.hh(2) ...
-            + (C2(1, 1) / Delt_t) * SV.T(1) + (C2(1, 2) / Delt_t) * SV.T(2) ...
-            + (C1(1, 1) / Delt_t) * SV.h(1) + (C1(1, 2) / Delt_t) * SV.h(2);
-
-        for i = 2:ModelSettings.NL
-            ARG1 = C2(i - 1, 2) / Delt_t;
-            ARG2 = C2(i, 1) / Delt_t;
-            ARG3 = C2(i, 2) / Delt_t;
-
-            ARG4 = C1(i - 1, 2) / Delt_t;
-            ARG5 = C1(i, 1) / Delt_t;
-            ARG6 = C1(i, 2) / Delt_t;
-
-            RHS(i) = -C7(i) + (C3(i - 1, 2) * P_g(i - 1) + C3(i, 1) * P_g(i) + C3(i, 2) * P_g(i + 1)) / Delt_t ...
-                - (ARG1 + C5_a(i - 1)) * SV.TT(i - 1) - (ARG2 + C5(i, 1)) * SV.TT(i) - (ARG3 + C5(i, 2)) * SV.TT(i + 1) ...
-                - (ARG4 + C4_a(i - 1)) * SV.hh(i - 1) - (ARG5 + C4(i, 1)) * SV.hh(i) - (ARG6 + C4(i, 2)) * SV.hh(i + 1) ...
-                + ARG1 * SV.T(i - 1) + ARG2 * SV.T(i) + ARG3 * SV.T(i + 1) ...
-                + ARG4 * SV.h(i - 1) + ARG5 * SV.h(i) + ARG6 * SV.h(i + 1);
-        end
-
-        RHS(n) = -C7(n) + (C3(n - 1, 2) * P_g(n - 1) + C3(n, 1) * P_g(n)) / Delt_t ...
-            - (C2(n - 1, 2) / Delt_t + C5_a(n - 1)) * SV.TT(n - 1) - (C2(n, 1) / Delt_t + C5(n, 1)) * SV.TT(n) ...
-            - (C1(n - 1, 2) / Delt_t + C4_a(n - 1)) * SV.hh(n - 1) - (C1(n, 1) / Delt_t + C4(n, 1)) * SV.hh(n) ...
-            + (C2(n - 1, 2) / Delt_t) * SV.T(n - 1) + (C2(n, 1) / Delt_t) * SV.T(n) ...
-            + (C1(n - 1, 2) / Delt_t) * SV.h(n - 1) + (C1(n, 1) / Delt_t) * SV.h(n);
-    else
-        ARG4 = C1(i - 1, 2) / Delt_t;
-        ARG5 = C1(i, 1) / Delt_t;
-        ARG6 = C1(i, 2) / Delt_t;
-
-        RHS(1) = -C7(1) + (C3(1, 1) * P_g(1) + C3(1, 2) * P_g(2)) / Delt_t ...
-            - (C1(1, 1) / Delt_t + C4(1, 1)) * SV.hh(1) - (C1(1, 2) / Delt_t + C4(1, 2)) * SV.hh(2) ...
-            + (C1(1, 1) / Delt_t) * SV.h(1) + (C1(1, 2) / Delt_t) * SV.h(2);
-        for i = 2:ModelSettings.NL
-            RHS(i) = -C7(i) + (C3(i - 1, 2) * P_g(i - 1) + C3(i, 1) * P_g(i) + C3(i, 2) * P_g(i + 1)) / Delt_t ...
-                - (ARG4 + C4(i - 1, 2)) * SV.hh(i - 1) - (ARG5 + C4(i, 1)) * SV.hh(i) - (ARG6 + C4(i, 2)) * SV.hh(i + 1) ...
-                + ARG4 * SV.h(i - 1) + ARG5 * SV.h(i) + ARG6 * SV.h(i + 1);
-        end
-        RHS(n) = -C7(n) + (C3(n - 1, 2) * P_g(n - 1) + C3(n, 1) * P_g(n)) / Delt_t ...
-            - (C1(n - 1, 2) / Delt_t + C4(n - 1, 2)) * SV.hh(n - 1) - (C1(n, 1) / Delt_t + C4(n, 1)) * SV.hh(n) ...
-            + (C1(n - 1, 2) / Delt_t) * SV.h(n - 1) + (C1(n, 1) / Delt_t) * SV.h(n);
-    end
-
-    for i = 1:ModelSettings.NN
-        for j = 1:ModelSettings.nD
-            C6(i, j) = C3(i, j) / Delt_t + C6(i, j);
-        end
-    end
-
-    AirMatrices.C6 = C6;
-
-    SAVE(1, 1, 3) = RHS(1);
-    SAVE(1, 2, 3) = C6(1, 1);
-    SAVE(1, 3, 3) = C6(1, 2);
-    SAVE(2, 1, 3) = RHS(n);
-    SAVE(2, 2, 3) = C6(n - 1, 2);
-    SAVE(2, 3, 3) = C6(n, 1);
+function [RHS, AirMatrices, SAVE] = Air_EQ(AirMatrices, Delt_t, P_g)
+
+    C1 = AirMatrices.C1
+    C2 = AirMatrices.C2
+    C3 = AirMatrices.C3
+    C4 = AirMatrices.C4
+    C4_a = AirMatrices.C4_a
+    C5 = AirMatrices.C5
+    C5_a = AirMatrices.C5_a
+    C6 = AirMatrices.C6
+    C7 = AirMatrices.C7
+
+    ModelSettings = io.getModelSettings();
+    n = ModelSettings.NN;
+
+    # Alias of SoilVariables
+    SV = SoilVariables;
+
+    if ModelSettings.Thmrlefc
+        RHS(1) = -C7(1) + (C3(1, 1) * P_g(1) + C3(1, 2) * P_g(2)) / Delt_t ...
+            - (C2(1, 1) / Delt_t + C5(1, 1)) * SV.TT(1) - (C2(1, 2) / Delt_t + C5(1, 2)) * SV.TT(2) ...
+            - (C1(1, 1) / Delt_t + C4(1, 1)) * SV.hh(1) - (C1(1, 2) / Delt_t + C4(1, 2)) * SV.hh(2) ...
+            + (C2(1, 1) / Delt_t) * SV.T(1) + (C2(1, 2) / Delt_t) * SV.T(2) ...
+            + (C1(1, 1) / Delt_t) * SV.h(1) + (C1(1, 2) / Delt_t) * SV.h(2);
+
+        for i = 2:ModelSettings.NL
+            ARG1 = C2(i - 1, 2) / Delt_t;
+            ARG2 = C2(i, 1) / Delt_t;
+            ARG3 = C2(i, 2) / Delt_t;
+
+            ARG4 = C1(i - 1, 2) / Delt_t;
+            ARG5 = C1(i, 1) / Delt_t;
+            ARG6 = C1(i, 2) / Delt_t;
+
+            RHS(i) = -C7(i) + (C3(i - 1, 2) * P_g(i - 1) + C3(i, 1) * P_g(i) + C3(i, 2) * P_g(i + 1)) / Delt_t ...
+                - (ARG1 + C5_a(i - 1)) * SV.TT(i - 1) - (ARG2 + C5(i, 1)) * SV.TT(i) - (ARG3 + C5(i, 2)) * SV.TT(i + 1) ...
+                - (ARG4 + C4_a(i - 1)) * SV.hh(i - 1) - (ARG5 + C4(i, 1)) * SV.hh(i) - (ARG6 + C4(i, 2)) * SV.hh(i + 1) ...
+                + ARG1 * SV.T(i - 1) + ARG2 * SV.T(i) + ARG3 * SV.T(i + 1) ...
+                + ARG4 * SV.h(i - 1) + ARG5 * SV.h(i) + ARG6 * SV.h(i + 1);
+        end
+
+        RHS(n) = -C7(n) + (C3(n - 1, 2) * P_g(n - 1) + C3(n, 1) * P_g(n)) / Delt_t ...
+            - (C2(n - 1, 2) / Delt_t + C5_a(n - 1)) * SV.TT(n - 1) - (C2(n, 1) / Delt_t + C5(n, 1)) * SV.TT(n) ...
+            - (C1(n - 1, 2) / Delt_t + C4_a(n - 1)) * SV.hh(n - 1) - (C1(n, 1) / Delt_t + C4(n, 1)) * SV.hh(n) ...
+            + (C2(n - 1, 2) / Delt_t) * SV.T(n - 1) + (C2(n, 1) / Delt_t) * SV.T(n) ...
+            + (C1(n - 1, 2) / Delt_t) * SV.h(n - 1) + (C1(n, 1) / Delt_t) * SV.h(n);
+    else
+        ARG4 = C1(i - 1, 2) / Delt_t;
+        ARG5 = C1(i, 1) / Delt_t;
+        ARG6 = C1(i, 2) / Delt_t;
+
+        RHS(1) = -C7(1) + (C3(1, 1) * P_g(1) + C3(1, 2) * P_g(2)) / Delt_t ...
+            - (C1(1, 1) / Delt_t + C4(1, 1)) * SV.hh(1) - (C1(1, 2) / Delt_t + C4(1, 2)) * SV.hh(2) ...
+            + (C1(1, 1) / Delt_t) * SV.h(1) + (C1(1, 2) / Delt_t) * SV.h(2);
+        for i = 2:ModelSettings.NL
+            RHS(i) = -C7(i) + (C3(i - 1, 2) * P_g(i - 1) + C3(i, 1) * P_g(i) + C3(i, 2) * P_g(i + 1)) / Delt_t ...
+                - (ARG4 + C4(i - 1, 2)) * SV.hh(i - 1) - (ARG5 + C4(i, 1)) * SV.hh(i) - (ARG6 + C4(i, 2)) * SV.hh(i + 1) ...
+                + ARG4 * SV.h(i - 1) + ARG5 * SV.h(i) + ARG6 * SV.h(i + 1);
+        end
+        RHS(n) = -C7(n) + (C3(n - 1, 2) * P_g(n - 1) + C3(n, 1) * P_g(n)) / Delt_t ...
+            - (C1(n - 1, 2) / Delt_t + C4(n - 1, 2)) * SV.hh(n - 1) - (C1(n, 1) / Delt_t + C4(n, 1)) * SV.hh(n) ...
+            + (C1(n - 1, 2) / Delt_t) * SV.h(n - 1) + (C1(n, 1) / Delt_t) * SV.h(n);
+    end
+
+    for i = 1:ModelSettings.NN
+        for j = 1:ModelSettings.nD
+            C6(i, j) = C3(i, j) / Delt_t + C6(i, j);
+        end
+    end
+
+    AirMatrices.C6 = C6;
+
+    SAVE(1, 1, 3) = RHS(1);
+    SAVE(1, 2, 3) = C6(1, 1);
+    SAVE(1, 3, 3) = C6(1, 2);
+    SAVE(2, 1, 3) = RHS(n);
+    SAVE(2, 2, 3) = C6(n - 1, 2);
+    SAVE(2, 3, 3) = C6(n, 1);
diff --git a/src/Air_BC.m b/src/+dryair/calculateBoundaryConditions.m
similarity index 97%
rename from src/Air_BC.m
rename to src/+dryair/calculateBoundaryConditions.m
index 3913b95a..befdf99f 100644
--- a/src/Air_BC.m
+++ b/src/+dryair/calculateBoundaryConditions.m
@@ -1,32 +1,32 @@
-function [RHS, AirMatrices] = Air_BC(BoundaryCondition, AirMatrices, ForcingData, RHS, KT)
-
-    TopPg = 100 .* (ForcingData.Pg_msr);
-    ModelSettings = io.getModelSettings();
-    n = ModelSettings.NN;
-    % Apply the bottom boundary condition called for by NBCPB
-    if BoundaryCondition.NBCPB == 1  % Bounded bottom with the water table
-        RHS(1) = BtmPg;
-        AirMatrices.C6(1, 1) = 1;
-        RHS(2) = RHS(2) - AirMatrices.C6(1, 2) * RHS(1);
-        AirMatrices.C6(1, 2) = 0;
-        AirMatrices.C6_a(1) = 0;
-    elseif NBCPB == 2 % The soil air is allowed to escape from the bottom
-        RHS(1) = RHS(1) + BoundaryCondition.BCPB;
-    end
-
-    % Apply the surface boundary condition called by NBCP
-    if BoundaryCondition.NBCP == 1    % Ponded infiltration with Bonded bottom
-        RHS(n) = BtmPg;
-        AirMatrices.C6(n, 1) = 1;
-        RHS(n - 1) = RHS(n - 1) - AirMatrices.C6(n - 1, 2) * RHS(n);
-        AirMatrices.C6(n - 1, 2) = 0;
-        AirMatrices.C6_a(n - 1) = 0;
-    elseif NBCP == 2          % Specified flux on the surface
-        RHS(n) = RHS(n) - BoundaryCondition.BCP;
-    else
-        RHS(n) = TopPg(KT);
-        AirMatrices.C6(n, 1) = 1;
-        RHS(n - 1) = RHS(n - 1) - AirMatrices.C6(n - 1, 2) * RHS(n);
-        AirMatrices.C6(n - 1, 2) = 0;
-        AirMatrices.C6_a(n - 1) = 0;
-    end
+function [RHS, AirMatrices] = Air_BC(BoundaryCondition, AirMatrices, ForcingData, RHS, KT)
+
+    TopPg = 100 .* (ForcingData.Pg_msr);
+    ModelSettings = io.getModelSettings();
+    n = ModelSettings.NN;
+    % Apply the bottom boundary condition called for by NBCPB
+    if BoundaryCondition.NBCPB == 1  % Bounded bottom with the water table
+        RHS(1) = BtmPg;
+        AirMatrices.C6(1, 1) = 1;
+        RHS(2) = RHS(2) - AirMatrices.C6(1, 2) * RHS(1);
+        AirMatrices.C6(1, 2) = 0;
+        AirMatrices.C6_a(1) = 0;
+    elseif NBCPB == 2 % The soil air is allowed to escape from the bottom
+        RHS(1) = RHS(1) + BoundaryCondition.BCPB;
+    end
+
+    % Apply the surface boundary condition called by NBCP
+    if BoundaryCondition.NBCP == 1    % Ponded infiltration with Bonded bottom
+        RHS(n) = BtmPg;
+        AirMatrices.C6(n, 1) = 1;
+        RHS(n - 1) = RHS(n - 1) - AirMatrices.C6(n - 1, 2) * RHS(n);
+        AirMatrices.C6(n - 1, 2) = 0;
+        AirMatrices.C6_a(n - 1) = 0;
+    elseif NBCP == 2          % Specified flux on the surface
+        RHS(n) = RHS(n) - BoundaryCondition.BCP;
+    else
+        RHS(n) = TopPg(KT);
+        AirMatrices.C6(n, 1) = 1;
+        RHS(n - 1) = RHS(n - 1) - AirMatrices.C6(n - 1, 2) * RHS(n);
+        AirMatrices.C6(n - 1, 2) = 0;
+        AirMatrices.C6_a(n - 1) = 0;
+    end
diff --git a/src/AirPARM.m b/src/+dryair/calculateDryAirParameters.m
similarity index 98%
rename from src/AirPARM.m
rename to src/+dryair/calculateDryAirParameters.m
index 481a1b32..ed0e5c8c 100644
--- a/src/AirPARM.m
+++ b/src/+dryair/calculateDryAirParameters.m
@@ -1,55 +1,55 @@
-function AirVariabes = AirPARM(SoilVariables, GasDispersivity, TransportCoefficient, InitialValues, GasDispersivity,...
-                               P_gg, Xah, XaT, Xaa, RHODA)
-
-    ModelSettings = io.getModelSettings();
-    Constants = io.define_constants();
-
-    AirVariabes.Cah = InitialValues.Cah;
-    AirVariabes.CaT = InitialValues.CaT;
-    AirVariabes.Caa = InitialValues.Caa;
-    AirVariabes.Kah = InitialValues.Kah;
-    AirVariabes.KaT = InitialValues.KaT;
-    AirVariabes.Kaa = InitialValues.Kaa;
-    AirVariabes.Vah = InitialValues.Vah;
-    AirVariabes.VaT = InitialValues.VaT;
-    AirVariabes.Vaa = InitialValues.Vaa;
-    AirVariabes.Cag = InitialValues.Cag;
-
-    for i = 1:ModelSettings.NL
-        for j = 1:ModelSettings.nD
-
-            KLhBAR = (SoilVariables.KfL_h(i, 1) + SoilVariables.KfL_h(i, 2)) / 2;
-            KLTBAR = (InitialValues.KL_T(i, 1) + InitialValues.KL_T(i, 2)) / 2;
-            DhDZ = (SoilVariables.hh(i + 1) + SoilVariables.hh_frez(i + 1) - SoilVariables.hh(i) - SoilVariables.hh_frez(i)) / ModelSettings.DeltZ(i);
-            DTDZ = (SoilVariables.TT(i + 1) - SoilVariables.TT(i)) / ModelSettings.DeltZ(i);
-            DPgDZ = (P_gg(i + 1) - P_gg(i)) / ModelSettings.DeltZ(i);
-            DTDBAR = (TransportCoefficient.D_Ta(i, 1) + TransportCoefficient.D_Ta(i, 2)) / 2;
-
-            if SoilVariables.KLa_Switch == 1
-                QL(i) = -(KLhBAR * (DhDZ + DPgDZ / Constants.Gamma_w) + (KLTBAR + DTDBAR) * DTDZ + KLhBAR);
-                QL_h(i) = -(KLhBAR * (DhDZ + DPgDZ / Constants.Gamma_w) + KLhBAR);
-                QL_a(i) = -(KLhBAR * (DPgDZ / Constants.Gamma_w));
-                QL_T(i) = -((KLTBAR + DTDBAR) * DTDZ);
-            else
-                QL(i) = -(KLhBAR * DhDZ + (KLTBAR + DTDBAR) * DTDZ + KLhBAR);
-                QL_h(i) = -(KLhBAR * DhDZ + KLhBAR);
-                QL_T(i) = -((KLTBAR + DTDBAR) * DTDZ);
-
-            end
-            MN = i + j - 1;
-
-            AirVariabes.Cah(i, j) = Xah(MN) * (SoilVariables.POR(i) + (Constants.Hc - 1) * SoilVariables.Theta_LL(i, j)) + (Constants.Hc - 1) * RHODA(MN) * SoilVariables.DTheta_LLh(i, j);
-            AirVariabes.CaT(i, j) = XaT(MN) * (SoilVariables.POR(i) + (Constants.Hc - 1) * SoilVariables.Theta_LL(i, j)) + (Constants.Hc - 1) * RHODA(MN) * SoilVariables.DTheta_LLT(i, j);
-            AirVariabes.Caa(i, j) = Xaa(MN) * (SoilVariables.POR(i) + (Constants.Hc - 1) * SoilVariables.Theta_LL(i, j));
-
-            AirVariabes.Kah(i, j) = Xah(MN) * (VaporVariables.D_V(i, j) + GasDispersivity.D_Vg(i)) + Constants.Hc * RHODA(MN) * SoilVariables.KfL_h(i, j);
-            AirVariabes.KaT(i, j) = XaT(MN) * (VaporVariables.D_V(i, j) + GasDispersivity.D_Vg(i)) + Constants.Hc * RHODA(MN) * (InitialValues.KL_T(i, j) + TransportCoefficient.D_Ta(i, j));
-            AirVariabes.Kaa(i, j) = Xaa(MN) * (VaporVariables.D_V(i, j) + GasDispersivity.D_Vg(i)) + RHODA(MN) * (GasDispersivity.Beta_g(i, j) + Constants.Hc * SoilVariables.KfL_h(i, j) / Constants.Gamma_w);
-
-            AirVariabes.Cag(i, j) = Constants.Hc * RHODA(MN) * SoilVariables.KfL_h(i, j);
-
-            AirVariabes.Vah(i, j) = -(GasDispersivity.V_A(i) + Constants.Hc * QL(i) / Constants.RHOL) * Xah(MN);
-            AirVariabes.VaT(i, j) = -(GasDispersivity.V_A(i) + Constants.Hc * QL(i) / Constants.RHOL) * XaT(MN);
-            AirVariabes.Vaa(i, j) = -(GasDispersivity.V_A(i) + Constants.Hc * QL(i) / Constants.RHOL) * Xaa(MN);
-        end
-    end
+function AirVariabes = AirPARM(SoilVariables, GasDispersivity, TransportCoefficient, InitialValues, GasDispersivity,...
+                               P_gg, Xah, XaT, Xaa, RHODA)
+
+    ModelSettings = io.getModelSettings();
+    Constants = io.define_constants();
+
+    AirVariabes.Cah = InitialValues.Cah;
+    AirVariabes.CaT = InitialValues.CaT;
+    AirVariabes.Caa = InitialValues.Caa;
+    AirVariabes.Kah = InitialValues.Kah;
+    AirVariabes.KaT = InitialValues.KaT;
+    AirVariabes.Kaa = InitialValues.Kaa;
+    AirVariabes.Vah = InitialValues.Vah;
+    AirVariabes.VaT = InitialValues.VaT;
+    AirVariabes.Vaa = InitialValues.Vaa;
+    AirVariabes.Cag = InitialValues.Cag;
+
+    for i = 1:ModelSettings.NL
+        for j = 1:ModelSettings.nD
+
+            KLhBAR = (SoilVariables.KfL_h(i, 1) + SoilVariables.KfL_h(i, 2)) / 2;
+            KLTBAR = (InitialValues.KL_T(i, 1) + InitialValues.KL_T(i, 2)) / 2;
+            DhDZ = (SoilVariables.hh(i + 1) + SoilVariables.hh_frez(i + 1) - SoilVariables.hh(i) - SoilVariables.hh_frez(i)) / ModelSettings.DeltZ(i);
+            DTDZ = (SoilVariables.TT(i + 1) - SoilVariables.TT(i)) / ModelSettings.DeltZ(i);
+            DPgDZ = (P_gg(i + 1) - P_gg(i)) / ModelSettings.DeltZ(i);
+            DTDBAR = (TransportCoefficient.D_Ta(i, 1) + TransportCoefficient.D_Ta(i, 2)) / 2;
+
+            if SoilVariables.KLa_Switch == 1
+                QL(i) = -(KLhBAR * (DhDZ + DPgDZ / Constants.Gamma_w) + (KLTBAR + DTDBAR) * DTDZ + KLhBAR);
+                QL_h(i) = -(KLhBAR * (DhDZ + DPgDZ / Constants.Gamma_w) + KLhBAR);
+                QL_a(i) = -(KLhBAR * (DPgDZ / Constants.Gamma_w));
+                QL_T(i) = -((KLTBAR + DTDBAR) * DTDZ);
+            else
+                QL(i) = -(KLhBAR * DhDZ + (KLTBAR + DTDBAR) * DTDZ + KLhBAR);
+                QL_h(i) = -(KLhBAR * DhDZ + KLhBAR);
+                QL_T(i) = -((KLTBAR + DTDBAR) * DTDZ);
+
+            end
+            MN = i + j - 1;
+
+            AirVariabes.Cah(i, j) = Xah(MN) * (SoilVariables.POR(i) + (Constants.Hc - 1) * SoilVariables.Theta_LL(i, j)) + (Constants.Hc - 1) * RHODA(MN) * SoilVariables.DTheta_LLh(i, j);
+            AirVariabes.CaT(i, j) = XaT(MN) * (SoilVariables.POR(i) + (Constants.Hc - 1) * SoilVariables.Theta_LL(i, j)) + (Constants.Hc - 1) * RHODA(MN) * SoilVariables.DTheta_LLT(i, j);
+            AirVariabes.Caa(i, j) = Xaa(MN) * (SoilVariables.POR(i) + (Constants.Hc - 1) * SoilVariables.Theta_LL(i, j));
+
+            AirVariabes.Kah(i, j) = Xah(MN) * (VaporVariables.D_V(i, j) + GasDispersivity.D_Vg(i)) + Constants.Hc * RHODA(MN) * SoilVariables.KfL_h(i, j);
+            AirVariabes.KaT(i, j) = XaT(MN) * (VaporVariables.D_V(i, j) + GasDispersivity.D_Vg(i)) + Constants.Hc * RHODA(MN) * (InitialValues.KL_T(i, j) + TransportCoefficient.D_Ta(i, j));
+            AirVariabes.Kaa(i, j) = Xaa(MN) * (VaporVariables.D_V(i, j) + GasDispersivity.D_Vg(i)) + RHODA(MN) * (GasDispersivity.Beta_g(i, j) + Constants.Hc * SoilVariables.KfL_h(i, j) / Constants.Gamma_w);
+
+            AirVariabes.Cag(i, j) = Constants.Hc * RHODA(MN) * SoilVariables.KfL_h(i, j);
+
+            AirVariabes.Vah(i, j) = -(GasDispersivity.V_A(i) + Constants.Hc * QL(i) / Constants.RHOL) * Xah(MN);
+            AirVariabes.VaT(i, j) = -(GasDispersivity.V_A(i) + Constants.Hc * QL(i) / Constants.RHOL) * XaT(MN);
+            AirVariabes.Vaa(i, j) = -(GasDispersivity.V_A(i) + Constants.Hc * QL(i) / Constants.RHOL) * Xaa(MN);
+        end
+    end
diff --git a/src/Air_MAT.m b/src/+dryair/calculateMatricCoefficients.m
similarity index 98%
rename from src/Air_MAT.m
rename to src/+dryair/calculateMatricCoefficients.m
index 90b9cc9c..23f19003 100644
--- a/src/Air_MAT.m
+++ b/src/+dryair/calculateMatricCoefficients.m
@@ -1,49 +1,49 @@
-function AirMatrices = Air_MAT(AirVariabes, InitialValues)
-
-    AirMatrices.C1 = InitialValues.C1;
-    AirMatrices.C2 = InitialValues.C2;
-    AirMatrices.C3 = InitialValues.C3;
-    AirMatrices.C4 = InitialValues.C4;
-    AirMatrices.C5 = InitialValues.C5;
-    AirMatrices.C6 = InitialValues.C6;
-
-    ModelSettings = io.getModelSettings();
-
-    for i = 1:ModelSettings.NL
-        AirMatrices.C1(i, 1) = AirMatrices.C1(i, 1) + AirVariabes.Cah(i, 1) * ModelSettings.DeltZ(i) / 2;
-        AirMatrices.C1(i + 1, 1) = AirMatrices.C1(i + 1, 1) + AirVariabes.Cah(i, 2) * ModelSettings.DeltZ(i) / 2;
-
-        AirMatrices.C2(i, 1) = AirMatrices.C2(i, 1) + AirVariabes.CaT(i, 1) * ModelSettings.DeltZ(i) / 2;
-        AirMatrices.C2(i + 1, 1) = AirMatrices.C2(i + 1, 1) + AirVariabes.CaT(i, 2) * ModelSettings.DeltZ(i) / 2;
-
-        AirMatrices.C3(i, 1) = AirMatrices.C3(i, 1) + AirVariabes.Caa(i, 1) * ModelSettings.DeltZ(i) / 2;
-        AirMatrices.C3(i + 1, 1) = AirMatrices.C3(i + 1, 1) + AirVariabes.Caa(i, 2) * ModelSettings.DeltZ(i) / 2;
-
-        C4ARG1 = (AirVariabes.Kah(i, 1) + AirVariabes.Kah(i, 2)) / (2 * ModelSettings.DeltZ(i));
-        C4ARG2_1 = AirVariabes.Vah(i, 1) / 3 + AirVariabes.Vah(i, 2) / 6;
-        C4ARG2_2 = AirVariabes.Vah(i, 1) / 6 + AirVariabes.Vah(i, 2) / 3;
-        AirMatrices.C4(i, 1) = AirMatrices.C4(i, 1) + C4ARG1 - C4ARG2_1;
-        AirMatrices.C4(i, 2) = AirMatrices.C4(i, 2) - C4ARG1 - C4ARG2_2;
-        AirMatrices.C4(i + 1, 1) = AirMatrices.C4(i + 1, 1) + C4ARG1 + C4ARG2_2;
-        AirMatrices.C4_a(i) = -C4ARG1 + C4ARG2_1;
-
-        C5ARG1 = (AirVariabes.KaT(i, 1) + AirVariabes.KaT(i, 2)) / (2 * ModelSettings.DeltZ(i));
-        C5ARG2_1 = AirVariabes.VaT(i, 1) / 3 + AirVariabes.VaT(i, 2) / 6;
-        C5ARG2_2 = AirVariabes.VaT(i, 1) / 6 + AirVariabes.VaT(i, 2) / 3;
-        AirMatrices.C5(i, 1) = AirMatrices.C5(i, 1) + C5ARG1 - C5ARG2_1;
-        AirMatrices.C5(i, 2) = AirMatrices.C5(i, 2) - C5ARG1 - C5ARG2_2;
-        AirMatrices.C5(i + 1, 1) = AirMatrices.C5(i + 1, 1) + C5ARG1 + C5ARG2_2;
-        AirMatrices.C5_a(i) = -C5ARG1 + C5ARG2_1;
-
-        C6ARG1 = (AirVariabes.Kaa(i, 1) + AirVariabes.Kaa(i, 2)) / (2 * ModelSettings.DeltZ(i));
-        C6ARG2_1 = AirVariabes.Vaa(i, 1) / 3 + AirVariabes.Vaa(i, 2) / 6;
-        C6ARG2_2 = AirVariabes.Vaa(i, 1) / 6 + AirVariabes.Vaa(i, 2) / 3;
-        AirMatrices.C6(i, 1) = AirMatrices.C6(i, 1) + C6ARG1 - C6ARG2_1;
-        AirMatrices.C6(i, 2) = AirMatrices.C6(i, 2) - C6ARG1 - C6ARG2_2;
-        AirMatrices.C6(i + 1, 1) = AirMatrices.C6(i + 1, 1) + C6ARG1 + C6ARG2_2;
-        AirMatrices.C6_a(i) = -C6ARG1 + C6ARG2_1;
-
-        C7ARG = (AirVariabes.Cag(i, 1) + AirVariabes.Cag(i, 2)) / 2;
-        AirMatrices.C7(i) = AirMatrices.C7(i) - C7ARG;
-        AirMatrices.C7(i + 1) = AirMatrices.C7(i + 1) + C7ARG;
-    end
+function AirMatrices = Air_MAT(AirVariabes, InitialValues)
+
+    AirMatrices.C1 = InitialValues.C1;
+    AirMatrices.C2 = InitialValues.C2;
+    AirMatrices.C3 = InitialValues.C3;
+    AirMatrices.C4 = InitialValues.C4;
+    AirMatrices.C5 = InitialValues.C5;
+    AirMatrices.C6 = InitialValues.C6;
+
+    ModelSettings = io.getModelSettings();
+
+    for i = 1:ModelSettings.NL
+        AirMatrices.C1(i, 1) = AirMatrices.C1(i, 1) + AirVariabes.Cah(i, 1) * ModelSettings.DeltZ(i) / 2;
+        AirMatrices.C1(i + 1, 1) = AirMatrices.C1(i + 1, 1) + AirVariabes.Cah(i, 2) * ModelSettings.DeltZ(i) / 2;
+
+        AirMatrices.C2(i, 1) = AirMatrices.C2(i, 1) + AirVariabes.CaT(i, 1) * ModelSettings.DeltZ(i) / 2;
+        AirMatrices.C2(i + 1, 1) = AirMatrices.C2(i + 1, 1) + AirVariabes.CaT(i, 2) * ModelSettings.DeltZ(i) / 2;
+
+        AirMatrices.C3(i, 1) = AirMatrices.C3(i, 1) + AirVariabes.Caa(i, 1) * ModelSettings.DeltZ(i) / 2;
+        AirMatrices.C3(i + 1, 1) = AirMatrices.C3(i + 1, 1) + AirVariabes.Caa(i, 2) * ModelSettings.DeltZ(i) / 2;
+
+        C4ARG1 = (AirVariabes.Kah(i, 1) + AirVariabes.Kah(i, 2)) / (2 * ModelSettings.DeltZ(i));
+        C4ARG2_1 = AirVariabes.Vah(i, 1) / 3 + AirVariabes.Vah(i, 2) / 6;
+        C4ARG2_2 = AirVariabes.Vah(i, 1) / 6 + AirVariabes.Vah(i, 2) / 3;
+        AirMatrices.C4(i, 1) = AirMatrices.C4(i, 1) + C4ARG1 - C4ARG2_1;
+        AirMatrices.C4(i, 2) = AirMatrices.C4(i, 2) - C4ARG1 - C4ARG2_2;
+        AirMatrices.C4(i + 1, 1) = AirMatrices.C4(i + 1, 1) + C4ARG1 + C4ARG2_2;
+        AirMatrices.C4_a(i) = -C4ARG1 + C4ARG2_1;
+
+        C5ARG1 = (AirVariabes.KaT(i, 1) + AirVariabes.KaT(i, 2)) / (2 * ModelSettings.DeltZ(i));
+        C5ARG2_1 = AirVariabes.VaT(i, 1) / 3 + AirVariabes.VaT(i, 2) / 6;
+        C5ARG2_2 = AirVariabes.VaT(i, 1) / 6 + AirVariabes.VaT(i, 2) / 3;
+        AirMatrices.C5(i, 1) = AirMatrices.C5(i, 1) + C5ARG1 - C5ARG2_1;
+        AirMatrices.C5(i, 2) = AirMatrices.C5(i, 2) - C5ARG1 - C5ARG2_2;
+        AirMatrices.C5(i + 1, 1) = AirMatrices.C5(i + 1, 1) + C5ARG1 + C5ARG2_2;
+        AirMatrices.C5_a(i) = -C5ARG1 + C5ARG2_1;
+
+        C6ARG1 = (AirVariabes.Kaa(i, 1) + AirVariabes.Kaa(i, 2)) / (2 * ModelSettings.DeltZ(i));
+        C6ARG2_1 = AirVariabes.Vaa(i, 1) / 3 + AirVariabes.Vaa(i, 2) / 6;
+        C6ARG2_2 = AirVariabes.Vaa(i, 1) / 6 + AirVariabes.Vaa(i, 2) / 3;
+        AirMatrices.C6(i, 1) = AirMatrices.C6(i, 1) + C6ARG1 - C6ARG2_1;
+        AirMatrices.C6(i, 2) = AirMatrices.C6(i, 2) - C6ARG1 - C6ARG2_2;
+        AirMatrices.C6(i + 1, 1) = AirMatrices.C6(i + 1, 1) + C6ARG1 + C6ARG2_2;
+        AirMatrices.C6_a(i) = -C6ARG1 + C6ARG2_1;
+
+        C7ARG = (AirVariabes.Cag(i, 1) + AirVariabes.Cag(i, 2)) / 2;
+        AirMatrices.C7(i) = AirMatrices.C7(i) - C7ARG;
+        AirMatrices.C7(i + 1) = AirMatrices.C7(i + 1) + C7ARG;
+    end
diff --git a/src/Air_sub.m b/src/+dryair/solveDryAirEquations.m
similarity index 100%
rename from src/Air_sub.m
rename to src/+dryair/solveDryAirEquations.m
diff --git a/src/Air_Solve.m b/src/+dryair/solveTridiagonalMatrixEquations.m
similarity index 97%
rename from src/Air_Solve.m
rename to src/+dryair/solveTridiagonalMatrixEquations.m
index cf1e6507..6239ac8b 100644
--- a/src/Air_Solve.m
+++ b/src/+dryair/solveTridiagonalMatrixEquations.m
@@ -1,17 +1,17 @@
-function [AirMatrices, P_gg, RHS] = Air_Solve(RHS, AirMatrices)
-
-    ModelSettings = io.getModelSettings();
-    RHS(1) = RHS(1) / AirMatrices.C6(1, 1);
-
-    for i = 2:ModelSettings.NN
-        AirMatrices.C6(i, 1) = AirMatrices.C6(i, 1) - AirMatrices.C6_a(i - 1) * AirMatrices.C6(i - 1, 2) / AirMatrices.C6(i - 1, 1);
-        RHS(i) = (RHS(i) - AirMatrices.C6_a(i - 1) * RHS(i - 1)) / AirMatrices.C6(i, 1);
-    end
-
-    for i = ModelSettings.NL:-1:1
-        RHS(i) = RHS(i) - AirMatrices.C6(i, 2) * RHS(i + 1) / AirMatrices.C6(i, 1);
-    end
-
-    for i = 1:ModelSettings.NN
-        P_gg(i) = RHS(i);
-    end
+function [AirMatrices, P_gg, RHS] = Air_Solve(RHS, AirMatrices)
+
+    ModelSettings = io.getModelSettings();
+    RHS(1) = RHS(1) / AirMatrices.C6(1, 1);
+
+    for i = 2:ModelSettings.NN
+        AirMatrices.C6(i, 1) = AirMatrices.C6(i, 1) - AirMatrices.C6_a(i - 1) * AirMatrices.C6(i - 1, 2) / AirMatrices.C6(i - 1, 1);
+        RHS(i) = (RHS(i) - AirMatrices.C6_a(i - 1) * RHS(i - 1)) / AirMatrices.C6(i, 1);
+    end
+
+    for i = ModelSettings.NL:-1:1
+        RHS(i) = RHS(i) - AirMatrices.C6(i, 2) * RHS(i + 1) / AirMatrices.C6(i, 1);
+    end
+
+    for i = 1:ModelSettings.NN
+        P_gg(i) = RHS(i);
+    end