New Glacier Module

src/CMakeLists.txt

@@ -62,12 +62,21 @@ set(MODULE_SRCS
         modules/PBSM3D.cpp
 		modules/snobal.cpp
 		modules/fsm.cpp
+		modules/katabatic_routing_glacier.cpp
 
 		modules/testing/mpi_example.cpp
 
 
 		CACHE INTERNAL "" FORCE )
 
+set(SUBMODULE_SRCS
+	# modules/submodules/water_flux.cpp
+        modules/submodules/Glacier.cpp
+		modules/submodules/katabatic_melt_energy.cpp
+		modules/submodules/melt_routing_glacier.cpp
+
+
+        CACHE INTERNAL "" FORCE )
 
 set(CHM_SRCS
 		#main.cpp needs to be added below so we can re use CHM_SRCS in the gtest build
@@ -77,6 +86,7 @@ set(CHM_SRCS
 		metdata.cpp
 
 		physics/Atmosphere.cpp
+        physics/PhysConst.cpp
 
 		mesh/triangulation.cpp
         mesh/ugrid_writer.cpp
@@ -110,6 +120,7 @@ set (HEADER_FILES
         modules/interp_met
 		modules/snobal
 		modules/testing
+        modules/submodules
 		math
         libmaw
         interpolation
@@ -234,6 +245,7 @@ add_executable(
 		${CHM_SRCS}
 		${FILTER_SRCS}
 		${MODULE_SRCS}
+        ${SUBMODULE_SRCS}
 )
 
 
@@ -412,17 +424,23 @@ if (BUILD_TESTS)
 	message(STATUS "Tests enabled. Run with make check")
 
 	set(TEST_SRCS
-			tests/test_station.cpp
-			tests/test_interpolation.cpp
-			tests/test_timeseries.cpp
-			tests/test_core.cpp
-			tests/test_variablestorage.cpp
-			tests/test_metdata.cpp
-			tests/test_netcdf.cpp
+        #tests/test_station.cpp
+        #	tests/test_interpolation.cpp
+        #		tests/test_timeseries.cpp
+        #	tests/test_core.cpp
+        #	tests/test_variablestorage.cpp
+        #	tests/test_metdata.cpp
+        #	tests/test_netcdf.cpp
 			#    test_mesh.cpp
 			tests/test_regexptokenizer.cpp
 			#    test_daily.cpp
-            tests/test_triangulation.cpp
+            #     tests/test_triangulation.cpp
+            tests/submoduletests/test_data_base.cpp
+			#tests/submoduletests/test_water_flux.cpp
+            tests/submoduletests/test_glacier.cpp
+			tests/submoduletests/test_katabatic_melt_energy.cpp
+			tests/submoduletests/test_melt_routing.cpp
+            tests/submoduletests/test_bisection.cpp
 			tests/main.cpp
     )
 
@@ -435,6 +453,7 @@ if (BUILD_TESTS)
 			${CHM_SRCS}
 			${FILTER_SRCS}
 			${MODULE_SRCS}
+            ${SUBMODULE_SRCS}
 			${TEST_SRCS}
 	)
 	set_target_properties(runUnitTests

src/core.cpp

@@ -1155,10 +1155,17 @@ void core::config_output(pt::ptree &value)
                 OGRCoordinateTransformation::DestroyCT(coordTrans);
             }
 
-            if(!_mesh->is_geographic())
-                out.face = _mesh->locate_face(out.x, out.y);
-            else
-                out.face = _mesh->locate_face(out.longitude, out.latitude);
+
+
+           // if(!_mesh->is_geographic())
+           // {
+           //     out.face = _mesh->locate_face(out.x, out.y);
+           // }
+           // else
+           // {
+           //     out.face = _mesh->locate_face(out.longitude, out.latitude);
+           // }
+           out.face = _mesh->face(0);
 
             if(out.face != nullptr)
             {

src/math/Bisection.hpp

@@ -0,0 +1,51 @@
+#pragma once
+
+#include <functional>
+#include <cmath>
+#include <limits>
+namespace math
+{
+    enum class Root { Found, DidNotConverge, NoRootGuaranteed };
+    struct Result
+    {
+        Root root;
+        double value;
+        size_t iterations;
+    };
+
+    inline Result bisection(
+        std::function<double(double)> f,
+        double low_guess,
+        double  high_guess,
+        double tol = 1e-10,
+        size_t max_iter = 1000
+    ) {
+        double fa = f(low_guess), fb = f(high_guess);
+
+        if (fa * fb > 0) 
+            return Result {
+                Root::NoRootGuaranteed,
+                std::numeric_limits<double>::quiet_NaN(),
+                0};
+
+        for (size_t i = 0; i < max_iter; ++i) {
+            double c = low_guess + (high_guess - low_guess) / 2;
+            double fc = f(c);
+
+            if (std::abs(fc) < tol || (high_guess - low_guess) / 2 < tol) 
+                return Result{Root::Found,c,i} ;
+
+            if (fa * fc < 0) {
+                high_guess = c;
+                fb = fc;
+            } else {
+                low_guess = c;
+                fa = fc;
+            }
+        }
+
+        return Result{Root::DidNotConverge,
+            std::numeric_limits<double>::quiet_NaN(),
+            max_iter};  // Did not converge
+    };
+};

src/mesh/triangulation.hpp

@@ -436,9 +436,10 @@ class face : public Fb
 
 //    void set_module_data(const std::string &module, face_info *fi);
 
-    template<typename T>
-    T& make_module_data(const std::string &module);
-
+    // Overload for module data constructor
+    template<typename T, typename... Args>
+    T& make_module_data(const std::string &module, Args... args);
+
     std::string _debug_name; //for debugging to find the elem that we want
     int _debug_ID; //also for debugging. ID == the position in the output order, starting at 0
     size_t cell_global_id;
@@ -1788,22 +1789,22 @@ timeseries::iterator face<Gt, Fb>::now()
 }
 
 
+// Overloaded for construtor that requires arguments
 template < class Gt, class Vb>
-template<typename T>
-T& face<Gt, Vb>::make_module_data(const std::string &module)
+template<typename T, typename... Args>
+T& face<Gt, Vb>::make_module_data(const std::string &module, Args... args)
 {
 
     //we don't already have this, make a new one.
     if(!_module_face_data[module])
     {
 //        T* data = new T;
-        _module_face_data[module] = std::make_unique<T>();
+        _module_face_data[module] = std::make_unique<T>(args...);
     }
 
     return get_module_data<T&>(module);
 }
 
-
 template < class Gt, class Fb>
 template < typename T>
 T& face<Gt, Fb>::get_module_data(const std::string &module)

src/modules/PBSM3D.cpp

@@ -423,7 +423,7 @@ void PBSM3D::run(mesh& domain)
         auto tol = [](double a, double b) -> bool { return fabs(a - b) < 1e-8; };
 
         // ice density
-        double rho_p = PhysConst::rho_ice;
+        double rho_p = PhysConst::rho_ice();
 #pragma omp for
         for (size_t i = 0; i < domain->size_local_faces(); i++)
         {
@@ -700,7 +700,7 @@ void PBSM3D::run(mesh& domain)
                         // c_4 = 0.5;
                         // g   = 9.81;
 
-                        return u2 * PhysConst::kappa / log(2.0 / (0.6131702345e-2 * ustar * ustar + .5 * lambda)) -
+                        return u2 * PhysConst::kappa() / log(2.0 / (0.6131702345e-2 * ustar * ustar + .5 * lambda)) -
                             ustar;
                     };
                     try
@@ -717,7 +717,7 @@ void PBSM3D::run(mesh& domain)
                 else
                 {
                     // follow PBSM (Pom & Li 2000; Alpine3D) and don't calculate the feedback of z0 on u*
-                    ustar = u2 * PhysConst::kappa / log(2.0 / 0.0002);
+                    ustar = u2 * PhysConst::kappa() / log(2.0 / 0.0002);
                 }
 
                 if (ustar >= u_star_saltation_threshold)
@@ -746,7 +746,7 @@ void PBSM3D::run(mesh& domain)
             {
                 // we still need a u* for spatial K estimation later
                 d.z0 = Snow::Z0_SNOW;
-                ustar = std::max(0.01, PhysConst::kappa * uref / log(Atmosphere::Z_U_R / d.z0));
+                ustar = std::max(0.01, PhysConst::kappa() * uref / log(Atmosphere::Z_U_R / d.z0));
             }
 
             // sanity checks
@@ -1153,7 +1153,7 @@ void PBSM3D::run(mesh& domain)
                     }
                 }
                 // Li and Pomeroy 2000
-                double l = PhysConst::kappa * (cz + d.z0) * l__max / (PhysConst::kappa * (cz + d.z0) + l__max);
+                double l = PhysConst::kappa() * (cz + d.z0) * l__max / (PhysConst::kappa() * (cz + d.z0) + l__max);
                 if (debug_output)
                     (*face)["l"_s] = l;
 

src/modules/Simple_Canopy.cpp

@@ -172,20 +172,20 @@ void Simple_Canopy::run(mesh_elem &face)
     // Canopy temperature is first approximated by the air temperature.
     double T1 = ta + mio::Cst::t_water_freezing_pt; // Canopy temperature (C to K)
 
-    double rho = air_pressure*1000/(PhysConst::Rgas*T1); // density of Air (pressure kPa to Pa = *1000)
+    double rho = air_pressure*1000/(PhysConst::RgasDry()*T1); // density of Air (pressure kPa to Pa = *1000)
 
     double U1 = U_R; // Wind speed (m/s) at height Z_vw [m] (top of canopy)
 
     // Aerodynamic resistance of canopy
-    ra = (log(Zref/Z0snow)*log(Zwind/Z0snow))/pow(PhysConst::kappa,2)/U1; // (s/m)
+    ra = (log(Zref/Z0snow)*log(Zwind/Z0snow))/pow(PhysConst::kappa(),2)/U1; // (s/m)
 
-    double deltaX = 0.622*PhysConst::Ls*Qs(air_pressure, T1)/(PhysConst::Rgas*(pow(T1,2))); // Must be (kg K-1)
+    double deltaX = 0.622*PhysConst::Ls()*Qs(air_pressure, T1)/(PhysConst::RgasDry()*(pow(T1,2))); // Must be (kg K-1)
 
     double q = (rh/100)*Qs(air_pressure, T1); // specific humidity (kg/kg)
 
     // snow surface temperature of snow in canopy
-    Ts = T1 + (Snow::emiss*(ilwr - PhysConst::sbc*pow(T1, 4.0)) + PhysConst::Ls*(q - Qs(air_pressure, T1))*rho/ra)/
-              (4.0*Snow::emiss*PhysConst::sbc*pow(T1, 3.0) + (PhysConst::Cp + PhysConst::Ls*deltaX)*rho/ra);
+    Ts = T1 + (Snow::emiss*(ilwr - PhysConst::sbc()*pow(T1, 4.0)) + PhysConst::Ls()*(q - Qs(air_pressure, T1))*rho/ra)/
+              (4.0*Snow::emiss*PhysConst::sbc()*pow(T1, 3.0) + (PhysConst::Cp() + PhysConst::Ls()*deltaX)*rho/ra);
 
     Ts -= mio::Cst::t_water_freezing_pt; // K to C
 
@@ -227,7 +227,7 @@ void Simple_Canopy::run(mesh_elem &face)
             Kstar_H = iswr * (1.0 - Alpha_c - Tauc * (1.0 - Albedo)); //  what is Kstar_H???
 
             // Incident long-wave at surface, "(W/m^2)"
-            Qlisn = ilwr * Vf_ + (1.0 - Vf_) * Vegetation::emiss_c * PhysConst::sbc * pow(T1, 4.0) + B_canopy * Kstar_H;
+            Qlisn = ilwr * Vf_ + (1.0 - Vf_) * Vegetation::emiss_c * PhysConst::sbc() * pow(T1, 4.0) + B_canopy * Kstar_H;
 
             // Incident short-wave at surface, "(W/m^2)"
             Qsisn = iswr * Tauc;
@@ -289,7 +289,7 @@ void Simple_Canopy::run(mesh_elem &face)
         Kd = iswr * (1.0 - Alpha_c - Tau_b_gap * (1.0 - Albedo));
 
         Qlisn = Vgap * ilwr + (1.0 - Vgap) * ((ilwr * Tau_b_gap +
-                                               (1.0 - Tau_b_gap) * Vegetation::emiss_c * PhysConst::sbc *
+                                               (1.0 - Tau_b_gap) * Vegetation::emiss_c * PhysConst::sbc() *
                                                pow(T1, 4.0f)) + B_canopy * Kd);
 
         Qsisn = cosxs * Qdfo * Tau_b_gap + Vgap * (iswr - Qdfo) + (1.0 - Vgap) * Tau_d * (iswr - Qdfo);
@@ -373,27 +373,27 @@ void Simple_Canopy::run(mesh_elem &face)
 
                 double Es = 611.15 * exp(22.452 * ta / (ta + 273.0));  // {sat pressure}
 
-                double SvDens = Es * PhysConst::M / (PhysConst::R * (ta + 273.0)); // {sat density}
+                double SvDens = Es * PhysConst::M() / (PhysConst::R() * (ta + 273.0)); // {sat density}
 
                 Lamb = 6.3e-4 * (ta + 273.0) + 0.0673;  // thermal conductivity of atmosphere
                 Nr = 2.0 * Snow::Radius * uVent / Atmosphere::KinVisc;  // Reynolds number
                 Nu = 1.79 + 0.606 * sqrt(Nr); // Nusselt number
                 SStar = M_PI * pow(Snow::Radius, 2) * (1.0 - Snow::AlbedoIce) *
                         iswr;  // SW to snow particle !!!! changed
                 A1 = Lamb * (ta + 273) * Nu;
-                B1 = PhysConst::Ls * PhysConst::M / (PhysConst::R * (ta + 273.0)) - 1.0;
+                B1 = PhysConst::Ls() * PhysConst::M() / (PhysConst::R() * (ta + 273.0)) - 1.0;
                 J = B1 / A1;
                 Sigma2 = rh / 100 - 1;
                 D = 2.06e-5 * pow((ta + 273.0) / 273.0, -1.75); // diffusivity of water vapour
                 C1 = 1.0 / (D * SvDens * Nu);
 
                 Alpha = 5.0;
-                Mpm = 4.0 / 3.0 * M_PI * PhysConst::rho_ice * pow(Snow::Radius, 3) *
+                Mpm = 4.0 / 3.0 * M_PI * PhysConst::rho_ice() * pow(Snow::Radius, 3) *
                       (1.0 + 3.0 / Alpha + 2.0 / pow(Alpha, 2));
 
                 // sublimation rate of single 'ideal' ice sphere:
 
-                double Vs = (2.0 * M_PI * Snow::Radius * Sigma2 - SStar * J) / (PhysConst::Ls * J + C1) / Mpm;
+                double Vs = (2.0 * M_PI * Snow::Radius * Sigma2 - SStar * J) / (PhysConst::Ls() * J + C1) / Mpm;
 
                 // snow exposure coefficient (Ce):
 
@@ -409,7 +409,7 @@ void Simple_Canopy::run(mesh_elem &face)
 
                 // calculate 'ice-bulb' temperature of intercepted snow:
 
-                double IceBulbT = ta - (Vi * PhysConst::Ls / 1e6 / PhysConst::Ci);
+                double IceBulbT = ta - (Vi * PhysConst::Ls() / 1e6 / PhysConst::Ci());
 
                 // determine whether canopy snow is unloaded:
 
@@ -431,8 +431,8 @@ void Simple_Canopy::run(mesh_elem &face)
                 // limit sublimation to canopy snow available and take sublimated snow away from canopy snow at timestep start
 
                 //Subl_Cpy = -data.Snow_load*Vi*Hs*Global::Interval*24*3600/Hs; // make W/m2 (original in CRHM)
-                Subl_Cpy = -data.Snow_load * Vi * PhysConst::Ls * global_param->dt() /
-                           PhysConst::Ls; // make W/m2 TODO: check Interval is same as dt() (in seconds
+                Subl_Cpy = -data.Snow_load * Vi * PhysConst::Ls() * global_param->dt() /
+                           PhysConst::Ls(); // make W/m2 TODO: check Interval is same as dt() (in seconds
                 // TODO: Hs/HS = 1 !!! (in CRHM, kept here for conistency...)
 
                 if (Subl_Cpy > data.Snow_load) {
@@ -719,4 +719,4 @@ void Simple_Canopy::load_checkpoint(mesh& domain, netcdf& chkpt)
 
     }
 
-}
+}