diff --git a/src/+dryair/assembleCoefficientMatrices.m b/src/+dryair/assembleCoefficientMatrices.m
index 0f55e318..4fd49003 100644
--- a/src/+dryair/assembleCoefficientMatrices.m
+++ b/src/+dryair/assembleCoefficientMatrices.m
@@ -1,4 +1,8 @@
-function [RHS, AirMatrices, SAVE] = Air_EQ(AirMatrices, Delt_t, P_g)
+function [RHS, AirMatrices, SAVE] = assembleCoefficientMatrices(AirMatrices, Delt_t, P_g)
+    %{
+        Assemble the coefficient matrices of Equation 4.32 STEMMUS Technical
+        Notes, page 44, for dry air equation.
+    %}
 
     C1 = AirMatrices.C1
     C2 = AirMatrices.C2
diff --git a/src/+dryair/calculateBoundaryConditions.m b/src/+dryair/calculateBoundaryConditions.m
index befdf99f..19adce0f 100644
--- a/src/+dryair/calculateBoundaryConditions.m
+++ b/src/+dryair/calculateBoundaryConditions.m
@@ -1,8 +1,12 @@
-function [RHS, AirMatrices] = Air_BC(BoundaryCondition, AirMatrices, ForcingData, RHS, KT)
+function [RHS, AirMatrices] = calculateBoundaryConditions(BoundaryCondition, AirMatrices, ForcingData, RHS, KT)
+    %{
+        Determine the boundary condition for solving the dry air equation.
+    %}
 
     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;
diff --git a/src/+dryair/calculateDryAirParameters.m b/src/+dryair/calculateDryAirParameters.m
index ed0e5c8c..03bed157 100644
--- a/src/+dryair/calculateDryAirParameters.m
+++ b/src/+dryair/calculateDryAirParameters.m
@@ -1,6 +1,9 @@
-function AirVariabes = AirPARM(SoilVariables, GasDispersivity, TransportCoefficient, InitialValues, GasDispersivity,...
-                               P_gg, Xah, XaT, Xaa, RHODA)
-
+function AirVariabes = calculateDryAirParameters(SoilVariables, GasDispersivity, TransportCoefficient, InitialValues, GasDispersivity,...
+                                                 P_gg, Xah, XaT, Xaa, RHODA)
+    %{
+        Calculate all the parameters related to dry air equation e.g., Equation
+        3.59-3.64, STEMMUS Technical Notes, page 27-28.
+    %}
     ModelSettings = io.getModelSettings();
     Constants = io.define_constants();
 
diff --git a/src/+dryair/calculateMatricCoefficients.m b/src/+dryair/calculateMatricCoefficients.m
index 23f19003..403ceca7 100644
--- a/src/+dryair/calculateMatricCoefficients.m
+++ b/src/+dryair/calculateMatricCoefficients.m
@@ -1,4 +1,9 @@
-function AirMatrices = Air_MAT(AirVariabes, InitialValues)
+function AirMatrices = calculateMatricCoefficients(AirVariabes, InitialValues)
+    %{
+        Calculate all the parameters related to matric coefficients e.g.,
+        c1-c7 as in Equation 4.32 STEMMUS Technical Notes, page 44, which is
+        an example for soil moisture equation, but for dry air equation.
+    %}
 
     AirMatrices.C1 = InitialValues.C1;
     AirMatrices.C2 = InitialValues.C2;
diff --git a/src/+dryair/solveDryAirEquations.m b/src/+dryair/solveDryAirEquations.m
index 3870ba19..bbe94e80 100644
--- a/src/+dryair/solveDryAirEquations.m
+++ b/src/+dryair/solveDryAirEquations.m
@@ -1,13 +1,18 @@
-function [RHS, SAVE, P_gg] = Air_sub(SoilVariables, GasDispersivity, TransportCoefficient, InitialValues, GasDispersivity,...
-                                     BoundaryCondition, P_gg, Xah, XaT, Xaa, RHODA, KT, Delt_t)
+function [RHS, SAVE, P_gg] = solveDryAirEquations(SoilVariables, GasDispersivity, TransportCoefficient, InitialValues, GasDispersivity,...
+                                                  BoundaryCondition, P_gg, Xah, XaT, Xaa, RHODA, KT, Delt_t)
+    %{
+        Solve the dry air equation with the Thomas algorithm to update the soil
+        air pressure 'P_gg', the finite difference time-stepping scheme is
+        exampled as for the soil moisture equation, which derived in 'STEMMUS
+        Technical Notes' section 4, Equation 4.32.
+    %}
+    AirVariabes = dryair.calculateDryAirParameters(SoilVariables, GasDispersivity, TransportCoefficient, InitialValues, GasDispersivity,...
+                                                   P_gg, Xah, XaT, Xaa, RHODA);
 
-    AirVariabes = AirPARM(SoilVariables, GasDispersivity, TransportCoefficient, InitialValues, GasDispersivity,...
-                          P_gg, Xah, XaT, Xaa, RHODA);
+    AirMatrices = dryair.calculateMatricCoefficients(AirVariabes, InitialValues);
 
-    AirMatrices = Air_MAT(AirVariabes, InitialValues);
+    [RHS, AirMatrices, SAVE] = dryair.assembleCoefficientMatrices(AirMatrices, Delt_t, P_g);
 
-    [RHS, AirMatrices, SAVE] = Air_EQ(AirMatrices, Delt_t, P_g);
+    [RHS, AirMatrices] = dryair.calculateBoundaryConditions(BoundaryCondition, AirMatrices, ForcingData, RHS, KT);
 
-    [RHS, AirMatrices] = Air_BC(BoundaryCondition, AirMatrices, ForcingData, RHS, KT);
-
-    [AirMatrices, P_gg, RHS] = Air_Solve(RHS, AirMatrices);
+    [AirMatrices, P_gg, RHS] = dryair.solveTridiagonalMatrixEquations(RHS, AirMatrices);
diff --git a/src/+dryair/solveTridiagonalMatrixEquations.m b/src/+dryair/solveTridiagonalMatrixEquations.m
index 6239ac8b..4afe27f9 100644
--- a/src/+dryair/solveTridiagonalMatrixEquations.m
+++ b/src/+dryair/solveTridiagonalMatrixEquations.m
@@ -1,4 +1,8 @@
-function [AirMatrices, P_gg, RHS] = Air_Solve(RHS, AirMatrices)
+function [AirMatrices, P_gg, RHS] = solveTridiagonalMatrixEquations(RHS, AirMatrices)
+    %{
+        Use Thomas algorithm to solve the tridiagonal matrix equations, which is
+        in the form of Equation 4.25 STEMMUS Technical Notes, page 41.
+    %}
 
     ModelSettings = io.getModelSettings();
     RHS(1) = RHS(1) / AirMatrices.C6(1, 1);
diff --git a/src/STEMMUS_SCOPE.m b/src/STEMMUS_SCOPE.m
index 3ed21a07..38c4de0b 100644
--- a/src/STEMMUS_SCOPE.m
+++ b/src/STEMMUS_SCOPE.m
@@ -804,7 +804,8 @@
         end
 
         if Soilairefc == 1
-            run Air_sub;
+            [RHS, SAVE, P_gg] = dryair.solveDryAirEquations(SoilVariables, GasDispersivity, TransportCoefficient, InitialValues, GasDispersivity,...
+                                                            BoundaryCondition, P_gg, Xah, XaT, Xaa, RHODA, KT, Delt_t);
         end
 
         if Thmrlefc == 1