Fix issue 181

src/+conductivity/+hydraulicConductivity/calculateDTheta.m

@@ -0,0 +1,43 @@
+function dtheta = calculateDTheta(subRoutine, heat_term, theta_s, theta_r, theta_m, gama_hh, theta_ll, theta_l, theta_uu, theta_u, hh, h, hh_frez, h_frez, phi_s, lamda, alpha, n, m)
+
+    % get soil constants
+    SoilConstants = io.getSoilConstants();
+
+    switch subRoutine
+        case 0
+            % for DTheta_LLh, heat_term = hh
+            % for DTheta_UUh, heat_term = hh + hh_frez
+            A = (-theta_r) / (abs(heat_term) * log(abs(SoilConstants.hd / SoilConstants.hm))) * (1 - (1 + abs(alpha * heat_term)^n)^(-m));
+            B = alpha * n * m * (theta_m - gama_hh * theta_r);
+            C = ((1 + abs(alpha * heat_term)^n)^(-m - 1));
+            D = (abs(alpha * heat_term)^(n - 1));
+            dtheta = A -  B * C * D;
+            if (hh + hh_frez) <= SoilConstants.hd
+                dtheta = 0;
+            end
+        case 1
+            % heat_term = hh or  hh + hh_frez
+            A = (theta_m - theta_r) * alpha * n;
+            B =  abs(alpha * heat_term)^(n - 1) * (-m);
+            C = (1 + abs(alpha * heat_term)^n)^(-m - 1);
+            dtheta = A * B * C;
+        case 2
+            A = theta_ll - theta_l;
+            B = hh + hh_frez - h - h_frez;
+            dtheta = A / B;
+        case 3
+            A = theta_uu - theta_u;
+            B = hh - h;
+            dtheta = A / B;
+        case 4
+            dtheta = theta_s / phi_s * (hh / phi_s)^(-1 * lamda - 1);
+        case 5
+            % this case differs with case 1 only regarding theta_m and theta_s
+            % this might be a bug, see issue 181
+            % heat_term = hh or  hh + hh_frez
+            A = (theta_s - theta_r) * alpha * n;
+            B =  abs(alpha * heat_term)^(n - 1) * (-m);
+            C = (1 + abs(alpha * heat_term)^n)^(-m - 1);
+            dtheta = A * B * C;
+    end
+end

src/+conductivity/+hydraulicConductivity/calculateDTheta_UUh.m

@@ -0,0 +1,65 @@
+function dtheta_uuh = calculateDTheta_UUh(theta_uu, theta_m, theta_ll, gama_hh, SoilVariables, VanGenuchten)
+    theta_s = VanGenuchten.Theta_s;
+    theta_r = VanGenuchten.Theta_r;
+    alpha = VanGenuchten.Alpha;
+    n = VanGenuchten.n;
+    m = VanGenuchten.m;
+
+    theta_u = SoilVariables.Theta_U;
+    theta_l = SoilVariables.Theta_L;
+    hh = SoilVariables.hh;
+    h = SoilVariables.h;
+    h_frez = SoilVariables.h_frez;
+    hh_frez = SoilVariables.hh_frez;
+    phi_s = SoilVariables.Phi_s;
+    lamda = SoilVariables.Lamda;
+
+    % get model settings
+    ModelSettings = io.getModelSettings();
+
+    heatTerm = hh + hh_frez;
+    if ModelSettings.SWCC == 1
+        if ModelSettings.SFCC == 1
+            if hh >= -1
+                if hh_frez >= 0
+                    dtheta_uuh = 0;
+                else
+                    if ModelSettings.Thmrlefc
+                        subRoutine = 0;
+                        dtheta_uuh = conductivity.hydraulicConductivity.calculateDTheta(subRoutine, heatTerm, theta_s, theta_r, theta_m, gama_hh, theta_ll, theta_l, theta_uu, theta_u, hh, h, hh_frez, h_frez, phi_s, lamda, alpha, n, m);
+                    elseif abs(hh - h) < 1e-3
+                        subRoutine = 1;
+                        dtheta_uuh = conductivity.hydraulicConductivity.calculateDTheta(subRoutine, heatTerm, theta_s, theta_r, theta_m, gama_hh, theta_ll, theta_l, theta_uu, theta_u, hh, h, hh_frez, h_frez, phi_s, lamda, alpha, n, m);
+                    else
+                        subRoutine = 2;
+                        dtheta_uuh = conductivity.hydraulicConductivity.calculateDTheta(subRoutine, heatTerm, theta_s, theta_r, theta_m, gama_hh, theta_ll, theta_l, theta_uu, theta_u, hh, h, hh_frez, h_frez, phi_s, lamda, alpha, n, m);
+                    end
+                end
+            elseif ModelSettings.Thmrlefc
+                subRoutine = 0;
+                dtheta_uuh = conductivity.hydraulicConductivity.calculateDTheta(subRoutine, heatTerm, theta_s, theta_r, theta_m, gama_hh, theta_ll, theta_l, theta_uu, theta_u, hh, h, hh_frez, h_frez, phi_s, lamda, alpha, n, m);
+            end
+        elseif hh >= -1e-6 || hh <= -1e7
+            dtheta_uuh = 0;
+        elseif ModelSettings.Thmrlefc || abs(hh - h) < 1e-3
+            subRoutine = 1;
+            dtheta_uuh = conductivity.hydraulicConductivity.calculateDTheta(subRoutine, heatTerm, theta_s, theta_r, theta_m, gama_hh, theta_ll, theta_l, theta_uu, theta_u, hh, h, hh_frez, h_frez, phi_s, lamda, alpha, n, m);
+        else
+            subRoutine = 3;
+            dtheta_uuh = conductivity.hydraulicConductivity.calculateDTheta(subRoutine, heatTerm, theta_s, theta_r, theta_m, gama_hh, theta_ll, theta_l, theta_uu, theta_u, hh, h, hh_frez, h_frez, phi_s, lamda, alpha, n, m);
+        end
+    else
+        if hh >= phi_s || hh <= -1e7
+            dtheta_uuh = 0;
+        elseif ModelSettings.Thmrlefc
+            subRoutine = 4;
+            dtheta_uuh = conductivity.hydraulicConductivity.calculateDTheta(subRoutine, heatTerm, theta_s, theta_r, theta_m, gama_hh, theta_ll, theta_l, theta_uu, theta_u, hh, h, hh_frez, h_frez, phi_s, lamda, alpha, n, m);
+        elseif abs(hh - h) < 1e-3
+            subRoutine = 5;
+            dtheta_uuh = conductivity.hydraulicConductivity.calculateDTheta(subRoutine, heatTerm, theta_s, theta_r, theta_m, gama_hh, theta_ll, theta_l, theta_uu, theta_u, hh, h, hh_frez, h_frez, phi_s, lamda, alpha, n, m);
+        else
+            subRoutine = 3;
+            dtheta_uuh = conductivity.hydraulicConductivity.calculateDTheta(subRoutine, heatTerm, theta_s, theta_r, theta_m, gama_hh, theta_ll, theta_l, theta_uu, theta_u, hh, h, hh_frez, h_frez, phi_s, lamda, alpha, n, m);
+        end
+    end
+end

src/+conductivity/+hydraulicConductivity/calculateGama_hh.m

@@ -0,0 +1,12 @@
+function gama_hh = calculateGama_hh(hh)
+    % get soil constants
+    SoilConstants = io.getSoilConstants();
+
+    if abs(hh) >= abs(SoilConstants.hd)
+        gama_hh = 0;
+    elseif abs(hh) >= abs(SoilConstants.hm)
+        gama_hh = log(abs(SoilConstants.hd) / abs(hh)) / log(abs(SoilConstants.hd) / abs(SoilConstants.hm));
+    else
+        gama_hh = 1;
+    end
+end

src/+conductivity/+hydraulicConductivity/calculateKL_h.m

@@ -0,0 +1,23 @@
+function kl_h = calculateKL_h(mu_w, se, Ks, m)
+
+    % load Constants
+    Constants = io.define_constants();
+
+    MU_WN = Constants.MU_W0 * exp(Constants.MU1 / (8.31441 * (20 + 133.3)));
+    CKT = MU_WN / mu_w;
+    if se == 0
+        kl_h = 0;
+    else
+        kl_h = CKT * Ks * (se^(0.5)) * (1 - (1 - se^(1 / m))^m)^2;
+    end
+    if kl_h <= 1E-20
+        kl_h = 1E-20;
+    end
+    if isnan(kl_h) == 1
+        kl_h = 0;
+        warning('\n case "isnan(kl_h) == 1", set "kl_h = 0" \r');
+    end
+    if ~isreal(kl_h)
+        warning('\n case "~isreal(kl_h)", dont know what to do! \r');
+    end
+end

src/+conductivity/+hydraulicConductivity/calculateSe.m

@@ -0,0 +1,57 @@
+function se = calculateSe(theta_ll, gama_hh, SoilVariables)
+    % get model settings
+    ModelSettings = io.getModelSettings();
+
+    % get soil constants
+    SoilConstants = io.getSoilConstants();
+
+    hh_frez = SoilVariables.hh_frez;
+    POR = SoilVariables.POR;
+    hh = SoilVariables.hh;
+    phi_s = SoilVariables.Phi_s;
+    se = SoilVariables.Se;
+    seCalculated = theta_ll / POR;
+
+    if ModelSettings.SWCC == 1
+        if ModelSettings.SFCC == 1
+            if hh >= -1
+                if hh_frez >= 0
+                    se = 1;
+                elseif Thmrlefc
+                    if (hh + hh_frez) <= SoilConstants.hd
+                        se = 0;
+                    else
+                        se = seCalculated;
+                    end
+                end
+            elseif ModelSettings.Thmrlefc
+                if (hh + hh_frez) <= SoilConstants.hd
+                    se = 0;
+                else
+                    se = seCalculated;
+                end
+            end
+        elseif theta_ll <= 0.06 || hh <= -1e7
+            se = 0;
+        else
+            se = seCalculated;
+        end
+    else
+        if hh + hh_frez <= -1e7 || hh <= -1e7
+            se = 0;
+        elseif hh + hh_frez >= phi_s
+            se = 1;
+        else
+            se = seCalculated;
+        end
+    end
+
+    if se > 1
+        se = 1;
+    elseif se < 0
+        se = 0;
+    end
+    if isnan(se) == 1
+        warning('\n case "isnan(se) == 1" happens. Dont know what to do! \r');
+    end
+end

src/+conductivity/+hydraulicConductivity/calculateTheta.m

@@ -0,0 +1,25 @@
+function theta = calculateTheta(subRoutine, theta_m, heat_term, gama_hh, theta_s, theta_r, lamda, phi_s, alpha, n, m)
+    % get soil constants
+    SoilConstants = io.getSoilConstants();
+    hd = SoilConstants.hd;
+
+    switch subRoutine
+        case 0
+            % heat_term = hh or hh + hh_frez
+            theta = gama_hh * theta_r + (theta_m - gama_hh * theta_r) / (1 + abs(alpha * heat_term)^n)^m;
+            if heat_term <= hd
+                theta = 0;
+            end
+        case 1
+            % heat_term = hh or hh + hh_frez
+            theta = theta_s * (heat_term / phi_s)^(-1 * lamda);
+            if heat_term <= -1e7
+                theta = Theta_r;
+            elseif heat_term >= phi_s
+                theta = theta_s;
+            end
+        case 2
+            % heat_term = hh or hh + hh_frez
+            theta = Theta_r + (theta_s - Theta_r) / (1 + abs(alpha * heat_term)^n)^m;
+    end
+end

src/+conductivity/+hydraulicConductivity/calculateTheta_II.m

@@ -0,0 +1,23 @@
+function theta_ii = calculateTheta_II(tt, xcap, hh, Theta_II)
+
+    % get model settings
+    ModelSettings = io.getModelSettings();
+
+    Tf1 = 273.15 + 1;
+    Tf2 = 273.15 - 3;
+
+    theta_ii = Theta_II;
+
+    if hh <= -1e7
+        theta_ii = 0;
+    end
+    if ModelSettings.SWCC == 1 && ModelSettings.SFCC ~= 1
+        if tt + 273.15 > Tf1
+            theta_ii = 0;
+        elseif tt + 273.15 >= Tf2
+            theta_ii = 0.5 * (1 - sin(pi() * (tt + 273.15 - 0.5 * Tf1 - 0.5 * Tf2) / (Tf1 - Tf2))) * xcap;
+        else
+            theta_ii = xcap;
+        end
+    end
+end

src/+conductivity/+hydraulicConductivity/calculateTheta_LL.m

@@ -0,0 +1,54 @@
+function theta_ll = calculateTheta_LL(theta_uu, theta_ii, theta_m, gama_hh, SoilVariables, VanGenuchten)
+    % get model settings
+    ModelSettings = io.getModelSettings();
+
+    % load Constants
+    Constants = io.define_constants();
+
+    hh = SoilVariables.hh;
+    hh_frez = SoilVariables.hh_frez;
+    phi_s = SoilVariables.Phi_s;
+    lamda = SoilVariables.Lamda;
+
+    theta_s = VanGenuchten.Theta_s;
+    theta_r = VanGenuchten.Theta_r;
+    alpha = VanGenuchten.Alpha;
+    n = VanGenuchten.n;
+    m = VanGenuchten.m;
+
+    heatTerm = hh + hh_frez;
+
+    % calculate theta_ll
+    theta_ll = SoilVariables.Theta_LL;
+    if ModelSettings.SWCC == 1
+        if ModelSettings.SFCC == 1
+            if hh >= -1 && hh_frez >= 0
+                theta_ll = theta_s;
+            elseif ModelSettings.Thmrlefc
+                subRoutine = 0;
+                theta_ll = conductivity.hydraulicConductivity.calculateTheta(subRoutine, theta_m, heatTerm, gama_hh, theta_s, theta_r, lamda, phi_s, alpha, n, m);
+            end
+        else
+            if hh >= -1e-6 && TT + 273.15 > (273.15 + 1)
+                theta_ll = theta_s;
+            elseif hh <= -1e7
+                theta_ll = theta_r;
+            elseif TT + 273.15 > (273.15 + 1)
+                subRoutine = 2;
+                theta_ll = conductivity.hydraulicConductivity.calculateTheta(subRoutine, theta_m, heatTerm, gama_hh, theta_s, theta_r, lamda, phi_s, alpha, n, m);
+            else
+                theta_ll = theta_uu - theta_ii * Constants.RHOI / Constants.RHOL;
+            end
+            if theta_ll <= 0.06
+                theta_ll = 0.06;
+            end
+        end
+    else
+        if hh >= phi_s
+            subRoutine = 1;
+            theta_ll = conductivity.hydraulicConductivity.calculateTheta(subRoutine, theta_m, heatTerm, gama_hh, theta_s, theta_r, lamda, phi_s, alpha, n, m);
+        elseif hh <= -1e7
+            theta_ll = theta_r;
+        end
+    end
+end

src/+conductivity/+hydraulicConductivity/calculateTheta_UU.m

@@ -0,0 +1,43 @@
+function theta_uu = calculateTheta_UU(theta_m, gama_hh, SoilVariables, VanGenuchten)
+
+    hh = SoilVariables.hh;
+    phi_s = SoilVariables.Phi_s;
+    lamda = SoilVariables.Lamda;
+
+    theta_s = VanGenuchten.Theta_s;
+    theta_r = VanGenuchten.Theta_r;
+    alpha = VanGenuchten.Alpha;
+    n = VanGenuchten.n;
+    m = VanGenuchten.m;
+
+    % get model settings
+    ModelSettings = io.getModelSettings();
+
+    % calculate theta_uu
+    if ModelSettings.SWCC == 1
+        if ModelSettings.SFCC == 1
+            if hh >= -1
+                theta_uu = theta_s;
+            elseif ModelSettings.Thmrlefc
+                if gama_hh == 0
+                    theta_uu = 0;
+                else
+                    subRoutine = 0;
+                    theta_uu = conductivity.hydraulicConductivity.calculateTheta(subRoutine, theta_m, hh, gama_hh, theta_s, theta_r, lamda, phi_s, alpha, n, m);
+                end
+            end
+        else
+            if hh >= -1e-6
+                theta_uu = theta_s;
+            elseif hh <= -1e7
+                theta_uu = theta_r;
+            else
+                subRoutine = 2;
+                theta_uu = conductivity.hydraulicConductivity.calculateTheta(subRoutine, theta_m, hh, gama_hh, theta_s, theta_r, lamda, phi_s, alpha, n, m);
+            end
+        end
+    else
+        subRoutine = 1;
+        theta_uu = conductivity.hydraulicConductivity.calculateTheta(subRoutine, theta_m, hh, gama_hh, theta_s, theta_r, lamda, phi_s, alpha, n, m);
+    end
+end

src/+conductivity/+hydraulicConductivity/calculateTheta_m.m

@@ -0,0 +1,15 @@
+function theta_m = calculateTheta_m(gama_hh, VanGenuchten, POR)
+
+    theta_s = VanGenuchten.Theta_s;
+    theta_r = VanGenuchten.Theta_r;
+    alpha = VanGenuchten.Alpha;
+    n = VanGenuchten.n;
+    m = VanGenuchten.m;
+
+    theta_m = gama_hh * theta_r + (theta_s - gama_hh * theta_r) * (1 + abs(alpha * (-1))^n)^m;
+    if theta_m >= POR
+        theta_m = POR;
+    elseif theta_m <= theta_s
+        theta_m = theta_s;
+    end
+end