diff --git a/.github/workflows/Build_and_Run_Standalone_Test.yml b/.github/workflows/Build_and_Run_Standalone_Test.yml index 1efecc60..146d4abf 100644 --- a/.github/workflows/Build_and_Run_Standalone_Test.yml +++ b/.github/workflows/Build_and_Run_Standalone_Test.yml @@ -57,15 +57,6 @@ jobs: make cd .. ./run_cfe.sh FORCINGPET - - - name : Build and Run AETROOTZONE option - run: | - git clone https://github.com/NOAA-OWP/SoilMoistureProfiles extern/SoilMoistureProfiles - cd build - cmake ../ -DAETROOTZONE=ON - make - cd .. - ./run_cfe.sh AETROOTZONE - name: Build and Run BMI Unit Test run: | diff --git a/.github/workflows/ngen_integration.yaml b/.github/workflows/ngen_integration.yaml index 59459457..b9942afd 100644 --- a/.github/workflows/ngen_integration.yaml +++ b/.github/workflows/ngen_integration.yaml @@ -36,7 +36,7 @@ jobs: - name: Build the CFE Library run: | - cmake -B cmake_build -S . -DNGEN=ON -DCMAKE_C_FLAGS='-g -Og -fsanitize=address -Werror' + cmake -B cmake_build -S . -DNGEN=ON make -C cmake_build - name: Save CFE to a Temp Directory @@ -88,7 +88,7 @@ jobs: - name: Run NGen Test for CFE (cfe1.0) Coupled with PET run: | mv ${{ steps.ngen_id1.outputs.build-dir }} ./ngen-build/ - inputfile='extern/cfe/cfe/realizations/realization_cfe_pet_surfgiuh.json' + inputfile='extern/cfe/cfe/realizations/realization_cfe_pet.json' ./ngen-build/ngen ./data/catchment_data.geojson "cat-27" ./data/nexus_data.geojson "nex-26" $inputfile - name: Run Ngen Test for CFE (cfe2.0) Coupled with PET diff --git a/CMakeLists.txt b/CMakeLists.txt index 0d050207..4ed24c8b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -7,6 +7,9 @@ set(Red "${Esc}[32m") set(CMAKE_CXX_STANDARD 14) set(CMAKE_CXX_STANDARD_REQUIRED ON) +set(CMAKE_C_COMPILER $ENV{CC}) +set(CMAKE_CXX_COMPILER $ENV{CXX}) + # module setup options option(BASE "BASE" OFF) @@ -43,6 +46,7 @@ endif() # set the project name project(cfebmi VERSION 1.0.0 DESCRIPTION "OWP CFE BMI Module Shared Library") +set(CMAKE_BUILD_TYPE Debug) IF(CMAKE_BUILD_TYPE MATCHES Debug) message("Debug build.") ENDIF(CMAKE_BUILD_TYPE MATCHES Debug) @@ -51,39 +55,49 @@ message(CMAKE_CXX_COMPILER " ${CMAKE_CXX_COMPILER}") message(CMAKE_C_COMPILER " ${CMAKE_C_COMPILER}") message("CMAKE_BUILD_TYPE = ${CMAKE_BUILD_TYPE}") +set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -O0") +set(CMAKE_C_FLAGS_DEBUG "${CMAKE_C_FLAGS_DEBUG} -O0") + # add the executable ## cfe + aorc + pet + smp if(AETROOTZONE) -add_executable(${exe_name} ./src/main_cfe_aorc_pet_rz_aet.cxx ./src/cfe.c ./src/bmi_cfe.c ./src/giuh.c - ./src/conceptual_reservoir.c ./src/nash_cascade.c ./extern/aorc_bmi/src/aorc.c - ./extern/aorc_bmi/src/bmi_aorc.c ./extern/evapotranspiration/src/pet.c - ./extern/evapotranspiration/src/bmi_pet.c) - -add_library(cfelib ./extern/SoilMoistureProfiles/src/bmi_soil_moisture_profile.cxx - ./extern/SoilMoistureProfiles/src/soil_moisture_profile.cxx - ./extern/SoilMoistureProfiles/include/bmi_soil_moisture_profile.hxx - ./extern/SoilMoistureProfiles/include/soil_moisture_profile.hxx) +add_executable(${exe_name} ./src/main_cfe_aorc_pet_rz_aet.cxx ./src/cfe.c ./src/bmi_cfe.c ./src/logger.c ./src/giuh.c ./src/conceptual_reservoir.c + ./extern/aorc_bmi/src/aorc.c ./extern/aorc_bmi/src/bmi_aorc.c ./extern/evapotranspiration/src/pet.c + ./extern/evapotranspiration/src/bmi_pet.c ./src/bmi_serialization.cxx ./src/nash_cascade.c) + +add_library(cfelib ./extern/SoilMoistureProfiles/src/bmi_soil_moisture_profile.cxx ./extern/SoilMoistureProfiles/src/soil_moisture_profile.cxx + ./extern/SoilMoistureProfiles/include/bmi_soil_moisture_profile.hxx ./extern/SoilMoistureProfiles/include/soil_moisture_profile.hxx) target_link_libraries(${exe_name} LINK_PUBLIC cfelib) elseif(FORCING) -add_executable(${exe_name} ./src/main_pass_forcings.c ./src/cfe.c ./src/bmi_cfe.c ./src/giuh.c ./src/conceptual_reservoir.c - ./src/nash_cascade.c ./extern/aorc_bmi/src/aorc.c ./extern/aorc_bmi/src/bmi_aorc.c) +add_executable(${exe_name} ./src/main_pass_forcings.c ./src/cfe.c ./src/bmi_cfe.c ./src/logger.c ./src/giuh.c ./src/conceptual_reservoir.c + ./extern/aorc_bmi/src/aorc.c ./extern/aorc_bmi/src/bmi_aorc.c ./src/bmi_serialization.cxx ./src/nash_cascade.c) elseif(FORCINGPET) -add_executable(${exe_name} ./src/main_cfe_aorc_pet.c ./src/cfe.c ./src/bmi_cfe.c ./src/giuh.c ./src/conceptual_reservoir.c - ./src/nash_cascade.c ./extern/aorc_bmi/src/aorc.c ./extern/aorc_bmi/src/bmi_aorc.c - ./extern/evapotranspiration/src/pet.c ./extern/evapotranspiration/src/bmi_pet.c) +add_executable(${exe_name} ./src/main_cfe_aorc_pet.c ./src/cfe.c ./src/bmi_cfe.c ./src/logger.c ./src/giuh.c ./src/conceptual_reservoir.c + ./extern/aorc_bmi/src/aorc.c ./extern/aorc_bmi/src/bmi_aorc.c ./extern/evapotranspiration/src/pet.c + ./extern/evapotranspiration/src/bmi_pet.c ./src/bmi_serialization.cxx ./src/nash_cascade.c) elseif(BASE) -add_executable(${exe_name} ./src/main.c ./src/cfe.c ./src/bmi_cfe.c ./src/giuh.c ./src/conceptual_reservoir.c - ./src/nash_cascade.c) +add_executable(${exe_name} ./src/main.c ./src/cfe.c ./src/bmi_cfe.c ./src/logger.c ./src/giuh.c ./src/conceptual_reservoir.c + ./src/bmi_serialization.cxx ./src/nash_cascade.c) elseif(UNITTEST) -add_executable(${exe_name} ./test/main_unit_test.c ./src/cfe.c ./src/bmi_cfe.c ./src/giuh.c ./src/conceptual_reservoir.c - ./src/nash_cascade.c) +add_executable(${exe_name} ./test/main_unit_test.c ./src/cfe.c ./src/bmi_cfe.c ./src/logger.c ./src/giuh.c ./src/conceptual_reservoir.c + ./src/bmi_serialization.cxx ./src/nash_cascade.c) endif() +# ----------------------------------------------------------------------------- +# Find the Boost library and configure usage +set(Boost_USE_STATIC_LIBS OFF) +set(Boost_USE_MULTITHREADED ON) +set(Boost_USE_STATIC_RUNTIME OFF) +find_package(Boost 1.79.0 REQUIRED COMPONENTS serialization) + + + if(NOT NGEN) target_link_libraries(${exe_name} PRIVATE m) target_include_directories(${exe_name} PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/include) +target_link_libraries(${exe_name} PRIVATE Boost::serialization) endif() @@ -96,9 +110,11 @@ set(CFE_LIB_DESC_CMAKE "OWP CFE BMI Module Shared Library") add_compile_definitions(BMI_ACTIVE) if(WIN32) - add_library(cfebmi ./src/bmi_cfe.c ./src/cfe.c ./src/giuh.c ./src/conceptual_reservoir.c ./src/nash_cascade.c) + add_library(cfebmi ./src/bmi_cfe.c ./src/logger.c ./src/cfe.c ./src/giuh.c ./src/conceptual_reservoir.c + ./src/bmi_serialization.cxx) else() - add_library(cfebmi SHARED ./src/bmi_cfe.c ./src/cfe.c ./src/giuh.c ./src/conceptual_reservoir.c ./src/nash_cascade.c) + add_library(cfebmi SHARED ./src/bmi_cfe.c ./src/logger.c ./src/cfe.c ./src/giuh.c ./src/conceptual_reservoir.c + ./src/bmi_serialization.cxx) endif() target_include_directories(cfebmi PRIVATE include) @@ -107,6 +123,8 @@ set_target_properties(cfebmi PROPERTIES VERSION ${PROJECT_VERSION}) set_target_properties(cfebmi PROPERTIES PUBLIC_HEADER ./include/bmi_cfe.h) +target_link_libraries(cfebmi PRIVATE Boost::serialization) + include(GNUInstallDirs) install(TARGETS cfebmi diff --git a/LICENSE b/LICENSE index d6456956..aec7cc70 100644 --- a/LICENSE +++ b/LICENSE @@ -1,5 +1,28 @@ +Copyright 2025 Raytheon Company - Apache License +Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: +1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. +2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation +and/or other materials provided with the distribution. + +Licensed under: https://opensource.org/license/bsd-2-clause + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +All rights reserved. Based on Government sponsored work under contract GS-35F-204GA. + +----------------- + + +“Software code created by U.S. Government employees is not subject to copyright +in the United States (17 U.S.C. §105). The United States/Department of Commerce +reserve all rights to seek and obtain copyright protection in countries other +than the United States for Software authored in its entirety by the Department +of Commerce. To this end, the Department of Commerce hereby grants to Recipient +a royalty-free, nonexclusive license to use, copy, and create derivative works +of the Software outside of the United States.” + + Apache License Version 2.0, January 2004 http://www.apache.org/licenses/ @@ -199,4 +222,4 @@ distributed under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and - limitations under the License. + diff --git a/bmi/bmi.h b/bmi/bmi.h index 9cc04d64..e6bb6d12 100755 --- a/bmi/bmi.h +++ b/bmi/bmi.h @@ -35,8 +35,8 @@ SOFTWARE. extern "C" { #endif -const static int BMI_SUCCESS = 0; -const static int BMI_FAILURE = 1; +#define BMI_SUCCESS (0) +#define BMI_FAILURE (1) #define BMI_MAX_UNITS_NAME (2048) #define BMI_MAX_TYPE_NAME (2048) diff --git a/bmi/bmi.hxx b/bmi/bmi.hxx new file mode 100644 index 00000000..87b290b7 --- /dev/null +++ b/bmi/bmi.hxx @@ -0,0 +1,84 @@ +// The Basic Model Interface (BMI) C++ specification. +// +// This language specification is derived from the Scientific +// Interface Definition Language (SIDL) file bmi.sidl located at +// https://github.com/csdms/bmi. + +#ifndef BMI_HXX +#define BMI_HXX +#include +#include + +namespace bmixx { + + //const int BMI_SUCCESS = 0; + // const int BMI_FAILURE = 1; + + const int MAX_COMPONENT_NAME = 2048; + const int MAX_VAR_NAME = 2048; + const int MAX_UNITS_NAME = 2048; + const int MAX_TYPE_NAME = 2048; + + class Bmi { + public: + // Model control functions. + virtual void Initialize(std::string config_file) = 0; + virtual void Update() = 0; + virtual void UpdateUntil(double time) = 0; + virtual void Finalize() = 0; + + // Model information functions. + virtual std::string GetComponentName() = 0; + virtual int GetInputItemCount() = 0; + virtual int GetOutputItemCount() = 0; + virtual std::vector GetInputVarNames() = 0; + virtual std::vector GetOutputVarNames() = 0; + + // Variable information functions + virtual int GetVarGrid(std::string name) = 0; + virtual std::string GetVarType(std::string name) = 0; + virtual std::string GetVarUnits(std::string name) = 0; + virtual int GetVarItemsize(std::string name) = 0; + virtual int GetVarNbytes(std::string name) = 0; + virtual std::string GetVarLocation(std::string name) = 0; + + virtual double GetCurrentTime() = 0; + virtual double GetStartTime() = 0; + virtual double GetEndTime() = 0; + virtual std::string GetTimeUnits() = 0; + virtual double GetTimeStep() = 0; + + // Variable getters + virtual void GetValue(std::string name, void *dest) = 0; + virtual void *GetValuePtr(std::string name) = 0; + virtual void GetValueAtIndices(std::string name, void *dest, int *inds, int count) = 0; + + // Variable setters + virtual void SetValue(std::string name, void *src) = 0; + virtual void SetValueAtIndices(std::string name, int *inds, int count, void *src) = 0; + + // Grid information functions + virtual int GetGridRank(const int grid) = 0; + virtual int GetGridSize(const int grid) = 0; + virtual std::string GetGridType(const int grid) = 0; + + virtual void GetGridShape(const int grid, int *shape) = 0; + virtual void GetGridSpacing(const int grid, double *spacing) = 0; + virtual void GetGridOrigin(const int grid, double *origin) = 0; + + virtual void GetGridX(const int grid, double *x) = 0; + virtual void GetGridY(const int grid, double *y) = 0; + virtual void GetGridZ(const int grid, double *z) = 0; + + virtual int GetGridNodeCount(const int grid) = 0; + virtual int GetGridEdgeCount(const int grid) = 0; + virtual int GetGridFaceCount(const int grid) = 0; + + virtual void GetGridEdgeNodes(const int grid, int *edge_nodes) = 0; + virtual void GetGridFaceEdges(const int grid, int *face_edges) = 0; + virtual void GetGridFaceNodes(const int grid, int *face_nodes) = 0; + virtual void GetGridNodesPerFace(const int grid, int *nodes_per_face) = 0; + }; +} + +#endif diff --git a/configs/cfe1.0/cfe_config_cat_87.txt b/configs/cfe1.0/cfe_config_cat_87.txt index 87dae955..ff6729b0 100644 --- a/configs/cfe1.0/cfe_config_cat_87.txt +++ b/configs/cfe1.0/cfe_config_cat_87.txt @@ -20,13 +20,12 @@ nash_storage=0.0,0.0 giuh_ordinates=0.06,0.51,0.28,0.12,0.03 num_timesteps=1 verbosity=2 -surface_runoff_scheme=GIUH -surface_water_partitioning_scheme=Xinanjiang +surface_partitioning_scheme=Xinanjiang a_Xinanjiang_inflection_point_parameter=1 b_Xinanjiang_shape_parameter=1 x_Xinanjiang_shape_parameter=1 urban_decimal_fraction=0.0 DEBUG=0 -#surface_water_partitioning_scheme=Schaake +#surface_partitioning_scheme=Schaake #ice_fraction=0 #ice_content_threshold=0.15 diff --git a/configs/cfe1.0/cfe_config_laramie_pass_aet_rz.txt b/configs/cfe1.0/cfe_config_laramie_pass_aet_rz.txt index b6c92437..93f5e0b6 100644 --- a/configs/cfe1.0/cfe_config_laramie_pass_aet_rz.txt +++ b/configs/cfe1.0/cfe_config_laramie_pass_aet_rz.txt @@ -20,9 +20,8 @@ nash_storage=0.0,0.0 giuh_ordinates=0.06,0.51,0.28,0.12,0.03 num_timesteps=1 verbosity=1 -surface_runoff_scheme=GIUH -surface_water_partitioning_scheme=Schaake -#surface_water_partitioning_scheme=Xinanjiang +surface_partitioning_scheme=Schaake +#surface_partitioning_scheme=Xinanjiang #a_Xinanjiang_inflection_point_parameter=1 #b_Xinanjiang_shape_parameter=1 #x_Xinanjiang_shape_parameter=1 diff --git a/include/bmi.h b/include/bmi.h index 9cc04d64..e6bb6d12 100755 --- a/include/bmi.h +++ b/include/bmi.h @@ -35,8 +35,8 @@ SOFTWARE. extern "C" { #endif -const static int BMI_SUCCESS = 0; -const static int BMI_FAILURE = 1; +#define BMI_SUCCESS (0) +#define BMI_FAILURE (1) #define BMI_MAX_UNITS_NAME (2048) #define BMI_MAX_TYPE_NAME (2048) diff --git a/include/bmi_cfe.h b/include/bmi_cfe.h index 5e9c5d26..05c79ad6 100644 --- a/include/bmi_cfe.h +++ b/include/bmi_cfe.h @@ -7,6 +7,7 @@ extern "C" { #include "./cfe.h" #include "bmi.h" +#include //-------------------------------------------------- // Experiment to simplify BMI implementation (SDP) @@ -111,6 +112,9 @@ struct cfe_state_struct { int verbosity; + char* serialized; + uint64_t serialized_length; // needs a permanent anchor for get_value_ptr + }; typedef struct cfe_state_struct cfe_state_struct; diff --git a/include/bmi_serialization.h b/include/bmi_serialization.h new file mode 100644 index 00000000..f0f6c600 --- /dev/null +++ b/include/bmi_serialization.h @@ -0,0 +1,18 @@ +#ifndef CFE_BMI_SERIALIZATION +#define CFE_BMI_SERIALIZATION + +#ifdef __cplusplus +extern "C" { +#endif + +#include "bmi.h" + +int free_serialized_cfe(Bmi* bmi); +int load_serialized_cfe(Bmi* bmi, const char* data); +int new_serialized_cfe(Bmi* bmi); + +#ifdef __cplusplus +} +#endif + +#endif diff --git a/include/cfe.h b/include/cfe.h index 045f7deb..3b93d6d3 100644 --- a/include/cfe.h +++ b/include/cfe.h @@ -17,7 +17,7 @@ // t-shirt approximation of the hydrologic routing funtionality of the National Water Model v 1.2, 2.0, and 2.1 // This code was developed to test the hypothesis that the National Water Model runoff generation, vadose zone -// dynamics, and conceptual groundwater model can be greatly simplified by acknowledging that it is truly a +// dynamics, and conceptual groundwater model can be greatly simplified by acknowledging that it is truly a // conceptual model. The hypothesis is supported by a number of observations made during a 2017-2018 deep dive // into the NWM code. Thesed are: // @@ -26,7 +26,7 @@ // function by Moore, 1985. The Schaake function is a single valued function of soil moisture deficit, // predicts 100% runoff when the soil is saturated, like the curve-number method, and is fundamentally simple. // 2. Run-on infiltration is strictly not calculated. Overland flow routing applies the Schaake function repeatedly -// to predict this phenomenon, which violates the underlying assumption of the PDM method that only rainfall +// to predict this phenomenon, which violates the underlying assumption of the PDM method that only rainfall // inputs affect soil moisture. // 3. The water-content based Richards' equation, applied using a coarse-discretization, can be replaced with a simple // conceptual reservoir because it never allows saturation or infiltration-excess runoff unless deactivated by @@ -118,8 +118,8 @@ struct massbal double vol_soil_to_gw ; // this should equal vol_to_gw double vol_soil_end ; double vol_et_from_soil ; - double vol_et_from_rain ; - double vol_et_to_atm ; + double vol_et_from_rain ; + double vol_et_to_atm ; double volin ; double volout ; double volend ; diff --git a/include/giuh.h b/include/giuh.h index 849f9b2f..82ce2f34 100644 --- a/include/giuh.h +++ b/include/giuh.h @@ -2,6 +2,6 @@ #define _GIUH_H extern double giuh_convolution_integral(double runoff_m, int num_giuh_ordinates, - double *giuh_ordinates, double *runoff_queue_m_per_timestep); + double *giuh_ordinates, double *runoff_queue_m_per_timestep); #endif diff --git a/include/logger.h b/include/logger.h new file mode 100644 index 00000000..877d9a8d --- /dev/null +++ b/include/logger.h @@ -0,0 +1,18 @@ +#ifndef LOGGER_H +#define LOGGER_H + +#include // for variable args: va_list + +typedef enum { + NONE = 0, + DEBUG = 1, + INFO = 2, + WARNING = 3, + SEVERE = 4, + FATAL = 5, +} LogLevel; + +// Public Methods +void Log(LogLevel messageLevel, const char* message, ...); + +#endif // LOGGER_H diff --git a/include/vecbuf.hpp b/include/vecbuf.hpp new file mode 100644 index 00000000..6db982cc --- /dev/null +++ b/include/vecbuf.hpp @@ -0,0 +1,115 @@ +#ifndef HPP_STRING_VECBUF +#define HPP_STRING_VECBUF +// https://gist.github.com/stephanlachnit/4a06f8475afd144e73235e2a2584b000 +// SPDX-FileCopyrightText: 2023 Stephan Lachnit +// SPDX-License-Identifier: MIT + +#include +#include +#include + +template> +class vecbuf : public std::basic_streambuf { +public: + using streambuf = std::basic_streambuf; + using char_type = typename streambuf::char_type; + using int_type = typename streambuf::int_type; + using traits_type = typename streambuf::traits_type; + using vector = std::vector; + using value_type = typename vector::value_type; + using size_type = typename vector::size_type; + + // Constructor for vecbuf with optional initial capacity + vecbuf(size_type capacity = 0) : vector_() { reserve(capacity); } + + // Forwarder for std::vector::shrink_to_fit() + constexpr void shrink_to_fit() { vector_.shrink_to_fit(); } + + // Forwarder for std::vector::clear() + constexpr void clear() { vector_.clear(); } + + // Forwarder for std::vector::reserve + constexpr void reserve(size_type capacity) { vector_.reserve(capacity); setp_from_vector(); } + + // Increase the capacity of the buffer by reserving the current_size + additional_capacity + constexpr void reserve_additional(size_type additional_capacity) { reserve(size() + additional_capacity); } + + // Forwarder for std::vector::data + constexpr const value_type* data() const { return vector_.data(); } + + // Forwarder for std::vector::size + constexpr size_type size() const { return vector_.size(); } + + // Forwarder for std::vector::capacity + constexpr size_type capacity() const { return vector_.capacity(); } + + // Implements std::basic_streambuf::xsputn + std::streamsize xsputn(const char_type* s, std::streamsize count) override { + try { + reserve_additional(count); + } + catch (const std::bad_alloc& error) { + // reserve did not work, use slow algorithm + return xsputn_slow(s, count); + } + // reserve worked, use fast algorithm + return xsputn_fast(s, count); + } + +protected: + // Calculates value to std::basic_streambuf::pbase from vector + constexpr value_type* pbase_from_vector() const { return const_cast(vector_.data()); } + + // Calculates value to std::basic_streambuf::pptr from vector + constexpr value_type* pptr_from_vector() const { return const_cast(vector_.data() + vector_.size()); } + + // Calculates value to std::basic_streambuf::epptr from vector + constexpr value_type* epptr_from_vector() const { return const_cast(vector_.data()) + vector_.capacity(); } + + // Sets the values for std::basic_streambuf::pbase, std::basic_streambuf::pptr and std::basic_streambuf::epptr from vector + constexpr void setp_from_vector() { streambuf::setp(pbase_from_vector(), epptr_from_vector()); streambuf::pbump(size()); } + +private: + // std::vector containing the data + vector vector_; + + // Fast implementation of std::basic_streambuf::xsputn if reserve_additional(count) succeeded + std::streamsize xsputn_fast(const char_type* s, std::streamsize count) { + // store current pptr (end of vector location) + auto* old_pptr = pptr_from_vector(); + // resize the vector, does not move since space already reserved + vector_.resize(vector_.size() + count); + // directly memcpy new content to old pptr (end of vector before it was resized) + traits_type::copy(old_pptr, s, count); + // reserve() already calls setp_from_vector(), only adjust pptr to new epptr + streambuf::pbump(count); + + return count; + } + + // Slow implementation of std::basic_streambuf::xsputn if reserve_additional(count) did not succeed, might calls std::basic_streambuf::overflow() + std::streamsize xsputn_slow(const char_type* s, std::streamsize count) { + // reserving entire vector failed, emplace char for char + std::streamsize written = 0; + while (written < count) { + try { + // copy one char, should throw eventually std::bad_alloc + vector_.emplace_back(s[written]); + } + catch (const std::bad_alloc& error) { + // try overflow(), if eof return, else continue + int_type c = this->overflow(traits_type::to_int_type(s[written])); + if (traits_type::eq_int_type(c, traits_type::eof())) { + return written; + } + } + // update pbase, pptr and epptr + setp_from_vector(); + written++; + } + return written; + } + +}; + +#endif diff --git a/realizations/realization_cfe_pet_surfgiuh.json b/realizations/realization_cfe_pet.json similarity index 100% rename from realizations/realization_cfe_pet_surfgiuh.json rename to realizations/realization_cfe_pet.json diff --git a/src/Makefile b/src/Makefile index b78e3ebe..4c17e172 100755 --- a/src/Makefile +++ b/src/Makefile @@ -13,7 +13,7 @@ LFLAGS = LIBS = -lm # define the C source files -SRCS = main.c cfe.c bmi_cfe.c +SRCS = main.c cfe.c bmi_cfe.c logger.c # define the C object files OBJS = $(SRCS:.c=.o) diff --git a/src/bmi_cfe.c b/src/bmi_cfe.c index 6242ec36..bd0f4edc 100644 --- a/src/bmi_cfe.c +++ b/src/bmi_cfe.c @@ -5,6 +5,9 @@ #include "bmi_cfe.h" #include #include +#include "logger.h" +#include "bmi_serialization.h" + #ifndef WATER_SPECIFIC_WEIGHT #define WATER_SPECIFIC_WEIGHT 9810 #define STANDARD_ATMOSPHERIC_PRESSURE_PASCALS 101325 @@ -28,7 +31,7 @@ #define PARAM_VAR_NAME_COUNT 18 // NOTE: If you update the params, also update the unit test in ../test/main_unit_test_bmi.c static const char *param_var_names[PARAM_VAR_NAME_COUNT] = { - "maxsmc", "satdk", "slope", "b", "Klf", + "maxsmc", "satdk", "slope", "b", "Klf", "Kn", "Cgw", "expon", "max_gw_storage", "satpsi","wltsmc","alpha_fc","refkdt", "a_Xinanjiang_inflection_point_parameter","b_Xinanjiang_shape_parameter","x_Xinanjiang_shape_parameter", @@ -177,8 +180,8 @@ Variable var_info[] = { // Root zone adjusted AET development -rlm -AJ // ------------------------------------------- { 91, "soil_moisture_profile", "double", 1}, - { 92, "soil_layer_depths_m", "double", 1}, - { 93, "max_rootzone_layer", "int", 1}, + { 92, "soil_layer_depths_m", "double", 1}, + { 93, "max_rootzone_layer", "int", 1}, //-------------------------------------------- { 94, "nwm_ponded_depth", "double", 1}, @@ -195,7 +198,7 @@ static const char *output_var_names[OUTPUT_VAR_NAME_COUNT] = { "DIRECT_RUNOFF", "NASH_LATERAL_RUNOFF", "DEEP_GW_TO_CHANNEL_FLUX", - "SOIL_TO_GW_FLUX", + "SOIL_TO_GW_FLUX", "Q_OUT", "POTENTIAL_ET", "ACTUAL_ET", @@ -203,7 +206,7 @@ static const char *output_var_names[OUTPUT_VAR_NAME_COUNT] = { "SOIL_STORAGE", "SOIL_STORAGE_CHANGE", "SURF_RUNOFF_SCHEME", - "NWM_PONDED_DEPTH" + "NWM_PONDED_DEPTH" }; static const char *output_var_types[OUTPUT_VAR_NAME_COUNT] = { @@ -218,10 +221,10 @@ static const char *output_var_types[OUTPUT_VAR_NAME_COUNT] = { "double", "double", "double", - "double", - "double", - "int", - "double" + "double", + "double", + "int", + "double" }; static const int output_var_item_count[OUTPUT_VAR_NAME_COUNT] = { @@ -232,14 +235,14 @@ static const int output_var_item_count[OUTPUT_VAR_NAME_COUNT] = { 1, 1, 1, - 1, - 1, + 1, + 1, 1, 1, 1, - 1, - 1, - 1 + 1, + 1, + 1 }; static const char *output_var_units[OUTPUT_VAR_NAME_COUNT] = { @@ -254,10 +257,10 @@ static const char *output_var_units[OUTPUT_VAR_NAME_COUNT] = { "m", "m", "m", - "m", - "m", - "none", - "m" + "m", + "m", + "1", + "m" }; static const int output_var_grids[OUTPUT_VAR_NAME_COUNT] = { @@ -273,9 +276,9 @@ static const int output_var_grids[OUTPUT_VAR_NAME_COUNT] = { 0, 0, 0, - 0, - 0, - 0 + 0, + 0, + 0 }; static const char *output_var_locations[OUTPUT_VAR_NAME_COUNT] = { @@ -316,9 +319,9 @@ static const char *input_var_types[INPUT_VAR_NAME_COUNT] = { static const char *input_var_units[INPUT_VAR_NAME_COUNT] = { "mm h-1", //"atmosphere_water__liquid_equivalent_precipitation_rate" "m s-1", //"water_potential_evaporation_flux" - "m", // ice fraction in meters - "none", // ice fraction [-] - "none" // soil moisture profile is in decimal fraction -rlm + "1", // ice fraction should be dimensionless + "1", // ice fraction [-] + "1" // soil moisture profile is in decimal fraction -rlm }; static const int input_var_item_count[INPUT_VAR_NAME_COUNT] = { @@ -332,17 +335,17 @@ static const int input_var_item_count[INPUT_VAR_NAME_COUNT] = { static const char input_var_grids[INPUT_VAR_NAME_COUNT] = { 0, 0, - 0, - 0, - 0 + 0, + 0, + 0 }; static const char *input_var_locations[INPUT_VAR_NAME_COUNT] = { "node", "node", - "node", - "node", - "node" + "node", + "node", + "node" }; static int Get_start_time (Bmi *self, double * time) @@ -362,7 +365,7 @@ static int Get_end_time (Bmi *self, double * time) cfe_state_struct *cfe; cfe = (cfe_state_struct *) self->data; Get_start_time(self, time); - + // see if forcings read in or via BMI (framework to set) //jmframe: the docs say that "If the model doesn’t define an end time, a large number" // so even if it is BMI, we should still use num_timesteps @@ -400,9 +403,8 @@ static int Get_time_units (Bmi *self, char * units) static int Get_current_time (Bmi *self, double * time) { Get_start_time(self, time); -#if CFE_DEBUG > 1 - printf("Current model time step: '%d'\n", ((cfe_state_struct *) self->data)->current_time_step); -#endif + Log(DEBUG, "Current model time step: '%d'\n", ((cfe_state_struct *) self->data)->current_time_step); + *time += (((cfe_state_struct *) self->data)->current_time_step * ((cfe_state_struct *) self->data)->time_step_size); return BMI_SUCCESS; } @@ -431,12 +433,11 @@ int read_init_config_cfe(const char* config_file, cfe_state_struct* model) // Note that this determines max line length including the ending return character, if present int count_result = read_file_line_counts_cfe(config_file, &config_line_count, &max_config_line_length); if (count_result == -1) { - printf("Invalid config file '%s'", config_file); + Log(WARNING, "Invalid config file '%s'", config_file); return BMI_FAILURE; } -#if CFE_DEBUG >= 1 - printf("Config file details - Line Count: %d | Max Line Length %d\n", config_line_count, max_config_line_length); -#endif + + Log(DEBUG, "Config file details - Line Count: %d | Max Line Length %d\n", config_line_count, max_config_line_length); FILE* fp = fopen(config_file, "r"); if (fp == NULL) @@ -546,9 +547,9 @@ int read_init_config_cfe(const char* config_file, cfe_state_struct* model) char *param_key, *param_value, *param_units; if (fgets(config_line, max_config_line_length + 1, fp) == NULL) return BMI_FAILURE; -#if CFE_DEBUG >= 3 - printf("Line value: ['%s']\n", config_line); -#endif + + Log(DEBUG, "Line value: ['%s']\n", config_line); + char* config_line_ptr = config_line; config_line_ptr = strsep(&config_line_ptr, "\n"); param_key = strsep(&config_line_ptr, "="); @@ -556,9 +557,7 @@ int read_init_config_cfe(const char* config_file, cfe_state_struct* model) param_value = strsep(&config_line_ptr, "["); param_units = strsep(&config_line_ptr, "]"); -#if CFE_DEBUG >= 1 - printf("Config Value - Param: '%s' | Value: '%s' | Units: '%s'\n", param_key, param_value, param_units); -#endif + Log(DEBUG, "Config Value - Param: '%s' | Value: '%s' | Units: '%s'\n", param_key, param_value, param_units); if (strcmp(param_key, "forcing_file") == 0) { model->forcing_file = strdup(param_value); @@ -574,9 +573,7 @@ int read_init_config_cfe(const char* config_file, cfe_state_struct* model) is_soil_params__depth_set = TRUE; // Check if units are present and print warning if missing from config file if ((param_units == NULL) || (strlen(param_units) < 1)) { -#if CFE_DEBUG >= 1 - printf ("WARNING: [units] expected for '%s' in config file \n", param_key); -#endif + Log(WARNING, "[units] expected for '%s' in config file \n", param_key); } continue; } @@ -590,9 +587,7 @@ int read_init_config_cfe(const char* config_file, cfe_state_struct* model) is_soil_params__satdk_set = TRUE; // Check if units are present and print warning if missing from config file if ((param_units == NULL) || (strlen(param_units) < 1)) { -#if CFE_DEBUG >= 1 - printf ("WARNING: [units] expected for '%s' in config file \n", param_key); -#endif + Log(WARNING, "[units] expected for '%s' in config file \n", param_key); } continue; } @@ -601,9 +596,7 @@ int read_init_config_cfe(const char* config_file, cfe_state_struct* model) is_soil_params__satpsi_set = TRUE; // Check if units are present and print warning if missing from config file if ((param_units == NULL) || (strlen(param_units) < 1)) { -#if CFE_DEBUG >= 1 - printf ("WARNING: [units] expected for '%s' in config file \n", param_key); -#endif + Log(WARNING, "[units] expected for '%s' in config file \n", param_key); } continue; } @@ -612,9 +605,7 @@ int read_init_config_cfe(const char* config_file, cfe_state_struct* model) is_soil_params__slop_set = TRUE; // Check if units are present and print warning if missing from config file if ((param_units == NULL) || (strlen(param_units) < 1)) { -#if CFE_DEBUG >= 1 - printf ("WARNING: [units] expected for '%s' in config file \n", param_key); -#endif + Log(WARNING, "[units] expected for '%s' in config file \n", param_key); } continue; } @@ -623,9 +614,7 @@ int read_init_config_cfe(const char* config_file, cfe_state_struct* model) is_soil_params__smcmax_set = TRUE; // Check if units are present and print warning if missing from config file if ((param_units == NULL) || (strlen(param_units) < 1)) { -#if CFE_DEBUG >= 1 - printf ("WARNING: [units] expected for '%s' in config file \n", param_key); -#endif + Log(WARNING, "[units] expected for '%s' in config file \n", param_key); } continue; } @@ -634,9 +623,7 @@ int read_init_config_cfe(const char* config_file, cfe_state_struct* model) is_soil_params__wltsmc_set = TRUE; // Check if units are present and print warning if missing from config file if ((param_units == NULL) || (strlen(param_units) < 1)) { -#if CFE_DEBUG >= 1 - printf ("WARNING: [units] expected for '%s' in config file \n", param_key); -#endif + Log(WARNING, "[units] expected for '%s' in config file \n", param_key); } continue; } @@ -655,9 +642,7 @@ int read_init_config_cfe(const char* config_file, cfe_state_struct* model) is_gw_max_set = TRUE; // Check if units are present and print warning if missing from config file if ((param_units == NULL) || (strlen(param_units) < 1)) { -#if CFE_DEBUG >= 1 - printf ("WARNING: [units] expected for '%s' in config file \n", param_key); -#endif + Log(SEVERE, "[units] expected for '%s' in config file \n", param_key); } // Also set the true storage if storage was already read and was a ratio, and so we were waiting for this /* if (is_gw_storage_set == TRUE && is_gw_storage_ratio == TRUE) { @@ -670,9 +655,7 @@ int read_init_config_cfe(const char* config_file, cfe_state_struct* model) is_Cgw_set = TRUE; // Check if units are present and print warning if missing from config file if ((param_units == NULL) || (strlen(param_units) < 1)) { -#if CFE_DEBUG >= 1 - printf ("WARNING: [units] expected for '%s' in config file \n", param_key); -#endif + Log(SEVERE, "[units] expected for '%s' in config file \n", param_key); } continue; } @@ -683,13 +666,13 @@ int read_init_config_cfe(const char* config_file, cfe_state_struct* model) } if (strcmp(param_key, "gw_storage") == 0) { model->gw_reservoir.gw_storage = strtod(param_value, NULL); - model->gw_reservoir.storage_m = model->gw_reservoir.gw_storage * model->gw_reservoir.storage_max_m; //edited by RLM to fix units from [m/m] to [m] + // YLiu: move the calculation of GW storage out of the FOR loop so that the position of gw_storage + // relative to that of max_gw_storage in the config file would not matter + //model->gw_reservoir.storage_m = model->gw_reservoir.gw_storage * model->gw_reservoir.storage_max_m; //edited by RLM to fix units from [m/m] to [m] is_gw_storage_set = TRUE; // Check if units are present and print warning if missing from config file if ((param_units == NULL) || (strlen(param_units) < 1)) { -#if CFE_DEBUG >= 1 - printf ("WARNING: [units] expected for '%s' in config file \n", param_key); -#endif + Log(WARNING, "[units] expected for '%s' in config file \n", param_key); } /* char* trailing_chars; gw_storage_literal = strtod(param_value, &trailing_chars); @@ -711,7 +694,7 @@ int read_init_config_cfe(const char* config_file, cfe_state_struct* model) is_alpha_fc_set = TRUE; continue; } - + if (strcmp(param_key, "soil_storage") == 0) { /* char* trailing_chars; double parsed_value = strtod(param_value, &trailing_chars); @@ -721,9 +704,7 @@ int read_init_config_cfe(const char* config_file, cfe_state_struct* model) is_soil_storage_set = TRUE; // Check if units are present and print warning if missing from config file if ((param_units == NULL) || (strlen(param_units) < 1)) { -#if CFE_DEBUG >= 1 - printf ("WARNING: [units] expected for '%s' in config file \n", param_key); -#endif + Log(WARNING, "[units] expected for '%s' in config file \n", param_key); } continue; } @@ -736,9 +717,7 @@ int read_init_config_cfe(const char* config_file, cfe_state_struct* model) is_K_nash_subsurface_set = TRUE; // Check if units are present and print warning if missing from config file if ((param_units == NULL) || (strlen(param_units) < 1)) { -#if CFE_DEBUG >= 1 - printf ("WARNING: [units] expected for '%s' in config file. [1/hour]\n", param_key); -#endif + Log(WARNING, "[units] expected for '%s' in config file. [1/hour]\n", param_key); } continue; } @@ -753,9 +732,7 @@ int read_init_config_cfe(const char* config_file, cfe_state_struct* model) continue; } if (strcmp(param_key, "giuh_ordinates") == 0) { -#if CFE_DEBUG >= 1 - printf("Found configured GIUH ordinate values ('%s')\n", param_value); -#endif + Log(DEBUG, "Found configured GIUH ordinate values ('%s')\n", param_value); giuh_originates_string_val = strdup(param_value); is_giuh_originates_string_val_set = TRUE; continue; @@ -770,21 +747,19 @@ int read_init_config_cfe(const char* config_file, cfe_state_struct* model) is_verbosity_set = TRUE; continue; } - + /*-------------------- Root zone AET development -rlm -----------------------*/ if (strcmp(param_key, "is_aet_rootzone") == 0) { if ( strcmp(param_value, "true")==0 || strcmp(param_value, "True")==0 || strcmp(param_value,"1")==0) is_aet_rootzone_set = TRUE; - + continue; } if (is_aet_rootzone_set == TRUE) { if (strcmp(param_key, "soil_layer_depths") == 0) { -#if CFE_DEBUG >= 1 - printf("Found configured soil depth values ('%s')\n", param_value); -#endif + Log(DEBUG, "Found configured soil depth values ('%s')\n", param_value); soil_layer_depths_string_val = strdup(param_value); is_soil_layer_depths_string_val_set = TRUE; continue; @@ -816,9 +791,7 @@ int read_init_config_cfe(const char* config_file, cfe_state_struct* model) model->nash_surface_params.K_nash = strtod(param_value, NULL); is_K_nash_surface_set = TRUE; if ((param_units == NULL) || (strlen(param_units) < 1)) { -#if CFE_DEBUG >= 1 - printf ("WARNING: [units] expected for '%s' in config file. [1/hour]\n", param_key); -#endif + Log(DEBUG, "[units] expected for '%s' in config file. [1/hour]\n", param_key); } continue; } @@ -831,9 +804,7 @@ int read_init_config_cfe(const char* config_file, cfe_state_struct* model) nash_storage_surface_string_val = strdup(param_value); is_nash_storage_surface_set = TRUE; if ((param_units == NULL) || (strlen(param_units) < 1)) { -#if CFE_DEBUG >= 1 - printf ("WARNING: [units] expected for '%s' in config file. [m]\n", param_key); -#endif + Log(WARNING, "[units] expected for '%s' in config file. [m]\n", param_key); } continue; } @@ -841,9 +812,7 @@ int read_init_config_cfe(const char* config_file, cfe_state_struct* model) model->nash_surface_params.K_infiltration = strtod(param_value, NULL); is_K_infiltration_nash_surface_set = TRUE; if ((param_units == NULL) || (strlen(param_units) < 1)) { -#if CFE_DEBUG >= 1 - printf ("WARNING: [units] expected for '%s' in config file. [1/hour]\n", param_key); -#endif + Log(WARNING, "[units] expected for '%s' in config file. [1/hour]\n", param_key); } continue; } @@ -851,9 +820,7 @@ int read_init_config_cfe(const char* config_file, cfe_state_struct* model) model->nash_surface_params.retention_depth = strtod(param_value, NULL); is_retention_depth_nash_surface_set = TRUE; if ((param_units == NULL) || (strlen(param_units) < 1)) { -#if CFE_DEBUG >= 1 - printf ("WARNING: [units] expected for '%s' in config file. [m]\n", param_key); -#endif + Log(WARNING, "[units] expected for '%s' in config file. [m]\n", param_key); } continue; } @@ -870,7 +837,7 @@ int read_init_config_cfe(const char* config_file, cfe_state_struct* model) is_infiltration_excess_method_set = TRUE; if (strcmp(param_key, "surface_partitioning_scheme") == 0) { - printf("WARNING: The name \"surface_partitioning_scheme\" in the config file will be deprecated; use the new name \"surface_water_partitioning_scheme\" instead. \n"); + Log(WARNING, "The name \"surface_partitioning_scheme\" in the config file will be deprecated; use the new name \"surface_water_partitioning_scheme\" instead. \n"); } continue; } @@ -898,7 +865,7 @@ int read_init_config_cfe(const char* config_file, cfe_state_struct* model) if (strcmp(param_key, "is_sft_coupled") == 0) { if ( strcmp(param_value, "true")==0 || strcmp(param_value, "True")==0 || strcmp(param_value,"1")==0) is_sft_coupled_set = TRUE; - + continue; } @@ -908,178 +875,132 @@ int read_init_config_cfe(const char* config_file, cfe_state_struct* model) is_ice_content_threshold_set = TRUE; // Check if units are present and print warning if missing from config file if ((param_units == NULL) || (strlen(param_units) < 1)) { -#if CFE_DEBUG >= 1 - printf ("WARNING: [units] expected for '%s' in config file \n", param_key); -#endif + Log(WARNING, "[units] expected for '%s' in config file \n", param_key); } } } } - + if (is_forcing_file_set == FALSE) { -#if CFE_DEBUG >= 1 - printf("Config param 'forcing_file' not found in config file\n"); -#endif - return BMI_FAILURE; + Log(SEVERE, "Config param 'forcing_file' not found in config file. Aborting...\n"); + return BMI_FAILURE; } if (is_soil_params__depth_set == FALSE) { -#if CFE_DEBUG >= 1 - printf("Config param 'soil_params.depth' not found in config file\n"); -#endif - return BMI_FAILURE; + Log(SEVERE, "Config param 'soil_params.depth' not found in config file. Aborting...\n"); + return BMI_FAILURE; } if (is_soil_params__bb_set == FALSE) { -#if CFE_DEBUG >= 1 - printf("Config param 'soil_params.bb' not found in config file\n"); -#endif - return BMI_FAILURE; + Log(SEVERE, "Config param 'soil_params.bb' not found in config file. Aborting...\n"); + return BMI_FAILURE; } - + if (is_soil_params__satdk_set == FALSE) { -#if CFE_DEBUG >= 1 - printf("Config param 'soil_params.satdk' not found in config file\n"); -#endif - return BMI_FAILURE; + Log(SEVERE, "Config param 'soil_params.satdk' not found in config file. Aborting...\n"); + return BMI_FAILURE; } if (is_soil_params__satpsi_set == FALSE) { -#if CFE_DEBUG >= 1 - printf("Config param 'soil_params.satpsi' not found in config file\n"); -#endif - return BMI_FAILURE; + Log(SEVERE, "Config param 'soil_params.satpsi' not found in config file. Aborting...\n"); + return BMI_FAILURE; } if (is_soil_params__slop_set == FALSE) { -#if CFE_DEBUG >= 1 - printf("Config param 'soil_params.slop' not found in config file\n"); -#endif - return BMI_FAILURE; + Log(WARNING, "Config param 'soil_params.slop' not found in config file. Aborting...\n"); + return BMI_FAILURE; } if (is_soil_params__smcmax_set == FALSE) { -#if CFE_DEBUG >= 1 - printf("Config param 'soil_params.smcmax' not found in config file\n"); -#endif - return BMI_FAILURE; + Log(SEVERE, "Config param 'soil_params.smcmax' not found in config file. Aborting...\n"); + return BMI_FAILURE; } if (is_soil_params__wltsmc_set == FALSE) { -#if CFE_DEBUG >= 1 - printf("Config param 'soil_params.wltsmc' not found in config file\n"); -#endif - return BMI_FAILURE; + Log(SEVERE, "Config param 'soil_params.wltsmc' not found in config file. Aborting...\n"); + return BMI_FAILURE; } if (is_soil_params__expon_set == FALSE) { -#if CFE_DEBUG >= 1 - printf("Config param 'soil_params.expon' not found in config file, defaulting to 1 (linear)\n"); -#endif - model->soil_reservoir.exponent_primary = 1.0; - //is_soil_params__expon_set == TRUE; - // Don't return BMI_FAILURE, this is a optional config - //return BMI_FAILURE; + Log(WARNING, "Config param 'soil_params.expon' not found in config file, defaulting to 1 (linear)\n"); + model->soil_reservoir.exponent_primary = 1.0; + //is_soil_params__expon_set == TRUE; + // Don't return BMI_FAILURE, this is a optional config + //return BMI_FAILURE; } if (is_soil_params__expon2_set == FALSE) { -#if CFE_DEBUG >= 1 - printf("Config param 'soil_params.expon_secondary' not found in config file, defaulting to 1 (linear)\n"); -#endif - model->soil_reservoir.exponent_secondary = 1.0; - // Don't return BMI_FAILURE, this is a optional config - //return BMI_FAILURE; + Log(WARNING, "Config param 'soil_params.expon_secondary' not found in config file, defaulting to 1 (linear)\n"); + model->soil_reservoir.exponent_secondary = 1.0; + // Don't return BMI_FAILURE, this is a optional config + //return BMI_FAILURE; } if (is_Cgw_set == FALSE) { -#if CFE_DEBUG >= 1 - printf("Config param 'Cgw' not found in config file. Aborting...\n"); -#endif - return BMI_FAILURE; + Log(SEVERE, "Config param 'Cgw' not found in config file. Aborting...\n"); + return BMI_FAILURE; } if (is_expon_set == FALSE) { -#if CFE_DEBUG >= 1 - printf("Config param 'expon' not found in config file. Aborting...\n"); -#endif - return BMI_FAILURE; + Log(WARNING, "Config param 'expon' not found in config file. Aborting...\n"); + return BMI_FAILURE; } if (is_alpha_fc_set == FALSE) { -#if CFE_DEBUG >= 1 - printf("Config param 'alpha_fc' not found in config file. Aborting...\n"); -#endif - return BMI_FAILURE; + Log(SEVERE, "Config param 'alpha_fc' not found in config file. Aborting...\n"); + return BMI_FAILURE; } if (is_soil_storage_set == FALSE) { -#if CFE_DEBUG >= 1 - printf("Config param 'soil_storage' not found in config file. Aborting...\n"); -#endif - return BMI_FAILURE; + Log(SEVERE, "Config param 'soil_storage' not found in config file. Aborting...\n"); + return BMI_FAILURE; } if (is_K_nash_subsurface_set == FALSE) { -#if CFE_DEBUG >= 1 - printf("Config param 'K_nash_subsurface' not found in config file. Aborting...\n"); -#endif - return BMI_FAILURE; + Log(SEVERE, "Config param 'K_nash' not found in config file. Aborting...\n"); + return BMI_FAILURE; } if (is_K_lf_set == FALSE) { -#if CFE_DEBUG >= 1 - printf("Config param 'K_lf' not found in config file. Aborting...\n"); -#endif - return BMI_FAILURE; + Log(SEVERE, "Config param 'K_lf' not found in config file. Aborting...\n"); + return BMI_FAILURE; } if (is_gw_max_set == FALSE) { -#if CFE_DEBUG >= 1 - printf("Config param 'max_gw_storage' not found in config file. Aborting...\n"); -#endif - return BMI_FAILURE; + Log(SEVERE, "Config param 'max_gw_storage' not found in config file. Aborting...\n"); + return BMI_FAILURE; } if (is_gw_storage_set == FALSE) { -#if CFE_DEBUG >= 1 - printf("Config param 'gw_storage' not found in config file. Aborting...\n"); -#endif - return BMI_FAILURE; + Log(SEVERE, "Config param 'gw_storage' not found in config file. Aborting...n"); + return BMI_FAILURE; + } + + // compute gw storage in meters + if ((is_gw_storage_set == TRUE) && (is_gw_max_set == TRUE)) { + model->gw_reservoir.storage_m = model->gw_reservoir.gw_storage * model->gw_reservoir.storage_max_m; } + if (is_num_timesteps_set == FALSE && strcmp(model->forcing_file, "BMI")) { -#if CFE_DEBUG >= 1 - printf("Config param 'num_timesteps' not found in config file. Aborting...\n"); -#endif - return BMI_FAILURE; + Log(SEVERE, "Config param 'num_timesteps' not found in config file. Aborting...\n"); + return BMI_FAILURE; } if (is_verbosity_set == FALSE) { -#if CFE_DEBUG >= 1 - printf("Config param 'verbosity' not found in config file (optional) \n"); - printf("setting verbosity to a low value\n"); -#endif - model->verbosity = 0; + Log(WARNING, "Config param 'verbosity' not found in config file\n"); + Log(WARNING, "setting verbosity to a low value\n"); + model->verbosity = 0; } if (is_infiltration_excess_method_set == FALSE) { -#if CFE_DEBUG >= 1 - printf("Config param \"surface_water_partitioning_scheme\" not found in config file. Aborting...\n"); -#endif - return BMI_FAILURE; + Log(WARNING, "Config param 'infiltration_excess_method' not found in config file. Aborting...\n"); + return BMI_FAILURE; } - /* xinanjiang_dev*/ +/* xinanjiang_dev*/ if(model->infiltration_excess_params_struct.surface_water_partitioning_scheme == Xinanjiang){ - if (is_a_Xinanjiang_inflection_point_parameter_set == FALSE) { -#if CFE_DEBUG >= 1 - printf("Config param 'a_Xinanjiang_inflection_point_parameter' not found in config file. Aborting...\n"); -#endif - return BMI_FAILURE; - } - if (is_b_Xinanjiang_shape_parameter_set == FALSE) { -#if CFE_DEBUG >= 1 - printf("Config param 'b_Xinanjiang_shape_parameter' not found in config file. Aborting...\n"); -#endif - return BMI_FAILURE; - } - if (is_x_Xinanjiang_shape_parameter_set == FALSE) { -#if CFE_DEBUG >= 1 - printf("Config param 'x_Xinanjiang_shape_parameter' not found in config file. Aborting...\n"); -#endif - return BMI_FAILURE; - } - if (is_urban_decimal_fraction_set == FALSE) { -#if CFE_DEBUG >= 1 - printf("Config param 'urban_decimal_fraction' not found in config file. Aborting...\n"); -#endif - return BMI_FAILURE; - } + if (is_a_Xinanjiang_inflection_point_parameter_set == FALSE) { + Log(SEVERE, "Config param 'a_Xinanjiang_inflection_point_parameter' not found in config file. Aborting...\n"); + return BMI_FAILURE; + } + if (is_b_Xinanjiang_shape_parameter_set == FALSE) { + Log(SEVERE, "Config param 'b_Xinanjiang_shape_parameter' not found in config file. Aborting...\n"); + return BMI_FAILURE; + } + if (is_x_Xinanjiang_shape_parameter_set == FALSE) { + Log(SEVERE, "Config param 'x_Xinanjiang_shape_parameter' not found in config file. Aborting...\n"); + return BMI_FAILURE; + } + if (is_urban_decimal_fraction_set == FALSE) { + Log(SEVERE, "Config param 'urban_decimal_fraction' not found in config file. Aborting...\n"); + return BMI_FAILURE; + } } - /*------------------- surface runoff scheme -AJK----------------------------- */ - if(is_surface_runoff_scheme_set == FALSE) { - model->surface_runoff_scheme = GIUH; + if(model->infiltration_excess_params_struct.surface_water_partitioning_scheme == Schaake){ + model->infiltration_excess_params_struct.Schaake_adjusted_magic_constant_by_soil_type = model->NWM_soil_params.refkdt * model->NWM_soil_params.satdk / 0.000002; + Log(INFO, "Schaake Magic Constant calculated\n"); } model->num_giuh_ordinates = 0; // initialize num_giuh_ordinates @@ -1091,23 +1012,17 @@ int read_init_config_cfe(const char* config_file, cfe_state_struct* model) // Handle GIUH ordinates, bailing if they were not provided if (is_giuh_originates_string_val_set == FALSE) { -#if CFE_DEBUG >= 1 - printf("GIUH ordinate string not set! Aborting...\n"); -#endif + Log(SEVERE, "GIUH ordinate string not set! Aborting...\n"); return BMI_FAILURE; } -#if CFE_DEBUG >= 1 - printf("GIUH ordinates string value found in config ('%s')\n", giuh_originates_string_val); -#endif + Log(INFO, "GIUH ordinates string value found in config ('%s')\n", giuh_originates_string_val); model->num_giuh_ordinates = count_delimited_values(giuh_originates_string_val, ","); -#if CFE_DEBUG >= 1 - printf("Counted number of GIUH ordinates (%d)\n", model->num_giuh_ordinates); -#endif + Log(INFO, "Counted number of GIUH ordinates (%d)\n", model->num_giuh_ordinates); if (model->num_giuh_ordinates < 1) { - printf("Number of GIUH ordinates less than 1 (%d). Aborting... \n", model->num_giuh_ordinates); + Log(SEVERE, "Number of GIUH ordinates less than 1 (%d). Aborting... \n", model->num_giuh_ordinates); return BMI_FAILURE; } @@ -1126,40 +1041,28 @@ int read_init_config_cfe(const char* config_file, cfe_state_struct* model) } else if(model->surface_runoff_scheme == NASH_CASCADE) { if (is_N_nash_surface_set == FALSE) { -#if CFE_DEBUG >= 1 - printf("Config param 'N_nash_surface' not found in config file. Aborting...\n"); -#endif - return BMI_FAILURE; + Log(SEVERE, "Config param 'N_nash_surface' not found in config file. Aborting...\n"); + return BMI_FAILURE; } if (is_K_nash_surface_set == FALSE) { -#if CFE_DEBUG >= 1 - printf("Config param 'K_nash_surface' not found in config file. Aborting...\n"); -#endif - return BMI_FAILURE; + Log(SEVERE, "Config param 'K_nash_surface' not found in config file. Aborting...\n"); + return BMI_FAILURE; } if (is_nsubsteps_nash_surface_set == FALSE) { -#if CFE_DEBUG >= 1 - printf("Config param 'nsubsteps_nash_surface' not found in config file, default value is 10.\n"); -#endif - model->nash_surface_params.nsubsteps = 10; // default value of the number of sub-timesteps + Log(INFO, "Config param 'nsubsteps_nash_surface' not found in config file, default value is 10.\n"); + model->nash_surface_params.nsubsteps = 10; // default value of the number of sub-timesteps } if (is_nash_storage_surface_set == FALSE) { -#if CFE_DEBUG >= 1 - printf("Config param 'nash_storage_surface' not found in config file. Aborting...\n"); -#endif - return BMI_FAILURE; + Log(SEVERE, "Config param 'nash_storage_surface' not found in config file. Aborting...\n"); + return BMI_FAILURE; } if (is_K_infiltration_nash_surface_set == FALSE) { -#if CFE_DEBUG >= 1 - printf("Config param 'Kinf_nash_surface' not found in config file, default value is 0.001 [1/hr] \n"); -#endif - model->nash_surface_params.K_infiltration = 0.001; // used in the runon infiltration + Log(WARNING, "Config param 'Kinf_nash_surface' not found in config file, default value is 0.001 [1/hr] \n"); + model->nash_surface_params.K_infiltration = 0.001; // used in the runon infiltration } if (is_retention_depth_nash_surface_set == FALSE) { -#if CFE_DEBUG >= 1 - printf("Config param 'retention_depth_nash_surface' not found in config file, default value is 1.0 [mm] \n"); -#endif - model->nash_surface_params.retention_depth = 0.001; // usually in the range of 1-5 mm + Log(WARNING, "Config param 'retention_depth_nash_surface' not found in config file, default value is 1.0 [mm] \n"); + model->nash_surface_params.retention_depth = 0.001; // usually in the range of 1-5 mm } @@ -1196,77 +1099,62 @@ int read_init_config_cfe(const char* config_file, cfe_state_struct* model) if(model->infiltration_excess_params_struct.surface_water_partitioning_scheme == Schaake) { model->infiltration_excess_params_struct.Schaake_adjusted_magic_constant_by_soil_type = model->NWM_soil_params.refkdt * model->NWM_soil_params.satdk / 0.000002; - -#if CFE_DEBUG >= 1 - printf("Schaake Magic Constant calculated\n"); -#endif + Log(INFO, "Schaake Magic Constant calculated\n"); } if (is_sft_coupled_set == TRUE && model->infiltration_excess_params_struct.surface_water_partitioning_scheme == Schaake) { if(!is_ice_content_threshold_set) { -#if CFE_DEBUG >= 1 - printf("is_sft_coupled and Schaake scheme are set to TRUE but param 'ice_fraction_threshold' not found in config file\n"); - exit(-9); -#endif + Log(FATAL, "is_sft_coupled and Schaake scheme are set to TRUE but param 'ice_fraction_threshold' not found in config file\n"); + exit(-9); } - + } - + // set sft_coupled flag to false if the parameter is not provided in the config file. model->soil_reservoir.is_sft_coupled = (is_sft_coupled_set == TRUE) ? TRUE : FALSE; - - + /*------------------- Root zone AET development -rlm----------------------------- */ if (is_aet_rootzone_set == TRUE ) { if (is_max_rootzone_layer_set == FALSE) { -#if CFE_DEBUG >= 1 - printf("Config param 'max_rootzone_layer' not found in config file. Aborting...\n"); -#endif + Log(WARNING, "Config param 'max_rootzone_layer' not found in config file\n"); return BMI_FAILURE; } if (is_soil_layer_depths_string_val_set == FALSE) { -#if CFE_DEBUG >= 1 - printf("Soil layer depths string/values not set! Aborting...\n"); -#endif + Log(WARNING, "Soil layer depths string/values not set!\n"); return BMI_FAILURE; } - -#if CFE_DEBUG >=1 - printf("Soil layer depths string values found in config ('%d')\n", is_soil_layer_depths_string_val_set); -#endif - + + Log(INFO, "Soil layer depths string values found in config ('%d')\n", is_soil_layer_depths_string_val_set); + model->soil_reservoir.n_soil_layers = lround(count_delimited_values(soil_layer_depths_string_val, ",")); - -#if CFE_DEBUG >= 1 - printf("Number of soil layer (%d)\n", model->soil_reservoir.n_soil_layers); -#endif - - if (model->soil_reservoir.n_soil_layers < 1) { - printf("Number of soil depths less than 1 (%d). Aborting... \n", model->soil_reservoir.n_soil_layers); - return BMI_FAILURE; - } - - model->soil_reservoir.soil_layer_depths_m = malloc(sizeof(double) * (model->soil_reservoir.n_soil_layers)); + Log(INFO, "n_soil_layers set in bmi_cfe.c: %d\n", model->soil_reservoir.n_soil_layers); + + Log(INFO, "Counted number of soil depths (%d)\n", model->soil_reservoir.n_soil_layers); + + if (model->soil_reservoir.n_soil_layers < 1) + return BMI_FAILURE; + + model->soil_reservoir.soil_layer_depths_m = malloc(sizeof(double) * (model->soil_reservoir.n_soil_layers + 1)); copy = soil_layer_depths_string_val; - - i = 0; + + i = 1; while((value = strsep(©, ",")) != NULL) { model->soil_reservoir.soil_layer_depths_m[i++] = strtod(value, NULL); } - + free(soil_layer_depths_string_val); - + //Check that the last depth read in from the cfe config file (soil_layer_depths) matches the total depth (soil_params.depth) //from the cfe config file. - if(model->NWM_soil_params.D != model->soil_reservoir.soil_layer_depths_m[model->soil_reservoir.n_soil_layers - 1]){ - printf("WARNING: soil_params.depth is not equal to the last soil layer depth in the CFE config file! Aborting...\n"); + if(model->NWM_soil_params.D != model->soil_reservoir.soil_layer_depths_m[model->soil_reservoir.n_soil_layers]){ + Log(SEVERE, "soil_params.depth is not equal to the last soil layer depth in the CFE config file!\n"); return BMI_FAILURE; } - + } - // Calculate the thickness (delta) of each soil layer +// Calculate the thickness (delta) of each soil layer if (is_aet_rootzone_set == TRUE ) { model->soil_reservoir.is_aet_rootzone = TRUE; model->soil_reservoir.delta_soil_layer_depth_m = malloc(sizeof(double) * (model->soil_reservoir.n_soil_layers)); @@ -1275,14 +1163,14 @@ int read_init_config_cfe(const char* config_file, cfe_state_struct* model) for (int i=0; i soil_reservoir.n_soil_layers; i++) { current_depth = model->soil_reservoir.soil_layer_depths_m[i]; if (current_depth <= previous_depth) - printf("WARNING: soil depths may be out of order. One or more soil layer depths is less than or equal to the previous layer. Check CFE config file.\n"); + Log(SEVERE, "Soil depths may be out of order. One or more soil layer depths is less than or equal to the previous layer. Check CFE config file.\n"); model->soil_reservoir.delta_soil_layer_depth_m[i] = current_depth - previous_depth; previous_depth = current_depth; } // Solve for the soil water content at field capacity via Clapp-Hornberger. See "Parameter Estimation for a Conceptual - // Functional Equivalent (CFE) Formulation of the National Water Model" equations 1-3 for detailed description. + // Functional Equivalent (CFE) Formulation of the National Water Model" equations 1-3 for detailed description. double base = (model->NWM_soil_params.alpha_fc * STANDARD_ATMOSPHERIC_PRESSURE_PASCALS)/(WATER_SPECIFIC_WEIGHT * model->NWM_soil_params.satpsi); - double exponent = -1/model->NWM_soil_params.bb; + double exponent = -1/model->NWM_soil_params.bb; model->soil_reservoir.soil_water_content_field_capacity = model->NWM_soil_params.smcmax * pow(base, exponent); } else { @@ -1292,9 +1180,7 @@ int read_init_config_cfe(const char* config_file, cfe_state_struct* model) /*--------------------END OF ROOT ZONE ADJUSTED AET DEVELOPMENT -rlm ------------------------------*/ -#if CFE_DEBUG >= 1 - printf("All CFE config params present\n"); -#endif + Log(INFO, "All CFE config params present\n"); // Now handle the Nash storage array properly (Nash Cascade for subsurface) if (is_nash_storage_subsurface_string_val_set == TRUE) { @@ -1328,9 +1214,7 @@ int read_init_config_cfe(const char* config_file, cfe_state_struct* model) model->nash_storage_subsurface[j] = 0.0; } fclose(fp); -#if CFE_DEBUG >= 1 - printf("Finished function parsing CFE config\n"); -#endif + Log(DEBUG, "Finished function parsing CFE config\n"); return BMI_SUCCESS; } @@ -1338,6 +1222,9 @@ int read_init_config_cfe(const char* config_file, cfe_state_struct* model) static int Initialize (Bmi *self, const char *file) { + // setup the logger + Log(INFO, "In CFE Initialize()\n"); + //FIXME, we can use the input file to help imply "framework" support or "standalone" //an empty init file string indicates things will come from set_value??? //what happens when both occur, that is we have a config file and framewrok @@ -1357,7 +1244,7 @@ static int Initialize (Bmi *self, const char *file) // LKC: removed alpha_fc and S_soil - both added to the soil structure in the model //double S_soil; - + //int is_S_soil_ratio; int config_read_result = read_init_config_cfe(file, cfe_bmi_data_ptr); @@ -1366,7 +1253,7 @@ static int Initialize (Bmi *self, const char *file) // time_step_size is set to 3600 in the "new_bmi_cfe" function. cfe_bmi_data_ptr->timestep_h = cfe_bmi_data_ptr->time_step_size / 3600.0; - + // JG NOTE: this is done in init_soil_reservoir //max_soil_storage = cfe_bmi_data_ptr->NWM_soil_params.D * cfe_bmi_data_ptr->NWM_soil_params.smcmax; @@ -1412,68 +1299,64 @@ static int Initialize (Bmi *self, const char *file) else { - cfe_bmi_data_ptr->is_forcing_from_bmi = FALSE; - - // Figure out the number of lines first (also char count) - int forcing_line_count, max_forcing_line_length; - int count_result = read_file_line_counts_cfe(cfe_bmi_data_ptr->forcing_file, &forcing_line_count, &max_forcing_line_length); - if (count_result == -1) { - printf("Configured forcing file '%s' could not be opened for reading\n", cfe_bmi_data_ptr->forcing_file); - return BMI_FAILURE; - } - if (forcing_line_count == 1) { - printf("Invalid header-only forcing file '%s'\n", cfe_bmi_data_ptr->forcing_file); - return BMI_FAILURE; - } - // Infer the number of time steps: assume a header, so equal to the number of lines minus 1 - cfe_bmi_data_ptr->num_timesteps = forcing_line_count - 1; - -#if CFE_DEBUG > 0 - printf("Counts - Lines: %d | Max Line: %d | Num Time Steps: %d\n", forcing_line_count, max_forcing_line_length, - cfe_bmi_data_ptr->num_timesteps); -#endif - - // Now initialize empty arrays that depend on number of time steps - cfe_bmi_data_ptr->forcing_data_precip_kg_per_m2 = malloc(sizeof(double) * (cfe_bmi_data_ptr->num_timesteps + 1)); - cfe_bmi_data_ptr->forcing_data_time = malloc(sizeof(long) * (cfe_bmi_data_ptr->num_timesteps + 1)); - - // Now open it again to read the forcings - FILE* ffp = fopen(cfe_bmi_data_ptr->forcing_file, "r"); - // Ensure still exists - if (ffp == NULL) { - printf("Forcing file '%s' disappeared!", cfe_bmi_data_ptr->forcing_file); - return BMI_FAILURE; - } - - // Read forcing file and parse forcings - char line_str[max_forcing_line_length + 1]; - long year, month, day, hour, minute; - double dsec; - // First read the header line - if (fgets(line_str, max_forcing_line_length + 1, ffp) == NULL) - return BMI_FAILURE; - - aorc_forcing_data_cfe forcings; - - for (i = 0; i < cfe_bmi_data_ptr->num_timesteps; i++) { - if (fgets(line_str, max_forcing_line_length + 1, ffp) == NULL) // read in a line of AORC data. - return BMI_FAILURE; - parse_aorc_line_cfe(line_str, &year, &month, &day, &hour, &minute, &dsec, &forcings); -#if CFE_DEBUG > 0 - printf("Forcing data: [%s]\n", line_str); - printf("Forcing details - s_time: %ld | precip: %f\n", forcings.time, forcings.precip_kg_per_m2); -#endif - cfe_bmi_data_ptr->forcing_data_precip_kg_per_m2[i] = forcings.precip_kg_per_m2; //* ((double)cfe_bmi_data_ptr->time_step_size); - cfe_bmi_data_ptr->forcing_data_time[i] = forcings.time; - - // TODO: make sure some kind of conversion isn't needed for the rain rate data - // assumed 1000 kg/m3 density of water. This result is mm/h; - //rain_rate[i] = (double) aorc_data.precip_kg_per_m2; - } + cfe_bmi_data_ptr->is_forcing_from_bmi = FALSE; + + // Figure out the number of lines first (also char count) + int forcing_line_count, max_forcing_line_length; + int count_result = read_file_line_counts_cfe(cfe_bmi_data_ptr->forcing_file, &forcing_line_count, &max_forcing_line_length); + if (count_result == -1) { + Log(WARNING, "Configured forcing file '%s' could not be opened for reading\n", cfe_bmi_data_ptr->forcing_file); + return BMI_FAILURE; + } + if (forcing_line_count == 1) { + Log(WARNING, "Invalid header-only forcing file '%s'\n", cfe_bmi_data_ptr->forcing_file); + return BMI_FAILURE; + } + // Infer the number of time steps: assume a header, so equal to the number of lines minus 1 + cfe_bmi_data_ptr->num_timesteps = forcing_line_count - 1; + + Log(INFO, "Counts - Lines: %d | Max Line: %d | Num Time Steps: %d\n", forcing_line_count, max_forcing_line_length, + cfe_bmi_data_ptr->num_timesteps); - cfe_bmi_data_ptr->epoch_start_time = cfe_bmi_data_ptr->forcing_data_time[0]; + // Now initialize empty arrays that depend on number of time steps + cfe_bmi_data_ptr->forcing_data_precip_kg_per_m2 = malloc(sizeof(double) * (cfe_bmi_data_ptr->num_timesteps + 1)); + cfe_bmi_data_ptr->forcing_data_time = malloc(sizeof(long) * (cfe_bmi_data_ptr->num_timesteps + 1)); + + // Now open it again to read the forcings + FILE* ffp = fopen(cfe_bmi_data_ptr->forcing_file, "r"); + // Ensure still exists + if (ffp == NULL) { + Log(WARNING, "Forcing file '%s' disappeared!", cfe_bmi_data_ptr->forcing_file); + return BMI_FAILURE; + } + + // Read forcing file and parse forcings + char line_str[max_forcing_line_length + 1]; + long year, month, day, hour, minute; + double dsec; + // First read the header line + if (fgets(line_str, max_forcing_line_length + 1, ffp) == NULL) + return BMI_FAILURE; + + aorc_forcing_data_cfe forcings; + + for (i = 0; i < cfe_bmi_data_ptr->num_timesteps; i++) { + if (fgets(line_str, max_forcing_line_length + 1, ffp) == NULL) // read in a line of AORC data. + return BMI_FAILURE; + parse_aorc_line_cfe(line_str, &year, &month, &day, &hour, &minute, &dsec, &forcings); + Log(INFO, "Forcing data: [%s]\n", line_str); + Log(INFO, "Forcing details - s_time: %ld | precip: %f\n", forcings.time, forcings.precip_kg_per_m2); + cfe_bmi_data_ptr->forcing_data_precip_kg_per_m2[i] = forcings.precip_kg_per_m2; //* ((double)cfe_bmi_data_ptr->time_step_size); + cfe_bmi_data_ptr->forcing_data_time[i] = forcings.time; + + // TODO: make sure some kind of conversion isn't needed for the rain rate data + // assumed 1000 kg/m3 density of water. This result is mm/h; + //rain_rate[i] = (double) aorc_data.precip_kg_per_m2; + } + + cfe_bmi_data_ptr->epoch_start_time = cfe_bmi_data_ptr->forcing_data_time[0]; } // end if is_forcing_from_bmi - + // divide by 1000 to convert from mm/h to m w/ 1h timestep as per t-shirt_0.99f cfe_bmi_data_ptr->timestep_rainfall_input_m = cfe_bmi_data_ptr->forcing_data_precip_kg_per_m2[0] / 1000; @@ -1512,10 +1395,9 @@ static int Initialize (Bmi *self, const char *file) cfe_bmi_data_ptr->soil_reservoir.smc_profile = malloc(sizeof(double)*cfe_bmi_data_ptr->soil_reservoir.n_soil_layers); else cfe_bmi_data_ptr->soil_reservoir.smc_profile = malloc(sizeof(double)*1); - -#if CFE_DEBUG > 0 - printf("At declaration of smc_profile size, soil_reservoir.n_soil_layers = %i\n", cfe_bmi_data_ptr->soil_reservoir.n_soil_layers); -#endif + + Log(INFO, "At declaration of smc_profile size, soil_reservoir.n_soil_layers = %i\n", cfe_bmi_data_ptr->soil_reservoir.n_soil_layers); + Log(INFO, "Success in CFE BMI Initialization\n"); return BMI_SUCCESS; } @@ -1530,14 +1412,14 @@ static int Update (Bmi *self) // double current_time, end_time; // self->get_current_time(self, ¤t_time); // self->get_end_time(self, &end_time); -// printf("end time: %lf\n", end_time); +// Log(DEBUG, "end time: %lf\n", end_time); // if (current_time >= end_time) { // return BMI_FAILURE; // } - + cfe_state_struct* cfe_ptr = ((cfe_state_struct *) self->data); - // Two modes to get forcing data... 0/FALSE) read from file, 1/TRUE) pass with bmi + // Two modes to get forcing data... 0/FALSE) read from file, 1/TRUE) pass with bmi if (cfe_ptr->is_forcing_from_bmi == TRUE) // BMI sets the precipitation to the aorc structure. // divide by 1000 to convert from mm/h to m w/ 1h timestep as per t-shirt_0.99f @@ -1551,10 +1433,10 @@ static int Update (Bmi *self) cfe_ptr->timestep_rainfall_input_m *= cfe_ptr->time_step_fraction; //Accumulate volume for mass balance cfe_ptr->vol_struct.volin += cfe_ptr->timestep_rainfall_input_m; - + run_cfe(cfe_ptr); - - // Advance the model time + + // Advance the model time cfe_ptr->current_time_step += 1; // compute NWM ponded depth, which is assumed to be the leftover water in the GIUH/NASH queue @@ -1578,18 +1460,22 @@ static int Update_until (Bmi *self, double t) // "the time argument can be a non-integral multiple of time steps" cfe_state_struct* cfe_ptr = ((cfe_state_struct *) self->data); - + double dt; double now; - if(self->get_time_step (self, &dt) == BMI_FAILURE) + if(self->get_time_step (self, &dt) == BMI_FAILURE) { + Log(SEVERE, "get_time_step failure. Aborting...\n"); return BMI_FAILURE; + } - if(self->get_current_time(self, &now) == BMI_FAILURE) + if(self->get_current_time(self, &now) == BMI_FAILURE) { + Log(SEVERE, "get_current_time failure. Aborting...\n"); return BMI_FAILURE; + } { - + int n; double frac; const double n_steps = (t - now) / dt; @@ -1598,8 +1484,8 @@ static int Update_until (Bmi *self, double t) } frac = n_steps - (int)n_steps; if (frac > 0){ - printf("WARNING: CFE trying to update a fraction of a timestep\n"); - + Log(WARNING, "CFE trying to update a fraction of a timestep\n"); + // change timestep to remaining fraction & call update() cfe_ptr->time_step_size = frac * dt; //Record the time step fraction so `Update` can adjust input if needed @@ -1609,7 +1495,7 @@ static int Update_until (Bmi *self, double t) cfe_ptr->time_step_fraction = 1.0; cfe_ptr->time_step_size = dt; } - + } return BMI_SUCCESS; @@ -1626,7 +1512,7 @@ static int Finalize (Bmi *self) free(model->forcing_data_precip_kg_per_m2); if( model->forcing_data_time != NULL ) free(model->forcing_data_time); - + if( model->giuh_ordinates != NULL ) free(model->giuh_ordinates); if( model->nash_storage_subsurface != NULL ) @@ -1635,7 +1521,7 @@ static int Finalize (Bmi *self) free(model->runoff_queue_m_per_timestep); if( model->flux_Qout_m != NULL ) free(model->flux_Qout_m); - + /* xinanjiang_dev: changing name to the more general "direct runoff" if( model->flux_Schaake_output_runoff_m != NULL ) free(model->flux_Schaake_output_runoff_m);*/ @@ -1653,6 +1539,10 @@ static int Finalize (Bmi *self) if( model->flux_perc_m != NULL ) free(model->flux_perc_m); + // clear serialized data + if (model->serialized != NULL) + free(model->serialized); + free(self->data); } @@ -1680,7 +1570,7 @@ static int Get_adjusted_index_for_variable(const char *name) if (strcmp(name, param_var_names[i]) == 0) return i + INPUT_VAR_NAME_COUNT + OUTPUT_VAR_NAME_COUNT; } - + Log(SEVERE, "adjusted index for variable not found. Aborting...\n"); return -1; } @@ -1703,7 +1593,7 @@ static int Get_var_grid(Bmi *self, const char *name, int *grid) } // If we get here, it means the variable name wasn't recognized grid[0] = '\0'; - + Log(SEVERE, "grid variable name wasn't recognized. Aborting...\n"); return BMI_FAILURE; } @@ -1734,8 +1624,24 @@ static int Get_var_type (Bmi *self, const char *name, char * type) } } + // serialization + if (strcmp(name, "serialization_create") == 0) { + strncpy(type, "uint64_t", BMI_MAX_TYPE_NAME); + return BMI_SUCCESS; + } else if (strcmp(name, "serialization_size") == 0) { + strncpy(type, "uint64_t", BMI_MAX_TYPE_NAME); + return BMI_SUCCESS; + } else if (strcmp(name, "serialization_state") == 0) { + strncpy(type, "char", BMI_MAX_TYPE_NAME); + return BMI_SUCCESS; + } else if (strcmp(name, "serialization_free") == 0) { + strncpy(type, "int", BMI_MAX_TYPE_NAME); + return BMI_SUCCESS; + } + // If we get here, it means the variable name wasn't recognized type[0] = '\0'; + Log(SEVERE, "variable type wasn't recognized. Aborting...\n"); return BMI_FAILURE; } @@ -1770,8 +1676,17 @@ static int Get_var_itemsize (Bmi *self, const char *name, int * size) *size = sizeof(long); return BMI_SUCCESS; } + else if (strcmp(type, "uint64_t") == 0) { + *size = sizeof(uint64_t); + return BMI_SUCCESS; + } + else if (strcmp(type, "char") == 0) { + *size = sizeof(char); + return BMI_SUCCESS; + } else { *size = 0; + Log(SEVERE, "variable itemsize wasn't recognized. Aborting...\n"); return BMI_FAILURE; } } @@ -1793,8 +1708,9 @@ static int Get_var_location (Bmi *self, const char *name, char * location) return BMI_SUCCESS; } } - // If we get here, it means the variable name wasn't recognized + // If we get here, it means the variable name location wasn't recognized location[0] = '\0'; + Log(SEVERE, "variable name location wasn't recognized. Aborting...\n"); return BMI_FAILURE; } @@ -1823,6 +1739,20 @@ static int Get_var_units (Bmi *self, const char *name, char * units) static int Get_var_nbytes (Bmi *self, const char *name, int * nbytes) { + // serialization + if (strcmp(name, "serialization_create") == 0 || strcmp(name, "serialization_size") == 0) { + *nbytes = sizeof(uint64_t); + return BMI_SUCCESS; + } else if (strcmp(name, "serialization_state") == 0) { + cfe_state_struct* model = (cfe_state_struct*)self->data; + if (model->serialized == NULL) + return BMI_FAILURE; + *nbytes = model->serialized_length; + return BMI_SUCCESS; + } else if (strcmp(name, "serialization_free") == 0) { + *nbytes = sizeof(int); + return BMI_SUCCESS; + } int item_size; int item_size_result = Get_var_itemsize(self, name, &item_size); if (item_size_result != BMI_SUCCESS) { @@ -1896,12 +1826,12 @@ static int Get_value_ptr (Bmi *self, const char *name, void **dest) if (strcmp (name, "Klf") == 0) { cfe_state_struct *cfe_ptr; cfe_ptr = (cfe_state_struct *) self->data; - // LKC: this seems to be a bug, since in init_soil_reservoir cfe_ptr->soil_reservoir.coeff_secondary = cfe_ptr->K_lf; + // LKC: this seems to be a bug, since in init_soil_reservoir cfe_ptr->soil_reservoir.coeff_secondary = cfe_ptr->K_lf; // Therefore we need to change cfe_ptr->K_lf here *dest = (void*)&cfe_ptr->K_lf; return BMI_SUCCESS; } - + if (strcmp (name, "Kn") == 0) { cfe_state_struct *cfe_ptr; @@ -1915,7 +1845,7 @@ static int Get_value_ptr (Bmi *self, const char *name, void **dest) *dest = (void*)&cfe_ptr->gw_reservoir.coeff_primary; return BMI_SUCCESS; } - + if (strcmp (name, "expon") == 0) { cfe_state_struct *cfe_ptr; cfe_ptr = (cfe_state_struct *) self->data; @@ -1934,42 +1864,42 @@ static int Get_value_ptr (Bmi *self, const char *name, void **dest) cfe_ptr = (cfe_state_struct *) self->data; *dest = (void*)&cfe_ptr->NWM_soil_params.satpsi; return BMI_SUCCESS; - } + } if (strcmp (name, "wltsmc") == 0) { cfe_state_struct *cfe_ptr; cfe_ptr = (cfe_state_struct *) self->data; *dest = (void*)&cfe_ptr->NWM_soil_params.wltsmc; return BMI_SUCCESS; - } + } if (strcmp (name, "alpha_fc") == 0) { cfe_state_struct *cfe_ptr; cfe_ptr = (cfe_state_struct *) self->data; *dest = (void*)&cfe_ptr->NWM_soil_params.alpha_fc; return BMI_SUCCESS; - } + } if (strcmp (name, "refkdt") == 0) { cfe_state_struct *cfe_ptr; cfe_ptr = (cfe_state_struct *) self->data; *dest = (void*)&cfe_ptr->NWM_soil_params.refkdt; return BMI_SUCCESS; - } - if (strcmp (name, "a_Xinanjiang_inflection_point_parameter") == 0) { + } + if (strcmp (name, "a_Xinanjiang_inflection_point_parameter") == 0) { cfe_state_struct *cfe_ptr; cfe_ptr = (cfe_state_struct *) self->data; *dest = (void*)&cfe_ptr->infiltration_excess_params_struct.a_Xinanjiang_inflection_point_parameter; return BMI_SUCCESS; - } + } - if (strcmp (name, "b_Xinanjiang_shape_parameter") == 0) { + if (strcmp (name, "b_Xinanjiang_shape_parameter") == 0) { cfe_state_struct *cfe_ptr; cfe_ptr = (cfe_state_struct *) self->data; *dest = (void*)&cfe_ptr->infiltration_excess_params_struct.b_Xinanjiang_shape_parameter; return BMI_SUCCESS; - } + } - if (strcmp (name, "x_Xinanjiang_shape_parameter") == 0) { + if (strcmp (name, "x_Xinanjiang_shape_parameter") == 0) { cfe_state_struct *cfe_ptr; cfe_ptr = (cfe_state_struct *) self->data; *dest = (void*)&cfe_ptr->infiltration_excess_params_struct.x_Xinanjiang_shape_parameter; @@ -2023,12 +1953,12 @@ static int Get_value_ptr (Bmi *self, const char *name, void **dest) *dest = (void *) ((cfe_state_struct *)(self->data))->flux_perc_m; return BMI_SUCCESS; } - + if (strcmp (name, "Q_OUT") == 0) { *dest = ((cfe_state_struct *)(self->data))->flux_Qout_m; return BMI_SUCCESS; } - + if (strcmp (name, "POTENTIAL_ET") == 0) { cfe_state_struct *cfe_ptr; cfe_ptr = (cfe_state_struct *) self->data; @@ -2042,28 +1972,28 @@ static int Get_value_ptr (Bmi *self, const char *name, void **dest) *dest = (void*)&cfe_ptr-> et_struct.actual_et_m_per_timestep; return BMI_SUCCESS; } - + if (strcmp (name, "GW_STORAGE") == 0) { cfe_state_struct *cfe_ptr; cfe_ptr = (cfe_state_struct *) self->data; *dest = (void*)&cfe_ptr-> gw_reservoir.storage_m; return BMI_SUCCESS; } - + if (strcmp (name, "SOIL_STORAGE") == 0) { cfe_state_struct *cfe_ptr; cfe_ptr = (cfe_state_struct *) self->data; *dest = (void*)&cfe_ptr-> soil_reservoir.storage_m; return BMI_SUCCESS; } - + if (strcmp (name, "SOIL_STORAGE_CHANGE") == 0) { cfe_state_struct *cfe_ptr; cfe_ptr = (cfe_state_struct *) self->data; *dest = (void*)&cfe_ptr->soil_reservoir.storage_change_m; return BMI_SUCCESS; } - + if (strcmp (name, "SURF_RUNOFF_SCHEME") == 0) { cfe_state_struct *cfe_ptr; cfe_ptr = (cfe_state_struct *) self->data; @@ -2086,7 +2016,6 @@ static int Get_value_ptr (Bmi *self, const char *name, void **dest) *dest = (void*)&cfe_ptr-> et_struct.potential_et_m_per_s; return BMI_SUCCESS; } - if (strcmp (name, "atmosphere_water__liquid_equivalent_precipitation_rate") == 0) { cfe_state_struct *cfe_ptr; cfe_ptr = (cfe_state_struct *) self->data; @@ -2100,7 +2029,6 @@ static int Get_value_ptr (Bmi *self, const char *name, void **dest) *dest = (void*)&cfe_ptr->soil_reservoir.ice_fraction_schaake; return BMI_SUCCESS; } - if (strcmp (name, "ice_fraction_xinanjiang") == 0) { cfe_state_struct *cfe_ptr; cfe_ptr = (cfe_state_struct *) self->data; @@ -2115,7 +2043,26 @@ static int Get_value_ptr (Bmi *self, const char *name, void **dest) *dest = (void *) ((cfe_state_struct *)(self->data))->soil_reservoir.smc_profile; return BMI_SUCCESS; } - //------------------------------------------------------------------- +//------------------------------------------------------------------- + + + /***********************************************************/ + /******** SERIALIZATION ******************************/ + /***********************************************************/ + if (strcmp(name, "serialization_size") == 0) { + cfe_state_struct* model = (cfe_state_struct*)self->data; + if (model->serialized == NULL) + return BMI_FAILURE; + *dest = (void*)&model->serialized_length; + return BMI_SUCCESS; + } + if (strcmp(name, "serialization_state") == 0) { + cfe_state_struct* model = (cfe_state_struct*)self->data; + if (model->serialized == NULL) + return BMI_FAILURE; + *dest = model->serialized; + return BMI_SUCCESS; + } return BMI_FAILURE; } @@ -2140,9 +2087,9 @@ static int Get_value_at_indices (Bmi *self, const char *name, void *dest, int *i return BMI_FAILURE; // For now, all variables are non-array scalar values, with only 1 item of type double - + // Thus, there is only ever one value to return (len must be 1) and it must always be from index 0 - if (len > 1 || inds[0] != 0) + if (len > 1 || inds[0] != 0) return BMI_FAILURE; void* ptr; @@ -2156,6 +2103,16 @@ static int Get_value_at_indices (Bmi *self, const char *name, void *dest, int *i static int Get_value (Bmi *self, const char *name, void *dest) { + // special cases for serialization + if (strcmp(name, "serialization_state") == 0) { + cfe_state_struct* model = (cfe_state_struct*)self->data; + if (model->serialized != NULL) { + memcpy(dest, model->serialized, model->serialized_length); + return BMI_SUCCESS; + } + return BMI_FAILURE; + } + // Use nested call to "by index" version // Here, for now at least, we know all the variables are scalar, so @@ -2194,6 +2151,19 @@ static int Set_value_at_indices (Bmi *self, const char *name, int * inds, int le static int Set_value (Bmi *self, const char *name, void *src) { + // serialization + if (strcmp(name, "serialization_state") == 0) { + return load_serialized_cfe(self, (char*)src); + } else if (strcmp(name, "serialization_create") == 0) { + return new_serialized_cfe(self); + } else if (strcmp(name, "serialization_free") == 0) { + return free_serialized_cfe(self); + } + // Avoid using set value, call instead set_value_at_index + // Use nested call to "by index" version + + // Here, for now at least, we know all the variables are scalar, so + int inds[] = {0}; void * dest = NULL; int nbytes = 0; @@ -2298,17 +2268,17 @@ static int Get_state_var_names (Bmi *self, char ** names) int n_state_vars = STATE_VAR_NAME_COUNT; int MAX_NAME_LEN = 512; //int MAX_NAME_LEN = BMI_MAX_VAR_NAME; - + for (int i = 0; i < n_state_vars; i++) { strncpy(names[i], var_info[i].name, MAX_NAME_LEN); - //-------------------------------- + //-------------------------------- // Option to print all the names //-------------------------------- - // if (i==0) printf(" State variable names:"); - // printf(" var name[%d] = %s\n", i, names[i]); + // if (i==0) Log(SEVERE, " State variable names:"); + // Log(SEVERE, " var name[%d] = %s\n", i, names[i]); } - + return BMI_SUCCESS; } @@ -2317,8 +2287,8 @@ static int Get_state_var_types (Bmi *self, char ** types) { //--------------------------------------------------- // Note: This pulls information from the var_info - // structure defined at the top, which helps to - // prevent implementation errors. + // structure defined at the top, which helps to + // prevent implementation errors. //--------------------------------------------------- // This is used for model state serialization, and // returns a string array of all state variable @@ -2330,19 +2300,19 @@ static int Get_state_var_types (Bmi *self, char ** types) return BMI_FAILURE; } - int n_state_vars = STATE_VAR_NAME_COUNT; + int n_state_vars = STATE_VAR_NAME_COUNT; int MAX_NAME_LEN = 512; //int MAX_NAME_LEN = BMI_MAX_VAR_NAME; for (int i = 0; i < n_state_vars; i++) { - strncpy(types[i], var_info[i].type, MAX_NAME_LEN); - //-------------------------------- + strncpy(types[i], var_info[i].type, MAX_NAME_LEN); + //-------------------------------- // Option to print all the types //-------------------------------- - // if (i==0) printf(" State var_types:"); - // printf(" var type[%d] = %s\n", i, types[i]); + // if (i==0) Log(SEVERE, " State var_types:"); + // Log(SEVERE, " var type[%d] = %s\n", i, types[i]); } - + return BMI_SUCCESS; } @@ -2351,8 +2321,8 @@ static int Get_state_var_sizes (Bmi *self, unsigned int size_list[]) { //--------------------------------------------------- // Note: This pulls information from the var_info - // structure defined at the top, which helps to - // prevent implementation errors. + // structure defined at the top, which helps to + // prevent implementation errors. //--------------------------------------------------- // This is used for model state serialization, and // returns a string array of all state variable @@ -2360,8 +2330,8 @@ static int Get_state_var_sizes (Bmi *self, unsigned int size_list[]) // Size is number of array elements (not bytes). // Just number of elements, even for n-dim arrays. //--------------------------------------------------- - if (!self) { - return BMI_FAILURE; + if (!self){ + return BMI_FAILURE; } cfe_state_struct *state; @@ -2389,7 +2359,7 @@ static int Get_state_var_sizes (Bmi *self, unsigned int size_list[]) for (int i = 0; i < n_state_vars; i++) { size_list[i] = var_info[i].size; } - + return BMI_SUCCESS; } @@ -2400,8 +2370,8 @@ static int Get_state_var_ptrs (Bmi *self, void *ptr_list[]) // Return array of pointers to state variables // in same order as defined in state struct. //---------------------------------------------- - if (!self) { - return BMI_FAILURE; + if (!self){ + return BMI_FAILURE; } cfe_state_struct *state; @@ -2417,33 +2387,33 @@ static int Get_state_var_ptrs (Bmi *self, void *ptr_list[]) //-------------------------------------------------- ptr_list[0] = &(state->timestep_rainfall_input_m); - ptr_list[1] = &(state->soil_reservoir_storage_deficit_m); + ptr_list[1] = &(state->soil_reservoir_storage_deficit_m); ptr_list[2] = &(state->infiltration_depth_m); ptr_list[3] = &(state->gw_reservoir_storage_deficit_m); ptr_list[4] = &(state->timestep_h); //-------------------------------- // Vars in soil reservoir struct //-------------------------------- - ptr_list[5] = &(state->soil_reservoir.is_exponential ); - ptr_list[6] = &(state->soil_reservoir.storage_max_m ); - ptr_list[7] = &(state->soil_reservoir.storage_m ); - ptr_list[8] = &(state->soil_reservoir.coeff_primary ); - ptr_list[9] = &(state->soil_reservoir.exponent_primary ); - ptr_list[10] = &(state->soil_reservoir.storage_threshold_primary_m ); - ptr_list[11] = &(state->soil_reservoir.storage_threshold_secondary_m ); - ptr_list[12] = &(state->soil_reservoir.coeff_secondary ); + ptr_list[5] = &(state->soil_reservoir.is_exponential ); + ptr_list[6] = &(state->soil_reservoir.storage_max_m ); + ptr_list[7] = &(state->soil_reservoir.storage_m ); + ptr_list[8] = &(state->soil_reservoir.coeff_primary ); + ptr_list[9] = &(state->soil_reservoir.exponent_primary ); + ptr_list[10] = &(state->soil_reservoir.storage_threshold_primary_m ); + ptr_list[11] = &(state->soil_reservoir.storage_threshold_secondary_m ); + ptr_list[12] = &(state->soil_reservoir.coeff_secondary ); ptr_list[13] = &(state->soil_reservoir.exponent_secondary ); ptr_list[14] = &(state->soil_reservoir.ice_fraction_schaake); ptr_list[15] = &(state->soil_reservoir.ice_fraction_xinanjiang); //------------------------------ // Vars in gw reservoir struct - //------------------------------ + //------------------------------ ptr_list[16] = &(state->gw_reservoir.is_exponential ); ptr_list[17] = &(state->gw_reservoir.storage_max_m ); ptr_list[18] = &(state->gw_reservoir.storage_m ); ptr_list[19] = &(state->gw_reservoir.coeff_primary ); ptr_list[20] = &(state->gw_reservoir.exponent_primary ); - ptr_list[21] = &(state->gw_reservoir.storage_threshold_primary_m ); + ptr_list[21] = &(state->gw_reservoir.storage_threshold_primary_m ); ptr_list[22] = &(state->gw_reservoir.storage_threshold_secondary_m ); ptr_list[23] = &(state->gw_reservoir.coeff_secondary ); ptr_list[24] = &(state->gw_reservoir.exponent_secondary ); @@ -2453,14 +2423,14 @@ static int Get_state_var_ptrs (Bmi *self, void *ptr_list[]) ptr_list[25] = &(state->NWM_soil_params.smcmax ); ptr_list[26] = &(state->NWM_soil_params.wltsmc); ptr_list[27] = &(state->NWM_soil_params.satdk); - ptr_list[28] = &(state->NWM_soil_params.satpsi); + ptr_list[28] = &(state->NWM_soil_params.satpsi); ptr_list[29] = &(state->NWM_soil_params.bb); ptr_list[30] = &(state->NWM_soil_params.slop); ptr_list[31] = &(state->NWM_soil_params.D); ptr_list[32] = &(state->NWM_soil_params.wilting_point_m); //-------------------- // Vars in et_struct - //-------------------- + //-------------------- ptr_list[33] = &(state->et_struct.potential_et_m_per_s ); ptr_list[34] = &(state->et_struct.potential_et_m_per_timestep ); ptr_list[35] = &(state->et_struct.actual_et_m_per_timestep ); @@ -2482,10 +2452,10 @@ static int Get_state_var_ptrs (Bmi *self, void *ptr_list[]) ptr_list[48] = &(state->vol_struct.vol_out_nash ); ptr_list[49] = &(state->vol_struct.vol_in_gw_start ); ptr_list[50] = &(state->vol_struct.vol_soil_start ); - //----------------------------------------- + //----------------------------------------- // More top-level, static allocation vars //----------------------------------------- - ptr_list[51] = &(state->epoch_start_time ); + ptr_list[51] = &(state->epoch_start_time ); ptr_list[52] = &(state->num_timesteps ); ptr_list[53] = &(state->current_time_step ); ptr_list[54] = &(state->time_step_size ); @@ -2512,7 +2482,7 @@ static int Get_state_var_ptrs (Bmi *self, void *ptr_list[]) ptr_list[70] = &(state->aorc.latitude ); ptr_list[71] = &(state->aorc.longitude ); ptr_list[72] = &(state->aorc.time ); - //------------------------------------------ + //------------------------------------------ // More top-level, dynamic allocation vars //---------------------------------------------------- // These vars ARE pointers to different-sized arrays @@ -2531,7 +2501,7 @@ static int Get_state_var_ptrs (Bmi *self, void *ptr_list[]) ptr_list[82] = state->flux_perc_m; ptr_list[83] = state->flux_lat_m; ptr_list[84] = state->flux_Qout_m; - ptr_list[85] = &(state->verbosity); + ptr_list[85] = &(state->verbosity); //--------------------------------------- // infiltration_excess_params_struct vars // xinanjiang or schaake flag [56] @@ -2546,7 +2516,7 @@ static int Get_state_var_ptrs (Bmi *self, void *ptr_list[]) ptr_list[90] = &(state->soil_reservoir.soil_layer_depths_m); ptr_list[91] = &(state->soil_reservoir.max_rootzone_layer); //------------------------------------------------------------- - + return BMI_SUCCESS; } @@ -2587,7 +2557,7 @@ static int Get_state_var_ptrs (Bmi *self, void *ptr_list[]) } int n_state_vars, i; - self->get_state_var_count(self, &n_state_vars); + self->get_state_var_count(self, &n_state_vars); unsigned int sizes[ n_state_vars ]; self->get_state_var_sizes(self, sizes); unsigned int size = sizes[ index ]; @@ -2607,7 +2577,7 @@ static int Get_state_var_ptrs (Bmi *self, void *ptr_list[]) // else if (index == 6){ // for (i=0; ititle[i] = *( (char *)src + i); } } -// else if (index == 7){ +// else if (index == 7){ // for (i=0; isubcat[i] = *( (char *)src + i); } } @@ -2619,9 +2589,9 @@ static int Get_state_var_ptrs (Bmi *self, void *ptr_list[]) else if (index == 1){ state->soil_reservoir_storage_deficit_m = *(double *)src; } else if (index == 2){ - state->infiltration_depth_m = *(double *)src; } + state->infiltration_depth_m = *(double *)src; } else if (index == 3){ - state->gw_reservoir_storage_deficit_m = *(double *)src; } + state->gw_reservoir_storage_deficit_m = *(double *)src; } else if (index == 4){ state->timestep_h = *(double *)src; } //---------------------------------------------------------------- @@ -2634,112 +2604,112 @@ static int Get_state_var_ptrs (Bmi *self, void *ptr_list[]) else if (index == 7){ state->soil_reservoir.storage_m = *(double *)src; } else if (index == 8){ - state->soil_reservoir.coeff_primary = *(double *)src; } + state->soil_reservoir.coeff_primary = *(double *)src; } else if (index == 9){ state->soil_reservoir.exponent_primary = *(double *)src; } else if (index == 10){ - state->soil_reservoir.storage_threshold_primary_m = *(double *)src; } + state->soil_reservoir.storage_threshold_primary_m = *(double *)src; } else if (index == 11){ state->soil_reservoir.storage_threshold_secondary_m = *(double *)src; } else if (index == 12){ - state->soil_reservoir.coeff_secondary = *(double *)src; } + state->soil_reservoir.coeff_secondary = *(double *)src; } else if (index == 13){ state->soil_reservoir.exponent_secondary = *(double *)src; } //---------------------------------------------------------------- // gw_reservoir vars - //---------------------------------------------------------------- + //---------------------------------------------------------------- else if (index == 14){ state->gw_reservoir.is_exponential = *(int *)src; } else if (index == 15){ - state->gw_reservoir.storage_max_m = *(double *)src; } + state->gw_reservoir.storage_max_m = *(double *)src; } else if (index == 16){ - state->gw_reservoir.storage_m = *(double *)src; } + state->gw_reservoir.storage_m = *(double *)src; } else if (index == 17){ - state->gw_reservoir.coeff_primary = *(double *)src; } + state->gw_reservoir.coeff_primary = *(double *)src; } else if (index == 18){ - state->gw_reservoir.exponent_primary = *(double *)src; } + state->gw_reservoir.exponent_primary = *(double *)src; } else if (index == 19){ state->gw_reservoir.storage_threshold_primary_m = *(double *)src; } else if (index == 20){ - state->gw_reservoir.storage_threshold_secondary_m = *(double *)src; } + state->gw_reservoir.storage_threshold_secondary_m = *(double *)src; } else if (index == 21){ - state->gw_reservoir.coeff_secondary = *(double *)src; } + state->gw_reservoir.coeff_secondary = *(double *)src; } else if (index == 22){ state->gw_reservoir.exponent_secondary = *(double *)src; } - //---------------------------------------------------------------- + //---------------------------------------------------------------- // NWM_soil_params vars - //---------------------------------------------------------------- + //---------------------------------------------------------------- else if (index == 23){ - state->NWM_soil_params.smcmax = *(double *)src; } + state->NWM_soil_params.smcmax = *(double *)src; } else if (index == 24){ - state->NWM_soil_params.wltsmc = *(double *)src; } + state->NWM_soil_params.wltsmc = *(double *)src; } else if (index == 25){ - state->NWM_soil_params.satdk = *(double *)src; } + state->NWM_soil_params.satdk = *(double *)src; } else if (index == 26){ - state->NWM_soil_params.satpsi= *(double *)src; } + state->NWM_soil_params.satpsi= *(double *)src; } else if (index == 27){ - state->NWM_soil_params.bb = *(double *)src; } + state->NWM_soil_params.bb = *(double *)src; } else if (index == 28){ - state->NWM_soil_params.slop = *(double *)src; } + state->NWM_soil_params.slop = *(double *)src; } else if (index == 29){ - state->NWM_soil_params.D = *(double *)src; } + state->NWM_soil_params.D = *(double *)src; } else if (index == 30){ state->NWM_soil_params.wilting_point_m = *(double *)src; } //---------------------------------------------------------------- // et_struct vars //---------------------------------------------------------------- else if (index == 31){ - state->et_struct.potential_et_m_per_s = *(double *)src; } + state->et_struct.potential_et_m_per_s = *(double *)src; } else if (index == 32){ - state->et_struct.potential_et_m_per_timestep = *(double *)src; } + state->et_struct.potential_et_m_per_timestep = *(double *)src; } else if (index == 33){ state->et_struct.actual_et_m_per_timestep = *(double *)src; } //---------------------------------------------------------------- // vol_struct vars //---------------------------------------------------------------- else if (index == 34){ - state->vol_struct.vol_runoff = *(double *)src; } + state->vol_struct.vol_runoff = *(double *)src; } else if (index == 35){ - state->vol_struct.vol_infilt = *(double *)src; } + state->vol_struct.vol_infilt = *(double *)src; } else if (index == 36){ state->vol_struct.vol_to_soil = *(double *)src; } else if (index == 37){ - state->vol_struct.vol_to_gw = *(double *)src; } + state->vol_struct.vol_to_gw = *(double *)src; } else if (index == 38){ - state->vol_struct.vol_soil_to_gw = *(double *)src; } + state->vol_struct.vol_soil_to_gw = *(double *)src; } else if (index == 39){ state->vol_struct.vol_soil_to_lat_flow = *(double *)src; } else if (index == 40){ - state->vol_struct.volstart = *(double *)src; } + state->vol_struct.volstart = *(double *)src; } else if (index == 41){ - state->vol_struct.volout = *(double *)src; } + state->vol_struct.volout = *(double *)src; } else if (index == 42){ state->vol_struct.volin = *(double *)src; } else if (index == 43){ - state->vol_struct.vol_from_gw = *(double *)src; } + state->vol_struct.vol_from_gw = *(double *)src; } else if (index == 44){ - state->vol_struct.vol_out_giuh = *(double *)src; } + state->vol_struct.vol_out_giuh = *(double *)src; } else if (index == 45){ state->vol_struct.vol_in_nash = *(double *)src; } else if (index == 46){ - state->vol_struct.vol_out_nash = *(double *)src; } + state->vol_struct.vol_out_nash = *(double *)src; } else if (index == 47){ - state->vol_struct.vol_in_gw_start = *(double *)src; } + state->vol_struct.vol_in_gw_start = *(double *)src; } else if (index == 48){ state->vol_struct.vol_soil_start = *(double *)src; } //---------------------------------------------------------------- // More top-level vars //---------------------------------------------------------------- else if (index == 49){ - state->epoch_start_time = *(long *)src; } + state->epoch_start_time = *(long *)src; } else if (index == 50){ - state->num_timesteps = *(int *)src; } //######## LONG? + state->num_timesteps = *(int *)src; } //######## LONG? else if (index == 51){ state->current_time_step = *(int *)src; } else if (index == 52){ - state->time_step_size = *(int *)src; } + state->time_step_size = *(int *)src; } else if (index == 53){ - state->is_forcing_from_bmi= *(int *)src; } + state->is_forcing_from_bmi= *(int *)src; } else if (index == 54){ // forcing_file is a string memcpy(state->forcing_file, src, size); } @@ -2747,38 +2717,38 @@ static int Get_state_var_ptrs (Bmi *self, void *ptr_list[]) else if (index == 55){ state->infiltration_excess_params_struct.surface_water_partitioning_scheme = *(int *)src; } else if (index == 56){ - state->N_nash = *(int *)src; } + state->N_nash = *(int *)src; } else if (index == 57){ state->K_lf = *(double *)src; } else if (index == 58){ - state->K_nash = *(double *)src; } + state->K_nash = *(double *)src; } else if (index == 59){ - state->num_giuh_ordinates = *(int *)src; } + state->num_giuh_ordinates = *(int *)src; } //---------------------------------------------------------------- // aorc forcing vars //---------------------------------------------------------------- else if (index == 60){ - state->aorc.precip_kg_per_m2 = *(double *)src; } + state->aorc.precip_kg_per_m2 = *(double *)src; } else if (index == 61){ - state->aorc.incoming_longwave_W_per_m2 = *(double *)src; } + state->aorc.incoming_longwave_W_per_m2 = *(double *)src; } else if (index == 62){ state->aorc.incoming_shortwave_W_per_m2 = *(double *)src; } else if (index == 63){ - state->aorc.surface_pressure_Pa = *(double *)src; } + state->aorc.surface_pressure_Pa = *(double *)src; } else if (index == 64){ - state->aorc.specific_humidity_2m_kg_per_kg= *(double *)src; } + state->aorc.specific_humidity_2m_kg_per_kg= *(double *)src; } else if (index == 65){ state->aorc.air_temperature_2m_K = *(double *)src; } else if (index == 66){ - state->aorc.u_wind_speed_10m_m_per_s = *(double *)src; } + state->aorc.u_wind_speed_10m_m_per_s = *(double *)src; } else if (index == 67){ - state->aorc.v_wind_speed_10m_m_per_s = *(double *)src; } + state->aorc.v_wind_speed_10m_m_per_s = *(double *)src; } else if (index == 68){ state->aorc.latitude = *(double *)src; } else if (index == 69){ - state->aorc.longitude= *(double *)src; } + state->aorc.longitude= *(double *)src; } else if (index == 70){ - state->aorc.time = *(long *)src; } + state->aorc.time = *(long *)src; } //---------------------------------------------------------------- // More top-level dynamically-allocated vars // (pointers to scalars or arrays) @@ -2788,12 +2758,12 @@ static int Get_state_var_ptrs (Bmi *self, void *ptr_list[]) // CORRECT: *( ((double *)src) + i) ?? // INCORRECT: *( (double *)(src + i)) // INCORRECT: *( (double *)src + i) ?? - // INCORRECT: *(double *)src + i + // INCORRECT: *(double *)src + i //--------------------------------------------- // Note: state->X is a pointer to an array // We don't need to change that pointer, // just the values in the array. - //--------------------------------------------- + //--------------------------------------------- else if (index == 71){ for (i=0; iforcing_data_precip_kg_per_m2[i] = *( ((double *)src) + i); } } @@ -2805,10 +2775,10 @@ static int Get_state_var_ptrs (Bmi *self, void *ptr_list[]) state->giuh_ordinates[i] = *( ((double *)src) + i); } } else if (index == 74){ for (i=0; inash_storage[i] = *( ((double *)src) + i); } } + state->nash_storage[i] = *( ((double *)src) + i); } } else if (index == 75){ for (i=0; irunoff_queue_m_per_timestep[i] = *( ((double *)src) + i); } } + state->runoff_queue_m_per_timestep[i] = *( ((double *)src) + i); } } else if (index == 76){ for (i=0; iinfiltration_excess_m[i] = *( ((double *)src) + i); } } @@ -2817,24 +2787,24 @@ static int Get_state_var_ptrs (Bmi *self, void *ptr_list[]) state->flux_direct_runoff_m[i] = *( ((double *)src) + i); } } else if (index == 78){ for (i=0; iflux_nash_lateral_runoff_m[i] = *( ((double *)src) + i); } } + state->flux_nash_lateral_runoff_m[i] = *( ((double *)src) + i); } } else if (index == 79){ for (i=0; iflux_from_deep_gw_to_chan_m[i] = *( ((double *)src) + i); } } + state->flux_from_deep_gw_to_chan_m[i] = *( ((double *)src) + i); } } else if (index == 80){ for (i=0; iflux_perc_m[i] = *( ((double *)src) + i); } } + state->flux_perc_m[i] = *( ((double *)src) + i); } } else if (index == 81){ for (i=0; iflux_lat_m[i] = *( ((double *)src) + i); } } + state->flux_lat_m[i] = *( ((double *)src) + i); } } else if (index == 82){ for (i=0; iflux_Qout_m[i] = *( ((double *)src) + i); } } + state->flux_Qout_m[i] = *( ((double *)src) + i); } } else if (index == 83){ // verbosity is not a pointer state->verbosity = *(int *)src; } //-------------------------------------------------------------------------- - // infiltration_excess_params_struct vars (includes xinanjiang AND schaake) + // direct_runoff_params_struc vars (includes xinanjiang AND schaake) //-------------------------------------------------------------------------- else if (index == 84){ state->infiltration_excess_params_struct.Schaake_adjusted_magic_constant_by_soil_type = *(double *)src; } @@ -3023,21 +2993,9 @@ int read_file_line_counts_cfe(const char* file_name, int* line_count, int* max_l cfe_state_struct *new_bmi_cfe(void) { cfe_state_struct *data; - data = (cfe_state_struct *) malloc(sizeof(cfe_state_struct)); + data = (cfe_state_struct *) calloc(1, sizeof(cfe_state_struct)); data->time_step_size = 3600; data->time_step_fraction = 1.0; - data->forcing_data_precip_kg_per_m2 = NULL; - data->forcing_data_time = NULL; - data->giuh_ordinates = NULL; - data->nash_storage_subsurface = NULL; - data->runoff_queue_m_per_timestep = NULL; - data->flux_Qout_m = NULL; - data->infiltration_excess_m = NULL; - data->flux_from_deep_gw_to_chan_m = NULL; - data->flux_direct_runoff_m = NULL; - data->flux_lat_m = NULL; - data->flux_nash_lateral_runoff_m = NULL; - data->flux_perc_m = NULL; return data; } @@ -3077,9 +3035,9 @@ Bmi* register_bmi_cfe(Bmi *model) { model->set_value = Set_value; model->set_value_at_indices = Set_value_at_indices; - model->get_grid_size = Get_grid_size; - model->get_grid_rank = Get_grid_rank; - model->get_grid_type = Get_grid_type; + model->get_grid_size = Get_grid_size; + model->get_grid_rank = Get_grid_rank; + model->get_grid_type = Get_grid_type; model->get_grid_shape = Get_grid_shape; // N/a for grid type scalar model->get_grid_spacing = Get_grid_spacing; // N/a for grid type scalar @@ -3117,7 +3075,7 @@ extern void run_cfe(cfe_state_struct* cfe_ptr){ &cfe_ptr->gw_reservoir_storage_deficit_m, // Set by CFE function after soil_resevroir calc &cfe_ptr->gw_reservoir, // Set in initialize and from config file cfe_ptr->flux_from_deep_gw_to_chan_m, // Set by CFE function after gw_reservoir calc - cfe_ptr->flux_direct_runoff_m, // Set in CFE by convolution_integral or Nash Cascade model + cfe_ptr->flux_direct_runoff_m, // Set in CFE by convolution_integral or Nash Cascade model cfe_ptr->num_giuh_ordinates, // Set by config file with func. count_delimited_values cfe_ptr->giuh_ordinates, // Set by configuration file. cfe_ptr->runoff_queue_m_per_timestep, // Set in initialize @@ -3152,9 +3110,9 @@ extern void init_soil_reservoir(cfe_state_struct* cfe_ptr) // This may need to be changed as follows later, but for now, use the constant value //double Omega = (alpha_fc * cfe->forcing_data_surface_pressure_Pa[0] / WATER_SPECIFIC_WEIGHT) - 0.5; double Omega = ( cfe_ptr->NWM_soil_params.alpha_fc * STANDARD_ATMOSPHERIC_PRESSURE_PASCALS / WATER_SPECIFIC_WEIGHT) - 0.5; - double lower_lim = pow(Omega, (1.0 - 1.0 / cfe_ptr->NWM_soil_params.bb)) / + double lower_lim = pow(Omega, (1.0 - 1.0 / cfe_ptr->NWM_soil_params.bb)) / (1.0 - 1.0 / cfe_ptr->NWM_soil_params.bb); - double upper_lim = pow(Omega + cfe_ptr->NWM_soil_params.D, + double upper_lim = pow(Omega + cfe_ptr->NWM_soil_params.D, (1.0 - 1.0 / cfe_ptr->NWM_soil_params.bb)) / (1.0 - 1.0 / cfe_ptr->NWM_soil_params.bb); @@ -3240,12 +3198,12 @@ extern void initialize_volume_trackers(cfe_state_struct* cfe_ptr) { /**************************************************************************/ /**************************************************************************/ /**************************************************************************/ -extern void print_cfe_flux_header() { - printf("# , hourly , infiltration excess, direct ,lateral, base, total, storage, ice fraction, ice fraction \n"); - printf("Time [h],rainfall [mm], excess [mm],runoff [mm],flow [mm],flow [mm],discharge [mm],storage [mm],schaake [mm],xinan [-]\n"); +extern void print_cfe_flux_header(){ + Log(INFO, "# , hourly , direct, giuh ,lateral, base, total, storage, ice fraction, ice fraction \n"); + Log(INFO, "Time [h],rainfall [mm],runoff [mm],runoff [mm],flow [mm],flow [mm],discharge [mm],storage [mm],schaake [mm],xinan [-]\n"); } extern void print_cfe_flux_at_timestep(cfe_state_struct* cfe_ptr) { - printf("%d, %lf, %lf, %lf, %lf, %lf, %lf, %lf, %lf, %lf\n", + Log(INFO, "%d, %lf, %lf, %lf, %lf, %lf, %lf, %lf, %lf, %lf\n", cfe_ptr->current_time_step, cfe_ptr->timestep_rainfall_input_m*1000.0, *cfe_ptr->infiltration_excess_m*1000.0, @@ -3268,15 +3226,15 @@ extern void mass_balance_check(cfe_state_struct* cfe_ptr){ double vol_surface_end = 0.0; // volume in the giuh or nash array at the end (on the surface) double vol_nash_subsurface_end = 0.0; // volume at the end in the nash cascade for subsurface double vol_soil_end; - + // the GIUH queue might have water in it at the end of the simulation, so sum it up. if (cfe_ptr->surface_runoff_scheme == GIUH) { for(i=0;inum_giuh_ordinates;i++) - vol_surface_end += cfe_ptr->runoff_queue_m_per_timestep[i]; + vol_surface_end += cfe_ptr->runoff_queue_m_per_timestep[i]; } else if (cfe_ptr->surface_runoff_scheme == NASH_CASCADE) { for(i=0;inash_surface_params.N_nash;i++) - vol_surface_end += cfe_ptr->nash_surface_params.nash_storage[i]; + vol_surface_end += cfe_ptr->nash_surface_params.nash_storage[i]; } cfe_ptr->vol_struct.vol_end_surface = vol_surface_end; // update the vol in the mass balance struct @@ -3284,7 +3242,7 @@ extern void mass_balance_check(cfe_state_struct* cfe_ptr){ for(i=0;iN_nash_subsurface;i++) vol_nash_subsurface_end += cfe_ptr->nash_storage_subsurface[i]; vol_soil_end=cfe_ptr->soil_reservoir.storage_m; - + double global_residual; /* xinanjiang_dev @@ -3295,83 +3253,87 @@ extern void mass_balance_check(cfe_state_struct* cfe_ptr){ double soil_residual; double nash_residual; double gw_residual; + global_residual = cfe_ptr->vol_struct.volstart + cfe_ptr->vol_struct.volin - cfe_ptr->vol_struct.volout - volend - cfe_ptr->vol_struct.vol_end_surface; - printf("========================= Simulation Summary ========================= \n"); - printf("********************* GLOBAL MASS BALANCE ********************* \n"); - printf(" Volume initial = %8.4lf m\n",cfe_ptr->vol_struct.volstart); - printf(" Volume input = %8.4lf m\n",cfe_ptr->vol_struct.volin); - printf(" Volume output = %8.4lf m\n",cfe_ptr->vol_struct.volout); - printf(" Final volume = %8.4lf m\n",volend); - printf(" Global residual = %6.4e m\n",global_residual); + Log(INFO, "========================= Simulation Summary ========================= \n"); + Log(INFO, "********************* GLOBAL MASS BALANCE ********************* \n"); + Log(INFO, " Volume initial = %8.4lf m\n",cfe_ptr->vol_struct.volstart); + Log(INFO, " Volume input = %8.4lf m\n",cfe_ptr->vol_struct.volin); + Log(INFO, " Volume output = %8.4lf m\n",cfe_ptr->vol_struct.volout); + Log(INFO, " Final volume = %8.4lf m\n",volend); + Log(INFO, " Global residual = %6.4e m\n",global_residual); if(cfe_ptr->vol_struct.volin>0.0) - printf(" Global percent error = %6.4e percent of inputs\n",global_residual/cfe_ptr->vol_struct.volin*100.0); + Log(INFO, " Global percent error = %6.4e percent of inputs\n",global_residual/cfe_ptr->vol_struct.volin*100.0); else - printf(" Global pct. err: %6.4e percent of initial\n",global_residual/cfe_ptr->vol_struct.volstart*100.0); + Log(INFO, " Global pct. err: %6.4e percent of initial\n",global_residual/cfe_ptr->vol_struct.volstart*100.0); if(!is_fabs_less_than_epsilon(global_residual,1.0e-12)) - printf("WARNING: GLOBAL MASS BALANCE CHECK FAILED\n"); + Log(WARNING, "GLOBAL MASS BALANCE CHECK FAILED\n"); direct_residual = cfe_ptr->vol_struct.volin - cfe_ptr->vol_struct.vol_runoff - cfe_ptr->vol_struct.vol_infilt - cfe_ptr->vol_struct.vol_et_from_rain; - printf("******************* PRECIPITATION MASS BALANCE *****************\n"); - printf(" Volume input = %8.4lf m\n",cfe_ptr->vol_struct.volin); - printf(" Infiltration Excess = %8.4lf m\n",cfe_ptr->vol_struct.vol_runoff); - printf(" Infiltration = %8.4lf m\n",cfe_ptr->vol_struct.vol_infilt); - printf(" Vol_et_from_rain = %8.4lf m\n",cfe_ptr->vol_struct.vol_et_from_rain); - printf(" Precip residual = %6.4e m\n",direct_residual); // should equal 0.0 - if(!is_fabs_less_than_epsilon(direct_residual,1.0e-12)) - printf("WARNING: DIRECT RUNOFF PARTITIONING MASS BALANCE CHECK FAILED\n"); - + Log(INFO, "******************* PRECIPITATION MASS BALANCE *****************\n"); + Log(INFO, " Volume input = %8.4lf m\n",cfe_ptr->vol_struct.volin); + Log(INFO, " Infiltration Excess = %8.4lf m\n",cfe_ptr->vol_struct.vol_runoff); + Log(INFO, " Infiltration = %8.4lf m\n",cfe_ptr->vol_struct.vol_infilt); + Log(INFO, " Vol_et_from_rain = %8.4lf m\n",cfe_ptr->vol_struct.vol_et_from_rain); + Log(INFO, " Precip residual = %6.4e m\n",direct_residual); // should equal 0.0 + if (!is_fabs_less_than_epsilon(direct_residual,1.0e-12)) + Log(SEVERE, "DIRECT RUNOFF PARTITIONING MASS BALANCE CHECK FAILEDS\n"); + giuh_residual = cfe_ptr->vol_struct.vol_runoff - cfe_ptr->vol_struct.vol_out_surface - cfe_ptr->vol_struct.vol_end_surface - cfe_ptr->vol_struct.vol_runon_infilt; - printf("********************* SURFACE MASS BALANCE *********************\n"); - printf(" Volume into surface = %8.4lf m\n",cfe_ptr->vol_struct.vol_runoff); - printf(" Volume out surface = %8.4lf m\n",cfe_ptr->vol_struct.vol_out_surface); - printf(" Volume end surface = %8.4lf m\n",cfe_ptr->vol_struct.vol_end_surface); - printf(" Runon infiltration = %8.4lf m\n",cfe_ptr->vol_struct.vol_runon_infilt); - printf(" Vol_et_from_retention_depth = %8.4lf m\n",cfe_ptr->vol_struct.vol_et_from_retention_depth); - printf(" Surface residual = %6.4e m\n", giuh_residual); // should equal zero + Log(INFO, "********************* SURFACE MASS BALANCE *********************\n"); + Log(INFO, " Volume into surface = %8.4lf m\n",cfe_ptr->vol_struct.vol_runoff); + Log(INFO, " Volume out surface = %8.4lf m\n",cfe_ptr->vol_struct.vol_out_surface); + Log(INFO, " Volume end surface = %8.4lf m\n",cfe_ptr->vol_struct.vol_end_surface); + Log(INFO, " Runon infiltration = %8.4lf m\n",cfe_ptr->vol_struct.vol_runon_infilt); + Log(INFO, " Vol_et_from_retention_depth = %8.4lf m\n",cfe_ptr->vol_struct.vol_et_from_retention_depth); + Log(INFO, " Surface residual = %6.4e m\n", giuh_residual); // should equal zero if(!is_fabs_less_than_epsilon(giuh_residual,1.0e-12)) - printf("WARNING: GIUH MASS BALANCE CHECK FAILED\n"); - + Log(SEVERE, "GIUH MASS BALANCE CHECK FAILED\n"); soil_residual = cfe_ptr->vol_struct.vol_soil_start + cfe_ptr->vol_struct.vol_infilt + cfe_ptr->vol_struct.vol_runon_infilt - cfe_ptr->vol_struct.vol_soil_to_lat_flow - vol_soil_end - cfe_ptr->vol_struct.vol_to_gw - cfe_ptr->vol_struct.vol_et_from_soil; - - printf("*********** SOIL WATER CONCEPTUAL RESERVOIR MASS BALANCE *******\n"); - printf(" Initial soil vol = %8.4lf m\n",cfe_ptr->vol_struct.vol_soil_start); - printf(" Volume into soil = %8.4lf m\n",cfe_ptr->vol_struct.vol_infilt); - printf(" Volume soil2latflow = %8.4lf m\n",cfe_ptr->vol_struct.vol_soil_to_lat_flow); - printf(" Volume soil to GW = %8.4lf m\n",cfe_ptr->vol_struct.vol_soil_to_gw); - printf(" Final volume soil = %8.4lf m\n",vol_soil_end); - printf(" Volume et from soil = %8.4lf m\n",cfe_ptr->vol_struct.vol_et_from_soil); - printf(" Volume soil residual = %6.4e m\n",soil_residual); + + Log(INFO, " SOIL WATER CONCEPTUAL RESERVOIR MASS BALANCE\n"); + Log(INFO, " init soil vol: %8.4lf m\n",cfe_ptr->vol_struct.vol_soil_start); + + + Log(INFO, "*********** SOIL WATER CONCEPTUAL RESERVOIR MASS BALANCE *******\n"); + Log(INFO, " Initial soil vol = %8.4lf m\n",cfe_ptr->vol_struct.vol_soil_start); + Log(INFO, " Volume into soil = %8.4lf m\n",cfe_ptr->vol_struct.vol_infilt); + Log(INFO, " Volume soil2latflow = %8.4lf m\n",cfe_ptr->vol_struct.vol_soil_to_lat_flow); + Log(INFO, " Volume soil to GW = %8.4lf m\n",cfe_ptr->vol_struct.vol_soil_to_gw); + Log(INFO, " Final volume soil = %8.4lf m\n",vol_soil_end); + Log(INFO, " Volume et from soil = %8.4lf m\n",cfe_ptr->vol_struct.vol_et_from_soil); + Log(INFO, " Volume soil residual = %6.4e m\n",soil_residual); if(!is_fabs_less_than_epsilon(soil_residual,1.0e-12)) - printf("WARNING: SOIL CONCEPTUAL RESERVOIR MASS BALANCE CHECK FAILED\n"); + Log(WARNING, "SOIL CONCEPTUAL RESERVOIR MASS BALANCE CHECK FAILED\n"); nash_residual = cfe_ptr->vol_struct.vol_in_nash - cfe_ptr->vol_struct.vol_out_nash - vol_nash_subsurface_end; - printf("************* NASH CASCADE (SUBSURFACE) MASS BALANCE ***********\n"); - printf(" Volume to Nash = %8.4lf m\n",cfe_ptr->vol_struct.vol_in_nash); - printf(" Volume from Nash = %8.4lf m\n",cfe_ptr->vol_struct.vol_out_nash); - printf(" Final vol. Nash = %8.4lf m\n",vol_nash_subsurface_end); - printf(" Nash casc resid. = %6.4e m\n",nash_residual); + Log(INFO, "************* NASH CASCADE (SUBSURFACE) MASS BALANCE ***********\n"); + Log(INFO, " Volume to Nash = %8.4lf m\n",cfe_ptr->vol_struct.vol_in_nash); + Log(INFO, " Volume from Nash = %8.4lf m\n",cfe_ptr->vol_struct.vol_out_nash); + Log(INFO, " Final vol. Nash = %8.4lf m\n",vol_nash_subsurface_end); + Log(INFO, " Nash casc resid. = %6.4e m\n",nash_residual); if(!is_fabs_less_than_epsilon(nash_residual,1.0e-12)) - printf("WARNING: NASH CASCADE CONCEPTUAL RESERVOIR MASS BALANCE CHECK FAILED\n"); - - + Log(SEVERE, "NASH CASCADE CONCEPTUAL RESERVOIR MASS BALANCE CHECK FAILED\n"); + + gw_residual = cfe_ptr->vol_struct.vol_in_gw_start + cfe_ptr->vol_struct.vol_to_gw - cfe_ptr->vol_struct.vol_from_gw - vol_in_gw_end; - printf("********** GROUNDWATER CONCEPTUAL RESERVOIR MASS BALANCE *******\n"); - printf(" Initial GW storage = %8.4lf m\n",cfe_ptr->vol_struct.vol_in_gw_start); - printf(" Volume to GW = %8.4lf m\n",cfe_ptr->vol_struct.vol_to_gw); - printf(" Volume from GW = %8.4lf m\n",cfe_ptr->vol_struct.vol_from_gw); - printf(" Final GW storage = %8.4lf m\n",vol_in_gw_end); - printf(" GW residual = %6.4e m\n",gw_residual); + Log(INFO, "********** GROUNDWATER CONCEPTUAL RESERVOIR MASS BALANCE *******\n"); + Log(INFO, " Initial GW storage = %8.4lf m\n",cfe_ptr->vol_struct.vol_in_gw_start); + Log(INFO, " Volume to GW = %8.4lf m\n",cfe_ptr->vol_struct.vol_to_gw); + Log(INFO, " Volume from GW = %8.4lf m\n",cfe_ptr->vol_struct.vol_from_gw); + Log(INFO, " Final GW storage = %8.4lf m\n",vol_in_gw_end); + Log(INFO, " GW residual = %6.4e m\n",gw_residual); if(!is_fabs_less_than_epsilon(gw_residual,1.0e-12)) - fprintf(stderr,"WARNING: GROUNDWATER CONCEPTUAL RESERVOIR MASS BALANCE CHECK FAILED\n"); + Log(WARNING, "GROUNDWATER CONCEPTUAL RESERVOIR MASS BALANCE CHECK FAILED\n"); } /**************************************************************************/ @@ -3392,59 +3354,59 @@ void parse_aorc_line_cfe(char *theString,long *year,long *month, long *day,long int start,end; int wordlen; char theWord[150]; - + //len=strlen(theString); - + start=0; /* begin at the beginning of theString */ get_word_cfe(theString,&start,&end,theWord,&wordlen); *year=atol(theWord); - + get_word_cfe(theString,&start,&end,theWord,&wordlen); *month=atol(theWord); - + get_word_cfe(theString,&start,&end,theWord,&wordlen); *day=atol(theWord); - + get_word_cfe(theString,&start,&end,theWord,&wordlen); *hour=atol(theWord); - + get_word_cfe(theString,&start,&end,theWord,&wordlen); *minute=atol(theWord); - + get_word_cfe(theString,&start,&end,theWord,&wordlen); *second=(double)atof(theWord); - + get_word_cfe(theString,&start,&end,theWord,&wordlen); aorc->precip_kg_per_m2=atof(theWord); //printf("%s, %s, %lf, %lf\n", theString, theWord, (double)atof(theWord), aorc->precip_kg_per_m2); - + get_word_cfe(theString,&start,&end,theWord,&wordlen); - aorc->incoming_longwave_W_per_m2=(double)atof(theWord); - + aorc->incoming_longwave_W_per_m2=(double)atof(theWord); + get_word_cfe(theString,&start,&end,theWord,&wordlen); - aorc->incoming_shortwave_W_per_m2=(double)atof(theWord); - + aorc->incoming_shortwave_W_per_m2=(double)atof(theWord); + get_word_cfe(theString,&start,&end,theWord,&wordlen); - aorc->surface_pressure_Pa=(double)atof(theWord); - + aorc->surface_pressure_Pa=(double)atof(theWord); + get_word_cfe(theString,&start,&end,theWord,&wordlen); aorc->specific_humidity_2m_kg_per_kg=(double)atof(theWord); - + get_word_cfe(theString,&start,&end,theWord,&wordlen); - aorc->air_temperature_2m_K=(double)atof(theWord); - + aorc->air_temperature_2m_K=(double)atof(theWord); + get_word_cfe(theString,&start,&end,theWord,&wordlen); - aorc->u_wind_speed_10m_m_per_s=(double)atof(theWord); - + aorc->u_wind_speed_10m_m_per_s=(double)atof(theWord); + get_word_cfe(theString,&start,&end,theWord,&wordlen); - aorc->v_wind_speed_10m_m_per_s=(double)atof(theWord); + aorc->v_wind_speed_10m_m_per_s=(double)atof(theWord); - //---------------------------------------------- + //---------------------------------------------- // Is this needed? It has not been set. (SDP) //---------------------------------------------- aorc->time = -9999.0; //###################################################### - + return; } @@ -3493,7 +3455,7 @@ void itwo_alloc_cfe(int ***array, int rows, int cols) { int i, frows, fcols; if ((rows == 0) || (cols == 0)) { - printf("Error: Attempting to allocate array of size 0\n"); + Log(SEVERE, "Attempting to allocate array of size 0\n"); } frows = rows + 1; /* added one for FORTRAN numbering */ @@ -3518,7 +3480,7 @@ void dtwo_alloc_cfe(double ***array, int rows, int cols) { int i, frows, fcols; if ((rows == 0) || (cols == 0)) { - printf("Error: Attempting to allocate array of size 0\n"); + Log(SEVERE, "Error: Attempting to allocate array of size 0\n"); } frows = rows + 1; /* added one for FORTRAN numbering */ @@ -3544,7 +3506,7 @@ void d_alloc_cfe(double **var, int size) { *var = (double *) malloc(size * sizeof(double)); if (*var == NULL) { - printf("Problem allocating memory for array in d_alloc.\n"); + Log(SEVERE, "Problem allocating memory for array in d_alloc.\n"); return; } else memset(*var, 0, size * sizeof(double)); @@ -3556,7 +3518,7 @@ void i_alloc_cfe(int **var, int size) { *var = (int *) malloc(size * sizeof(int)); if (*var == NULL) { - printf("Problem allocating memory in i_alloc\n"); + Log(SEVERE, "Problem allocating memory in i_alloc\n"); return; } else memset(*var, 0, size * sizeof(int)); @@ -3621,7 +3583,7 @@ double greg_2_jul_cfe(long year, ** Copied from Algorithm 199 in Collected algorithms of the CACM ** Author: Robert G. Tantzen, Translator: Nat Howard ** -** Modified by FLO 4/99 to account for nagging round off error +** Modified by FLO 4/99 to account for nagging round off error ** */ void calc_date_cfe(double jd, long *y, long *m, long *d, long *h, long *mi, @@ -3673,3 +3635,4 @@ int dayofweek(double j) { j += 0.5; return (int) (j + 1) % 7; } + diff --git a/src/bmi_serialization.cxx b/src/bmi_serialization.cxx new file mode 100644 index 00000000..a2dcabeb --- /dev/null +++ b/src/bmi_serialization.cxx @@ -0,0 +1,162 @@ +extern "C" { +#include "bmi_serialization.h" +#include "bmi_cfe.h" +#include "logger.h" +} +#include "vecbuf.hpp" + +#include + +#include +#include +#include + + +class CfeSerializer { + public: + CfeSerializer(Bmi* bmi); + cfe_state_struct* data; + private: + friend class boost::serialization::access; + template + void serialize(Archive& ar, const unsigned int version); +}; + + +CfeSerializer::CfeSerializer(Bmi* bmi) { + void* data = bmi->data; + this->data = (cfe_state_struct*)data; +} + + +template +void CfeSerializer::serialize(Archive& ar, const unsigned int version) { + using boost::serialization::make_array; + cfe_state_struct* state = this->data; + + // base + ar & state->current_time_step; // update + ar & state->timestep_rainfall_input_m; // update + ar & state->nwm_ponded_depth_m; // update + ar & state->soil_reservoir_storage_deficit_m; // cfe + ar & *state->infiltration_excess_m; // cfe; stored as pointer on struct + ar & state->infiltration_depth_m; // cfe + ar & *state->flux_perc_m; // cfe; stored as pointer on struct + ar & *state->flux_lat_m; // cfe; stored as pointer on struct + ar & state->gw_reservoir_storage_deficit_m; // cfe + ar & *state->flux_from_deep_gw_to_chan_m; // cfe; stored as pointer on struct + ar & *state->flux_direct_runoff_m; // cfe; stored as pointer on struct + ar & *state->flux_nash_lateral_runoff_m; // cfe, nash_cascade_subsurface; stored as pointer on struct + ar & *state->flux_Qout_m; // cfe; stored as pointer on struct + ar & make_array(state->nash_storage_subsurface, // nash_cascade_subsurface + state->N_nash_subsurface); // config + + // conceptual_reservoir + ar & state->soil_reservoir.storage_m; // cfe, et_from_soil + ar & state->soil_reservoir.storage_change_m; // cfe + ar & make_array(state->soil_reservoir.smc_profile, // et_from_soil + state->soil_reservoir.n_soil_layers); // config + ar & state->soil_reservoir.ice_fraction_schaake; // cfe + ar & state->soil_reservoir.ice_fraction_xinanjiang; // cfe + ar & state->gw_reservoir.storage_m; // cfe + + // nash_cascade_parameters + if (state->nash_surface_params.nash_storage != NULL) { // will be null or malloc'd based on config + ar & make_array(state->nash_surface_params.nash_storage, // et_from_retention_depth, nash_cascade_surface + state->nash_surface_params.N_nash); // config + } + + // evapotransporation_struct + ar & state->et_struct.potential_et_m_per_timestep; // cfe + ar & state->et_struct.reduced_potential_et_m_per_timestep; // cfe, et_from_soil, et_from_rainfall + ar & state->et_struct.actual_et_from_rain_m_per_timestep; // cfe + ar & state->et_struct.actual_et_from_soil_m_per_timestep; // cfe, et_from_soil + ar & state->et_struct.actual_et_m_per_timestep; // cfe + ar & state->et_struct.potential_et_m_per_s; // cfe + + // massball, all set in cfe + ar & state->vol_struct.volin; // set in update + ar & state->vol_struct.vol_et_from_rain; + ar & state->vol_struct.vol_et_to_atm; + ar & state->vol_struct.volout; + ar & state->vol_struct.vol_et_from_retention_depth; + ar & state->vol_struct.vol_out_surface; + ar & state->vol_struct.vol_et_from_soil; + ar & state->vol_struct.vol_runoff; + ar & state->vol_struct.vol_infilt; + ar & state->vol_struct.vol_to_soil; + ar & state->vol_struct.vol_to_gw; + ar & state->vol_struct.vol_soil_to_gw; + ar & state->vol_struct.vol_soil_to_lat_flow; + ar & state->vol_struct.vol_from_gw; + ar & state->vol_struct.vol_out_surface; + ar & state->vol_struct.vol_in_nash; + ar & state->vol_struct.vol_out_nash; + + if (state->surface_runoff_scheme == GIUH) { // config + ar & make_array(state->runoff_queue_m_per_timestep, // giuh_convolution_integral + state->num_giuh_ordinates + 1); // config + } + /// state->nash_surface_params.runon_infiltration not used for input or GetValue + // else if (state->surface_runoff_scheme == NASH_CASCADE) { + // ar & state->nash_surface_params.runon_infiltration; + // } +} + + +extern "C" { + +int free_serialized_cfe(Bmi* bmi) { + cfe_state_struct* model = (cfe_state_struct*)bmi->data; + if (model->serialized != NULL) { + free(model->serialized); + model->serialized = NULL; + } + model->serialized_length = 0; + return BMI_SUCCESS; +} + +int load_serialized_cfe(Bmi* bmi, const char* data) { + CfeSerializer serializer(bmi); + std::istringstream stream(data); + boost::archive::binary_iarchive archive(stream); + try { + archive >> serializer; + return BMI_SUCCESS; + } catch (const std::exception &e) { + Log(LogLevel::SEVERE, "Deserializing CFE encountered an error: %s", e.what()); + return BMI_FAILURE; + } +} + +int new_serialized_cfe(Bmi* bmi) { + CfeSerializer serializer(bmi); + vecbuf stream; + boost::archive::binary_oarchive archive(stream); + try { + archive << serializer; + } catch (const std::exception &e) { + Log(LogLevel::SEVERE, "Serializing CFE encountered an error: %s", e.what()); + return BMI_FAILURE; + } + // copy serialized data into cfe data + cfe_state_struct* model = (cfe_state_struct*)bmi->data; + // clear previous data if it exists + if (model->serialized != NULL) { + free(model->serialized); + } + // set size and allocate memory + model->serialized_length = stream.size(); + model->serialized = (char*)malloc(model->serialized_length); + // make sure memory could be allocated + if (model->serialized == NULL) { + Log(LogLevel::SEVERE, "CFE failed to allocate memory for serializing."); + model->serialized_length = 0; + return BMI_FAILURE; + } + // copy stream data to new allocation + memcpy(model->serialized, stream.data(), model->serialized_length); + return BMI_SUCCESS; +} + +} diff --git a/src/cfe.c b/src/cfe.c index 7e0fa015..62eb014b 100644 --- a/src/cfe.c +++ b/src/cfe.c @@ -1,8 +1,12 @@ #include "cfe.h" +#include "logger.h" +#include #define max(a,b) ({ __typeof__ (a) _a = (a); __typeof__ (b) _b = (b); _a > _b ? _a : _b; }) #define min(a,b) ({ __typeof__ (a) _a = (a); __typeof__ (b) _b = (b); _a < _b ? _a : _b; }) +static bool groundwaterFluxIssueLogged = false; + // CFE STATE SPACE FUNCTION // ####################################################################### // Adapted version of Conceptual Functional Equivalent model re-written in state-space form July, 2021 //#################################################################################################### @@ -49,12 +53,12 @@ extern void cfe( double flux_lat_m = *flux_lat_m_ptr; // water moved from soil reservoir to lateral flow Nash cascad this timestep [m] double gw_reservoir_storage_deficit_m = *gw_reservoir_storage_deficit_m_ptr; // deficit in gw reservoir storage [m] double flux_from_deep_gw_to_chan_m = *flux_from_deep_gw_to_chan_m_ptr; // water moved from gw reservoir to catchment outlet nexus this timestep [m] - double direct_runoff_m = *direct_runoff_m_ptr; // water leaving GIUH or NASH cascade reservoir to outlet this timestep [m] + double direct_runoff_m = *direct_runoff_m_ptr; // water leaving GIUH or NASH cascade reservoir to outlet this timestep [m] double nash_lateral_runoff_m = *nash_lateral_runoff_m_ptr; // water leaving lateral subsurface flow Nash cascade this timestep [m] double Qout_m = *Qout_m_ptr; // the total runoff this timestep (Surface+Nash+GW) [m] // LOCAL VARIABLES, the values of which are not important to describe the model state. They are like notes on scrap paper. - + double diff=0.0; double primary_flux=0.0; // pointers to these variables passed to conceptual nonlinear reservoir which has two outlets, primary & secondary double secondary_flux=0.0; // pointers to these variables passed to conceptual nonlinear reservoir which has two outlets, primary & secondary @@ -64,14 +68,14 @@ extern void cfe( // store current soil storage_m in a temp. variable double storage_temp_m = soil_reservoir_struct->storage_m; - + //################################################## // partition rainfall using Schaake function //################################################## evap_struct->potential_et_m_per_timestep = evap_struct->potential_et_m_per_s * time_step_size; evap_struct->reduced_potential_et_m_per_timestep = evap_struct->potential_et_m_per_s * time_step_size; - + evap_struct->actual_et_from_rain_m_per_timestep = 0; if (timestep_rainfall_input_m > 0) { et_from_rainfall(×tep_rainfall_input_m,evap_struct); @@ -114,16 +118,16 @@ extern void cfe( // LKC: This needs to be calcualted here after et_from_soil since soil_reservoir_struct->storage_m changes soil_reservoir_storage_deficit_m=(NWM_soil_params_struct.smcmax*NWM_soil_params_struct.D-soil_reservoir_struct->storage_m); - + // NEW FLO - if(0.0 < timestep_rainfall_input_m) + if(0.0 < timestep_rainfall_input_m) { if (infiltration_excess_params_struct.surface_water_partitioning_scheme == Schaake) { // to ensure that ice_fraction_schaake is set to 0.0 for uncoupled SFT if(!soil_reservoir_struct->is_sft_coupled) soil_reservoir_struct->ice_fraction_schaake = 0.0; - + Schaake_partitioning_scheme(timestep_h, soil_reservoir_struct->storage_threshold_primary_m, infiltration_excess_params_struct.Schaake_adjusted_magic_constant_by_soil_type, soil_reservoir_storage_deficit_m, timestep_rainfall_input_m, NWM_soil_params_struct.smcmax, @@ -136,7 +140,7 @@ extern void cfe( // to ensure that ice_fraction_xinan is set to 0.0 for uncoupled SFT if(!soil_reservoir_struct->is_sft_coupled) soil_reservoir_struct->ice_fraction_xinanjiang = 0.0; - + Xinanjiang_partitioning_scheme(timestep_rainfall_input_m, soil_reservoir_struct->storage_threshold_primary_m, soil_reservoir_struct->storage_max_m, soil_reservoir_struct->storage_m, &infiltration_excess_params_struct, &infiltration_excess_m, &infiltration_depth_m, @@ -144,8 +148,8 @@ extern void cfe( } else { - fprintf(stderr,"Problem, must specify one of Schaake of Xinanjiang partitioning scheme.\n"); - fprintf(stderr,"Program terminating.\n"); + Log(FATAL, "Problem, must specify one of Schaake of Xinanjiang partitioning scheme.\n"); + Log(FATAL, "Program terminating.\n"); exit(-1); // note -1 is arbitrary #############BOMB################ NEW FLO } } @@ -166,30 +170,30 @@ extern void cfe( // LKC: should the deficit of soil be set to zero here - if the condition is reached, soil is full soil_reservoir_storage_deficit_m = 0; } - - - + + + #ifdef DEBUG /* xinanjiang_dev - printf("After Schaake function: rain:%8.5lf mm runoff:%8.5lf mm infiltration:%8.5lf mm residual:%e m\n", + Log(DEBUG, "After Schaake function: rain:%8.5lf mm runoff:%8.5lf mm infiltration:%8.5lf mm residual:%e m\n", timestep_rainfall_input_m,Schaake_output_runoff_m*1000.0,infiltration_depth_m*1000.0, timestep_rainfall_input_m-Schaake_output_runoff_m-infiltration_depth_m); */ - printf("After direct runoff function: rain:%8.5lf mm runoff:%8.5lf mm infiltration:%8.5lf mm residual:%e m\n", + Log(DEBUG, "After direct runoff function: rain:%8.5lf mm runoff:%8.5lf mm infiltration:%8.5lf mm residual:%e m\n", timestep_rainfall_input_m,infiltration_excess_m*1000.0,infiltration_depth_m*1000.0, timestep_rainfall_input_m-infiltration_excess_m-infiltration_depth_m); #endif massbal_struct->vol_runoff += infiltration_excess_m; // edit FLO massbal_struct->vol_infilt += infiltration_depth_m; // edit FLO - + // put infiltration flux into soil conceptual reservoir. If not enough room // limit amount transferred to deficit //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - + if(flux_perc_m>soil_reservoir_storage_deficit_m) { diff=flux_perc_m-soil_reservoir_storage_deficit_m; // the amount that there is not capacity ffor - infiltration_depth_m=soil_reservoir_storage_deficit_m; + infiltration_depth_m=soil_reservoir_storage_deficit_m; massbal_struct->vol_runoff+=diff; // send excess water back to GIUH runoff edit FLO massbal_struct->vol_infilt-=diff; // correct overprediction of infilt. edit FLO infiltration_excess_m+=diff; // bug found by Nels. This was missing and fixes it. @@ -198,25 +202,25 @@ extern void cfe( soil_reservoir_storage_deficit_m = 0; } - massbal_struct->vol_to_soil += infiltration_depth_m; + massbal_struct->vol_to_soil += infiltration_depth_m; soil_reservoir_struct->storage_m += infiltration_depth_m; // put the infiltrated water in the soil. - + // calculate fluxes from the soil storage into the deep groundwater (percolation) and to lateral subsurface flow //-------------------------------------------------------------------------------------------------------------- - + conceptual_reservoir_flux_calc(soil_reservoir_struct,&percolation_flux,&lateral_flux); flux_perc_m=percolation_flux; // m/h <<<<<<<<<<< flux of percolation from soil to g.w. reservoir >>>>>>>>> - + flux_lat_m=lateral_flux; // m/h <<<<<<<<<<< flux into the lateral flow Nash cascade >>>>>>>> - + // calculate flux of base flow from deep groundwater reservoir to channel //-------------------------------------------------------------------------- //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ gw_reservoir_storage_deficit_m= gw_reservoir_struct->storage_max_m-gw_reservoir_struct->storage_m; // limit amount transferred to deficit iff there is insuffienct avail. storage - if(flux_perc_m>gw_reservoir_storage_deficit_m) { + if (flux_perc_m>gw_reservoir_storage_deficit_m) { diff=flux_perc_m-gw_reservoir_storage_deficit_m; flux_perc_m=gw_reservoir_storage_deficit_m; massbal_struct->vol_runoff+=diff; // send excess water back to GIUH runoff @@ -226,7 +230,7 @@ extern void cfe( massbal_struct->vol_to_gw +=flux_perc_m; massbal_struct->vol_soil_to_gw +=flux_perc_m; - + gw_reservoir_struct->storage_m += flux_perc_m; soil_reservoir_struct->storage_m -= flux_perc_m; soil_reservoir_struct->storage_m -= flux_lat_m; @@ -235,26 +239,34 @@ extern void cfe( // get change in storage (used in the soil moisture coupler) soil_reservoir_struct->storage_change_m = soil_reservoir_struct->storage_m - storage_temp_m ; - + conceptual_reservoir_flux_calc(gw_reservoir_struct,&primary_flux,&secondary_flux); - - + + flux_from_deep_gw_to_chan_m=primary_flux; // m/h <<<<<<<<<< BASE FLOW FLUX >>>>>>>>> if(flux_from_deep_gw_to_chan_m > gw_reservoir_struct->storage_m) { flux_from_deep_gw_to_chan_m=gw_reservoir_struct->storage_m; // TODO: set a flag when flux larger than storage - printf("WARNING: Groundwater flux larger than storage \n"); + // Adding specific warning message requested by OWP (9/11/2025) + if (!groundwaterFluxIssueLogged) { + Log(WARNING, + "Groundwater flux larger than storage. " + "While this will not cause a run failure, " + "parameter combinations may not be realistic " + "and a mass balance error will occur.\n"); + groundwaterFluxIssueLogged = true; + } } - + massbal_struct->vol_from_gw+=flux_from_deep_gw_to_chan_m; - + // in the instance of calling the gw reservoir the secondary flux should be zero- verify - if(is_fabs_less_than_epsilon(secondary_flux,1.0e-09)==FALSE) printf("problem with nonzero flux point 1\n"); - + if(is_fabs_less_than_epsilon(secondary_flux,1.0e-09)==FALSE) Log(WARNING, "problem with nonzero flux point 1\n"); + // adjust state of deep groundwater conceptual nonlinear reservoir //----------------------------------------------------------------- - + gw_reservoir_struct->storage_m -= flux_from_deep_gw_to_chan_m; // Solve the convolution integral ffor this time step @@ -283,7 +295,7 @@ extern void cfe( massbal_struct->vol_out_nash += nash_lateral_runoff_m; #ifdef DEBUG - fprintf(out_debug_fptr,"%d %lf %lf %lf\n",tstep,flux_lat_m,nash_lateral_runoff_m,flux_from_deep_gw_to_chan_m); + Log(DEBUG,"%d %lf %lf %lf\n",tstep,flux_lat_m,nash_lateral_runoff_m,flux_from_deep_gw_to_chan_m); #endif Qout_m = direct_runoff_m + nash_lateral_runoff_m + flux_from_deep_gw_to_chan_m; @@ -321,7 +333,7 @@ void Schaake_partitioning_scheme(double timestep_h, double field_capacity_m, dou This subtroutine takes water_input_depth_m and partitions it into infiltration_excess_m and infiltration_depth_m using the scheme from Schaake et al. 1996. ! -------------------------------------------------------------------------------- -! ! modified by FLO April 2020 to eliminate reference to ice processes, +! ! modified by FLO April 2020 to eliminate reference to ice processes, ! ! and to de-obfuscate and use descriptive and dimensionally consistent variable names. ! -------------------------------------------------------------------------------- IMPLICIT NONE @@ -386,32 +398,32 @@ void Schaake_partitioning_scheme(double timestep_h, double field_capacity_m, dou } - // impermeable fraction due to frozen soil - double factor = 1.0; - //field_capacity_m = SMCREF in NOAH_MP - // ice_content_threshold = frzk in NOAH_MP - //double frzk = 0.15; // Ice content above which soil is impermeable - - if (ice_fraction_schaake > 1.0E-2) { - int cv_frz = 3; // should this be moved to config file as well, probably not. - double field_capacity = field_capacity_m/soil_depth; // divide by the reservior depth to get SMCREF [m/m] (unitless) - double frz_fact = smcmax/field_capacity * (0.412 / 0.468); - double frzx = ice_content_threshold * frz_fact; - - double acrt = cv_frz * frzx / ice_fraction_schaake; - double sum1 =1; - - for (int i1=1; i1 < cv_frz; i1++) { - int k = 1; - for (int i2=i1+1; i2 1.0E-2) { + int cv_frz = 3; // should this be moved to config file as well, probably not. + double field_capacity = field_capacity_m/soil_depth; // divide by the reservior depth to get SMCREF [m/m] (unitless) + double frz_fact = smcmax/field_capacity * (0.412 / 0.468); + double frzx = ice_content_threshold * frz_fact; + + double acrt = cv_frz * frzx / ice_fraction_schaake; + double sum1 =1; + + for (int i1=1; i1 < cv_frz; i1++) { + int k = 1; + for (int i2=i1+1; i2 1.0e-06) printf("Massball err. warning: %f\n", + + if(fabs(water_input_depth_m - (*infiltration_depth_m) - (*infiltration_excess_m)) > 1.0e-06) Log(DEBUG, "Massball err. warning: %f\n", water_input_depth_m - (*infiltration_depth_m) - (*infiltration_excess_m)); -#endif + return; } @@ -582,17 +595,17 @@ void et_from_rainfall(double *timestep_rainfall_input_m, struct evapotranspirati iff it is raining, take PET from rainfall first. Wet veg. is efficient evaporator. */ - if (*timestep_rainfall_input_m >0.0) { - - if (*timestep_rainfall_input_m > et_struct->potential_et_m_per_timestep) { + if (*timestep_rainfall_input_m >0.0){ + if (*timestep_rainfall_input_m > et_struct->potential_et_m_per_timestep){ + et_struct->actual_et_from_rain_m_per_timestep = et_struct->potential_et_m_per_timestep; *timestep_rainfall_input_m -= et_struct->actual_et_from_rain_m_per_timestep; } - else { + else{ // LKC: This was incorrectly set to potential instead of actual - et_struct->actual_et_from_rain_m_per_timestep = *timestep_rainfall_input_m; + et_struct->actual_et_from_rain_m_per_timestep = *timestep_rainfall_input_m; *timestep_rainfall_input_m = 0.0; } @@ -629,15 +642,15 @@ void et_from_retention_depth(struct nash_cascade_parameters *nash_surface_params void et_from_soil(struct conceptual_reservoir *soil_res, struct evapotranspiration_structure *et_struct, struct NWM_soil_parameters *soil_parms) { /* - take AET from soil moisture storage, + take AET from soil moisture storage, using Budyko type function to limit PET if wiltingactual_et_from_soil_m_per_timestep = 0; // if rootzone-based AET turned ON @@ -651,7 +664,7 @@ void et_from_soil(struct conceptual_reservoir *soil_res, struct evapotranspirati // wilting point, actual et from this timestep is 0 and no water is removed. if ( soil_res->smc_profile[soil_res->max_rootzone_layer] <= soil_parms->wltsmc ) { et_struct->actual_et_from_soil_m_per_timestep = 0; - } + } // calculate the amount of moisture removed by evapotranspiration for the bottom layer of the root zone else if (soil_res->smc_profile[soil_res->max_rootzone_layer] >= soil_res->soil_water_content_field_capacity) { et_struct->actual_et_from_soil_m_per_timestep = min(et_struct->reduced_potential_et_m_per_timestep, layer_storage_m); @@ -660,32 +673,32 @@ void et_from_soil(struct conceptual_reservoir *soil_res, struct evapotranspirati Budyko_numerator = soil_res->smc_profile[soil_res->max_rootzone_layer] - soil_parms->wltsmc; Budyko_denominator = soil_res->soil_water_content_field_capacity - soil_parms->wltsmc; Budyko = Budyko_numerator / Budyko_denominator; - + et_struct->actual_et_from_soil_m_per_timestep = min(Budyko * et_struct->reduced_potential_et_m_per_timestep, layer_storage_m); } - + // Reduce remaining PET and remove moisture from soil profile equal to the calculated AET (actual_et_from_soil_m_per_timestep - et_struct->reduced_potential_et_m_per_timestep -= et_struct->actual_et_from_soil_m_per_timestep; - soil_res->smc_profile[soil_res->max_rootzone_layer] -= (et_struct->actual_et_from_soil_m_per_timestep / + et_struct->reduced_potential_et_m_per_timestep -= et_struct->actual_et_from_soil_m_per_timestep; + soil_res->smc_profile[soil_res->max_rootzone_layer] -= (et_struct->actual_et_from_soil_m_per_timestep / soil_res->delta_soil_layer_depth_m[soil_res->max_rootzone_layer]); soil_res->storage_m -= et_struct->actual_et_from_soil_m_per_timestep; } - else if (et_struct->reduced_potential_et_m_per_timestep > 0) { - - if (soil_res->storage_m >= soil_res->storage_threshold_primary_m) { + else if (et_struct->reduced_potential_et_m_per_timestep > 0){ + + if (soil_res->storage_m >= soil_res->storage_threshold_primary_m){ et_struct->actual_et_from_soil_m_per_timestep = min(et_struct->reduced_potential_et_m_per_timestep, soil_res->storage_m); - } - else if (soil_res->storage_m > soil_parms->wilting_point_m && soil_res->storage_m < soil_res->storage_threshold_primary_m) { - + } + else if (soil_res->storage_m > soil_parms->wilting_point_m && soil_res->storage_m < soil_res->storage_threshold_primary_m){ + Budyko_numerator = soil_res->storage_m - soil_parms->wilting_point_m; Budyko_denominator = soil_res->storage_threshold_primary_m - soil_parms->wilting_point_m; Budyko = Budyko_numerator / Budyko_denominator; - // LKC: Include check to guarantee EAT is not larger than soil storage + // LKC: Include check to guarantee EAT is not larger than soil storage et_struct->actual_et_from_soil_m_per_timestep = min(Budyko * et_struct->reduced_potential_et_m_per_timestep, soil_res->storage_m); - + } - + soil_res->storage_m -= et_struct->actual_et_from_soil_m_per_timestep; et_struct->reduced_potential_et_m_per_timestep = et_struct->reduced_potential_et_m_per_timestep - et_struct->actual_et_from_soil_m_per_timestep; } @@ -695,6 +708,7 @@ void et_from_soil(struct conceptual_reservoir *soil_res, struct evapotranspirati extern int is_fabs_less_than_epsilon(double a,double epsilon) // returns true if fabs(a)coeff_primary * (exp_term - 1.0); *secondary_flux_m = 0.0; + + // cap so flux never exceeds storage, but commented out here per OWP request (9/11/2025) + //if (*primary_flux_m > da_reservoir->storage_m) { + // *primary_flux_m = da_reservoir->storage_m; + //} return; } diff --git a/src/giuh.c b/src/giuh.c index 243e8ec3..44aeb850 100644 --- a/src/giuh.c +++ b/src/giuh.c @@ -14,26 +14,26 @@ extern double giuh_convolution_integral(double runoff_m,int num_giuh_ordinates, // This function solves the convolution integral involving N // GIUH ordinates. //############################################################## - double runoff_m_current_timestep; + double runoff_m_now; int N,i; - N = num_giuh_ordinates; - runoff_queue_m_per_timestep[N] = 0.0; + N=num_giuh_ordinates; + runoff_queue_m_per_timestep[N]=0.0; for(i=0;i +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// Define Environment Variable Names +#define MODULE_NAME "CFE" +#define EV_MODULE_LOGLEVEL "CFE_LOGLEVEL" // This modules log level +#define EV_MODULE_LOGFILEPATH "CFE_LOGFILEPATH" // This modules log full log filename +#define EV_NGEN_LOGFILEPATH "NGEN_LOG_FILE_PATH" // ngen log file + +#define EV_EWTS_LOGGING "NGEN_EWTS_LOGGING" // Enable/disable of Error Warning and Trapping System + +#define DS "/" // Directory separator +#define LOG_DIR_NGENCERF "/ngencerf/data" // ngenCERF log directory +#define LOG_DIR_DEFAULT "run-logs" // Default subdir if ngen log file env not found +#define LOG_FILE_EXT "log" +#define LOG_MODULE_NAME_LEN 8 // Width of module name for log entries + +bool loggingEnabled = true; +bool openedAppendMode = true; +bool loggerInitialized = false; +FILE* logFile = NULL; +char logFilePath[1024] = ""; +LogLevel logLevel = INFO; +char moduleName[LOG_MODULE_NAME_LEN+1] = ""; + +void TrimString(const char *input, char *output, size_t outputSize) { + if (!input || !output || outputSize == 0) { + if (output && outputSize > 0) output[0] = '\0'; + return; + } + + // Skip leading whitespace + while (isspace((unsigned char)*input)) { + input++; + } + + size_t len = strlen(input); + if (len == 0) { + output[0] = '\0'; + return; + } + + // Find the end of the string and move backward past trailing whitespace + const char *end = input + len - 1; + while (end > input && isspace((unsigned char)*end)) { + end--; + } + + size_t trimmedLen = end - input + 1; + if (trimmedLen >= outputSize) { + trimmedLen = outputSize - 1; + } + + strncpy(output, input, trimmedLen); + output[trimmedLen] = '\0'; +} + +void SetLogModuleName(void) { + + // Copy MODULE_NAME to a string variable + char src[20]; + strncpy(src, MODULE_NAME, sizeof(src)); + src[sizeof(src) - 1] = '\0'; // ensure null termination + + // Convert to uppercase and copy, up to width or null terminator + size_t i = 0; + for (; i < LOG_MODULE_NAME_LEN && src[i] != '\0'; i++) { + moduleName[i] = toupper((unsigned char)src[i]); + } + + // Pad with spaces if needed + for (; i < LOG_MODULE_NAME_LEN; i++) { + moduleName[i] = ' '; + } + moduleName[LOG_MODULE_NAME_LEN] = '\0'; // null-terminate +} + +void CreateTimestamp(char *buffer, size_t size, bool appendMS, bool iso) { + struct timeval tv; + gettimeofday(&tv, NULL); + + struct tm utc_tm; + gmtime_r(&tv.tv_sec, &utc_tm); + + char base[32]; + if (iso) { + strftime(base, sizeof(base), "%Y-%m-%dT%H:%M:%S", &utc_tm); + } else { + strftime(base, sizeof(base), "%Y%m%dT%H%M%S", &utc_tm); + } + + if (appendMS) { + snprintf(buffer, size, "%s.%03ld", base, tv.tv_usec / 1000); + } else { + snprintf(buffer, size, "%s", base); + } +} + +bool DirectoryExists(const char *path) { + struct stat info; + return (stat(path, &info) == 0 && (info.st_mode & S_IFDIR)); +} + +bool CreateDirectory(const char *path) { + if (!DirectoryExists(path)) { + char cmd[512]; + snprintf(cmd, sizeof(cmd), "mkdir -p \"%s\"", path); + int status = system(cmd); + if (status == -1 || (WIFEXITED(status) && WEXITSTATUS(status) != 0)) { + fprintf(stderr, "Failed to create directory: %s\n", path); + return false; + } + } + return true; +} + +bool LogFileReady(bool appendMode) { + if (logFile && !ferror(logFile)) { + fseek(logFile, 0, SEEK_END); + return true; + } else if (strlen(logFilePath) > 0) { + logFile = fopen(logFilePath, appendMode ? "a" : "w"); + if (logFile != NULL) { + openedAppendMode = appendMode; + return true; + } + return false; + } + return false; +} + +/** + * Set the log file path name using the following pattern + * - Use the module log file if available (unset when first run by ngen), otherwise + * - Use ngen log file if available, otherwise + * - Use /ngencerf/data/run-logs//_ if available, othrewise + * - Use ~/run-logs//_ + * - Onced opened, save the full log path to the modules log environment variable so + * it is only opened once for each ngen run (vs for each catchment) + */ +void SetupLogFile(void) { + logFilePath[0] = '\0'; + bool appendEntries = true; + bool moduleLogEnvExists = false; + + const char *envVar = getenv(EV_MODULE_LOGFILEPATH); + if (envVar && envVar[0] != '\0') { + strncpy(logFilePath, envVar, sizeof(logFilePath) - 1); + moduleLogEnvExists = true; + } else { + envVar = getenv(EV_NGEN_LOGFILEPATH); + if (envVar && envVar[0] != '\0') { + strncpy(logFilePath, envVar, sizeof(logFilePath) - 1); + } else { + appendEntries = false; + char logFileDir[512]; + + if (DirectoryExists(LOG_DIR_NGENCERF)) { + snprintf(logFileDir, sizeof(logFileDir), "%s%s%s", LOG_DIR_NGENCERF, DS, LOG_DIR_DEFAULT); + } else { + const char *home = getenv("HOME"); // Get users home directly pathname + if (home) { + snprintf(logFileDir, sizeof(logFileDir), "%s%s%s", home, DS, LOG_DIR_DEFAULT); + } + else { + snprintf(logFileDir, sizeof(logFileDir), "~%s%s", DS, LOG_DIR_DEFAULT); + } + } + + if (CreateDirectory(logFileDir)) { + const char *envUsername = getenv("USER"); + if (envUsername) { + strncat(logFileDir, DS, sizeof(logFileDir) - strlen(logFileDir) - 1); + strncat(logFileDir, envUsername, sizeof(logFileDir) - strlen(logFileDir) - 1); + } else { + char dateStr[32]; + CreateTimestamp(dateStr, sizeof(dateStr), 0, 0); + strncat(logFileDir, DS, sizeof(logFileDir) - strlen(logFileDir) - 1); + strncat(logFileDir, dateStr, sizeof(logFileDir) - strlen(logFileDir) - 1); + } + + if (CreateDirectory(logFileDir)) { + char timestamp[32]; + CreateTimestamp(timestamp, sizeof(timestamp), 0, 0); + snprintf(logFilePath, sizeof(logFilePath), "%s%s%s_%s.%s", + logFileDir, DS, MODULE_NAME, timestamp, LOG_FILE_EXT); + } + } + } + } + + if (LogFileReady(appendEntries)) { + if (!moduleLogEnvExists) setenv(EV_MODULE_LOGFILEPATH, logFilePath, 1); + printf("Module %s Log File: %s\n", MODULE_NAME, logFilePath); + fflush(stdout); // Force flushing of stdout + LogLevel saveLevel = logLevel; + logLevel = INFO; // Ensure this INFO message is always logged + Log(logLevel, "Opened log file %s in %s mode\n", logFilePath, (openedAppendMode ? "Append" : "Truncate")); + logLevel = saveLevel; + } else { + printf("Unable to open log file "); + if (strlen(logFilePath) > 0) { + printf("%s (Perhaps check permissions)\n", logFilePath); + } + printf("Log entries will be written to stdout\n"); + fflush(stdout); // Force flushing of stdout + } +} + +const char* ConvertLogLevelToString(LogLevel level) { + switch (level) { + case DEBUG: return "DEBUG "; + case INFO: return "INFO "; + case WARNING: return "WARNING"; + case SEVERE: return "SEVERE "; + case FATAL: return "FATAL "; + default: return "NONE "; + } +} + +LogLevel ConvertStringToLogLevel(const char* str) { + char trimmedStr[20] = {0}; + TrimString(str, trimmedStr, sizeof(trimmedStr)); + if (strcasecmp(trimmedStr, "DEBUG") == 0) return DEBUG; + if (strcasecmp(trimmedStr, "INFO") == 0) return INFO; + if (strcasecmp(trimmedStr, "WARNING") == 0) return WARNING; + if (strcasecmp(trimmedStr, "SEVERE") == 0) return SEVERE; + if (strcasecmp(trimmedStr, "FATAL") == 0) return FATAL; + return NONE; +} + +void SetLoggingFlag(void) { + const char* ewtsStr = getenv(EV_EWTS_LOGGING); + if (ewtsStr && ewtsStr[0] != '\0') { + char trimmedStr[20] = {0}; + TrimString(ewtsStr, trimmedStr, sizeof(trimmedStr)); + bool loggingEnabled = ((strcasecmp(trimmedStr, "ENABLED") == 0))? true:false; + } + printf("%s Logging %s\n", MODULE_NAME, ((loggingEnabled)?"ENABLED":"DISABLED")); + fflush(stdout); // Force flushing of stdout +} + +void SetLogLevel(void) { + const char* envLogLevel = getenv(EV_MODULE_LOGLEVEL); + if (envLogLevel && envLogLevel[0] != '\0') { + logLevel = ConvertStringToLogLevel(envLogLevel); + } + char llStr[10]; + char msg[100]; + TrimString(ConvertLogLevelToString(logLevel), llStr, sizeof(llStr)); + snprintf(msg, sizeof(msg), "Log level set to %s\n", llStr); + printf("%s %s", MODULE_NAME, msg); + fflush(stdout); // Force flushing of stdout + LogLevel saveLevel = logLevel; + logLevel = INFO; // Ensure this INFO message is always logged + Log(logLevel, msg); + logLevel = saveLevel; +} + +void SetLogPreferences(void) { + + if (!loggerInitialized) { + loggerInitialized = true; // Only call this once + + SetLoggingFlag(); // Based on the environment variable + if (loggingEnabled) { + SetLogModuleName(); // Set the module name used in log entries + SetupLogFile(); + SetLogLevel(); + // FOR TESTING ONLY! Uncomment next line to generate test log entries + // WriteLogTestMsgs(); + } + } +} + +void TrimToOneNewline(char *str, size_t max_len) { + size_t len = strlen(str); + + // Strip trailing whitespace including \n, \r, spaces, tabs + while (len > 0 && isspace((unsigned char)str[len - 1])) { + str[--len] = '\0'; + } + + // Ensure room to add a newline + if (len + 1 < max_len) { + str[len] = '\n'; + str[len + 1] = '\0'; + } +} + +LogLevel GetLogLevel(void) { + return logLevel; +} + +bool IsLoggingEnabled(void) { + return loggingEnabled; +} + +void Log(LogLevel messageLevel, const char* message, ...) { + if (!loggerInitialized) SetLogPreferences(); // Initialize logger on first call + + if (loggingEnabled && (messageLevel >= logLevel)) { + + va_list arglist; + va_start(arglist, message); + va_list arglist_copy; + va_copy(arglist_copy, arglist); + + int length = vsnprintf(NULL, 0, message, arglist_copy); + va_end(arglist_copy); + + char *buffer = malloc(length + 1); + if (!buffer) return; // Always good to check + + vsnprintf(buffer, length + 1, message, arglist); + va_end(arglist); + + char timestamp[32]; + CreateTimestamp(timestamp, sizeof(timestamp), 1, 1); + + char logPrefix[128]; + snprintf(logPrefix, sizeof(logPrefix), "%s %s %s", timestamp, moduleName, ConvertLogLevelToString(messageLevel)); + + TrimToOneNewline(buffer, sizeof(buffer)); + char *line = strtok(buffer, "\n"); + if (LogFileReady(true)) { + while (line != NULL) { + fprintf(logFile, "%s %s\n", logPrefix, line); + line = strtok(NULL, "\n"); + } + fflush(logFile); + } else { + while (line != NULL) { + printf("%s %s\n", logPrefix, line); + line = strtok(NULL, "\n"); + } + fflush(stdout); // Force flushing of stdout + } + + free(buffer); + } +} + +void WriteLogTestMsgs() { + + LogLevel savedLogLevel = logLevel; + + setenv(EV_MODULE_LOGLEVEL, "SEVERE", 1); Log(SEVERE, "Sample Log for LogLevel::SEVERE"); + setenv(EV_MODULE_LOGLEVEL, "FATAL", 1); Log(FATAL, "Sample Log for LogLevel::FATAL"); + setenv(EV_MODULE_LOGLEVEL, "WARNING", 1); Log(WARNING, "Sample Log for LogLevel::WARNING"); + setenv(EV_MODULE_LOGLEVEL, "INFO", 1); Log(INFO, "Sample Log for LogLevel::INFO"); + setenv(EV_MODULE_LOGLEVEL, "DEBUG", 1); Log(DEBUG, "Sample Log for LogLevel::DEBUG"); + + const char* multiline_log = + "First line of multiline log:\n" + "Second line of multiline log\n" + "Third line of multiline log\n" + "Fourth line of multiline log"; + Log(DEBUG, multiline_log); + + setenv(EV_MODULE_LOGLEVEL, ConvertLogLevelToString(savedLogLevel), 1); +} diff --git a/src/main_cfe_aorc_pet_rz_aet.cxx b/src/main_cfe_aorc_pet_rz_aet.cxx index 54f5f838..4cf76b9b 100644 --- a/src/main_cfe_aorc_pet_rz_aet.cxx +++ b/src/main_cfe_aorc_pet_rz_aet.cxx @@ -12,6 +12,7 @@ #include "../extern/evapotranspiration/include/pet.h" #include "../extern/evapotranspiration/include/bmi_pet.h" +#include "../bmi/bmi.hxx" #include "../extern/SoilMoistureProfiles/include/bmi_soil_moisture_profile.hxx" #include "../extern/SoilMoistureProfiles/include/soil_moisture_profile.hxx" diff --git a/src/main_pass_forcings.c b/src/main_pass_forcings.c index 80c01031..ce774bf5 100644 --- a/src/main_pass_forcings.c +++ b/src/main_pass_forcings.c @@ -128,7 +128,7 @@ int /************************************************************************ Finalize both the CFE and AORC bmi models ************************************************************************/ - printf("Finalizing CFE and AORC models\n"); + printf("Finalizing BFE and AORC models\n"); cfe_bmi_model->finalize(cfe_bmi_model); aorc_bmi_model->finalize(aorc_bmi_model);