refactored cfe source file

CMakeLists.txt

@@ -47,22 +47,23 @@ set(CMAKE_C_FLAGS_DEBUG "${CMAKE_C_FLAGS_DEBUG} -O0")
 
 ## 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 ./forcing_code/src/aorc.c ./forcing_code/src/bmi_aorc.c ./forcing_code/src/pet.c ./forcing_code/src/bmi_pet.c)
+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 ./forcing_code/src/aorc.c ./forcing_code/src/bmi_aorc.c ./forcing_code/src/pet.c ./forcing_code/src/bmi_pet.c)
 
 add_library(cfelib ./smc_coupler/src/bmi_soil_moisture_profile.cxx ./smc_coupler/src/soil_moisture_profile.cxx ./smc_coupler/include/bmi_soil_moisture_profile.hxx ./smc_coupler/include/soil_moisture_profile.hxx)
 elseif(FORCING)
-add_executable(${exe_name} ./src/main_pass_forcings.c ./src/cfe.c ./src/bmi_cfe.c ./src/giuh.c ./forcing_code/src/aorc.c ./forcing_code/src/bmi_aorc.c)
+add_executable(${exe_name} ./src/main_pass_forcings.c ./src/cfe.c ./src/bmi_cfe.c ./src/giuh.c ./src/conceptual_reservoir.c ./forcing_code/src/aorc.c ./forcing_code/src/bmi_aorc.c)
 elseif(FORCINGPET)
-add_executable(${exe_name} ./src/main_cfe_aorc_pet.c ./src/cfe.c ./src/bmi_cfe.c ./src/giuh.c ./forcing_code/src/aorc.c ./forcing_code/src/bmi_aorc.c ./forcing_code/src/pet.c ./forcing_code/src/bmi_pet.c)
+add_executable(${exe_name} ./src/main_cfe_aorc_pet.c ./src/cfe.c ./src/bmi_cfe.c ./src/giuh.c ./src/conceptual_reservoir.c ./forcing_code/src/aorc.c ./forcing_code/src/bmi_aorc.c ./forcing_code/src/pet.c ./forcing_code/src/bmi_pet.c)
 else()
-add_executable(${exe_name} ./src/main.c ./src/cfe.c ./src/bmi_cfe.c ./src/giuh.c)
+add_executable(${exe_name} ./src/main.c ./src/cfe.c ./src/bmi_cfe.c ./src/giuh.c ./src/conceptual_reservoir.c)
 
 endif()
 if(AETROOTZONE)
 target_link_libraries(${exe_name} LINK_PUBLIC cfelib)
 endif()
 target_link_libraries(${exe_name} PRIVATE m)
 
+target_include_directories(${exe_name} PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/include)
 
 unset(BASE CACHE)
 unset(FORCING CACHE)

README.md

@@ -19,7 +19,6 @@ There are four examples for running CFE as described below. They assume you have
 ````
 git clone https://github.com/NOAA-OWP/cfe.git
 cd cfe
-git checkout ajk/sft_aet_giuh_merge (update after PR merge)
 git clone https://github.com/NOAA-OWP/SoilMoistureProfiles.git smc_coupler (needed if AETROOTZONE=ON)
 mkdir build && cd build
 cmake ../ [-DBASE=ON,-DFORCING=ON,-DFORCINGPET=ON,-DAETROOTZONE=ON] (pick one option, e.g. `cmake ../ -DFORCING=ON`)
@@ -80,21 +79,21 @@ A [configs/](./configs/) directory contains primiary configuration text files fo
 ### 1. Read local forcing file
 To compile and run CFE with locally read forcing data, run the following from the command line:
 
-1. `gcc -lm ./src/main.c ./src/cfe.c ./src/bmi_cfe.c ./src/giuh.c -o cfe_base`. This will generate an executable called `cfe_base`.
+1. `gcc -lm -Iinclude ./src/main.c ./src/cfe.c ./src/bmi_cfe.c ./src/giuh.c ./src/conceptual_reservoir.c -o cfe_base`. This will generate an executable called `cfe_base`.
 2.  Then run the model with example forcing data: `./cfe_base ./configs/cat_58_bmi_config_cfe.txt`  
 
 
 ### 2. CFE Model gets forcings passed from BMI
 
 CFE was designed to read its own forcing file, but we have added an option to get forcings passed in through BMI using its `set_value` functionality. To demonstrate this functionality we have included the BMI-enabled AORC forcing read module. Follow the steps below:  
 
-1. `gcc -lm ./src/main_pass_forcings.c ./src/cfe.c ./src/bmi_cfe.c ./src/giuh.c ./forcing_code/src/aorc.c ./forcing_code/src/bmi_aorc.c  -o cfe_forcing`. This generates an executable called `cfe_forcing`. 
+1. `gcc -lm -Iinclude ./src/main_pass_forcings.c ./src/cfe.c ./src/bmi_cfe.c ./src/giuh.c ./src/conceptual_reservoir.c ./forcing_code/src/aorc.c ./forcing_code/src/bmi_aorc.c  -o cfe_forcing`. This generates an executable called `cfe_forcing`. 
 2. To run this executable you must pass the path to the corresponding configuration files for **BOTH** CFE and AORC (in that order): `./cfe_forcing ./configs/cat_89_bmi_config_cfe_pass.txt ./configs/cat_89_bmi_config_aorc.txt`
 
 ### 3. CFE Model gets forcings AND potential evapotranspiration passed from BMI
 CFE can remove mass from the modeled system through evapotranspiration (directly from precipitation and from the soil using the Budyko function). Follow the steps below:  
 
-1. `gcc -lm ./src/main_cfe_aorc_pet.c ./forcing_code/src/pet.c ./forcing_code/src/bmi_pet.c ./src/cfe.c ./src/bmi_cfe.c ./src/giuh.c ./forcing_code/src/aorc.c ./forcing_code/src/bmi_aorc.c -o cfe_forcingpet`. This generates an executable called `cfe_forcingpet`.
+1. `gcc -lm -Iinclude ./src/main_cfe_aorc_pet.c ./forcing_code/src/pet.c ./forcing_code/src/bmi_pet.c ./src/cfe.c ./src/bmi_cfe.c ./src/giuh.c ./src/conceptual_reservoir.c ./forcing_code/src/aorc.c ./forcing_code/src/bmi_aorc.c -o cfe_forcingpet`. This generates an executable called `cfe_forcingpet`.
 2. To run this executable you must pass the path to the corresponding configuration files for CFE, PET and AORC (in that order):  `./cfe_forcingpet ./configs/cat_89_bmi_config_cfe_pass.txt ./configs/cat_89_bmi_config_aorc.txt ./configs/cat_89_bmi_config_pet_pass.txt`
 
 ### 4. CFE rootzone-based example couples C and C++ modules and can't be built without cmake

include/cfe.h

@@ -7,7 +7,7 @@
 #include <float.h>
 #include <string.h>
 #include <assert.h>
-
+#include "conceptual_reservoir.h"
 #include "giuh.h"
 
 #define TRUE 1
@@ -69,40 +69,6 @@
 // define data structures
 //--------------------------
 
-struct conceptual_reservoir
-{
-// this data structure describes a nonlinear reservoir having two outlets, one primary with an activation
-// threshold that may be zero, and a secondary outlet with a threshold that may be zero
-// this will also simulate a linear reservoir by setting the exponent parameter to 1.0 iff is_exponential==FALSE
-// iff is_exponential==TRUE, then it uses the exponential discharge function from the NWM V2.0 forumulation
-// as the primary discharge with a zero threshold, and does not calculate a secondary discharge.
-//--------------------------------------------------------------------------------------------------
-int    is_exponential;  // set this true TRUE to use the exponential form of the discharge equation
-double gw_storage;   // Initial Storage - LKC: added since I need to keep track of it when changing parameters
-double storage_max_m;   // maximum storage in this reservoir
-double storage_m;       // state variable.
-double storage_change_m; // storage change in the current step
-double coeff_primary;    // the primary outlet
-double exponent_primary;
-double storage_threshold_primary_m;
-double storage_threshold_secondary_m;
-double coeff_secondary;
-double exponent_secondary;
-double ice_fraction_schaake, ice_fraction_xinan;
-int   is_sft_coupled; // boolean - true if SFT is ON otherwise OFF (default is OFF)
-
-//---Root zone adjusted AET development -rlm -ahmad -------------
-double *smc_profile; //soil moisture content profile
-int n_soil_layers; // number of soil layers
-double *soil_layer_depths_m; // soil layer depths defined in the config file in units of [m]
-int aet_root_zone; // boolean - true if aet_root_zone is ON otherwise OFF (default is OFF)
-int max_root_zone_layer;  // maximum root zone layer is used to identify the maximum layer to remove water from for AET
-double *delta_soil_layer_depth_m; // used to calculate the total soil moisture in each layer
-double soil_water_content_field_capacity;  // water content [m/m] at field capacity.  Used in AET routine 
-
-//---------------------------------------------------------------
-};
-
 struct NWM_soil_parameters
 {
 // using same variable names as used in NWM.  <sorry>
@@ -187,9 +153,6 @@ extern void Xinanjiang_partitioning_scheme(double water_input_depth_m, double fi
 					   struct direct_runoff_parameters_structure *parms, double *surface_runoff_depth_m,
 					   double *infiltration_depth_m, double ice_fraction_xinan);
 
-extern void conceptual_reservoir_flux_calc(struct conceptual_reservoir *da_reservoir,
-                                           double *primary_flux_m, double *secondary_flux_m);
-
 extern double convolution_integral(double runoff_m, int num_giuh_ordinates, 
                                    double *giuh_ordinates, double *runoff_queue_m_per_timestep);
 

include/conceptual_reservoir.h

@@ -0,0 +1,47 @@
+#ifndef _CONCEPTUAL_RESERVOIR_H
+#define _CONCEPTUAL_RESERVOIR_H
+
+#include <stdio.h>
+#include <math.h>
+#include <stdlib.h>
+
+#define TRUE 1
+#define FALSE 0
+
+struct conceptual_reservoir {
+  // this data structure describes a nonlinear reservoir having two outlets, one primary with an activation
+  // threshold that may be zero, and a secondary outlet with a threshold that may be zero
+  // this will also simulate a linear reservoir by setting the exponent parameter to 1.0 iff is_exponential==FALSE
+  // iff is_exponential==TRUE, then it uses the exponential discharge function from the NWM V2.0 forumulation
+  // as the primary discharge with a zero threshold, and does not calculate a secondary discharge.
+  //--------------------------------------------------------------------------------------------------
+  int    is_exponential;  // set this true TRUE to use the exponential form of the discharge equation
+  double gw_storage;   // Initial Storage - LKC: added since I need to keep track of it when changing parameters
+  double storage_max_m;   // maximum storage in this reservoir
+  double storage_m;       // state variable.
+  double storage_change_m; // storage change in the current step
+  double coeff_primary;    // the primary outlet
+  double exponent_primary;
+  double storage_threshold_primary_m;
+  double storage_threshold_secondary_m;
+  double coeff_secondary;
+  double exponent_secondary;
+  double ice_fraction_schaake, ice_fraction_xinan;
+  int   is_sft_coupled; // boolean - true if SFT is ON otherwise OFF (default is OFF)
+
+  //---Root zone adjusted AET development -rlm -ahmad -------------
+  double *smc_profile; //soil moisture content profile
+  int n_soil_layers; // number of soil layers
+  double *soil_layer_depths_m; // soil layer depths defined in the config file in units of [m]
+  int aet_root_zone; // boolean - true if aet_root_zone is ON otherwise OFF (default is OFF)
+  int max_root_zone_layer;  // maximum root zone layer is used to identify the maximum layer to remove water from for AET
+  double *delta_soil_layer_depth_m; // used to calculate the total soil moisture in each layer
+  double soil_water_content_field_capacity;  // water content [m/m] at field capacity.  Used in AET routine 
+
+  //---------------------------------------------------------------
+};
+
+extern void conceptual_reservoir_flux_calc(struct conceptual_reservoir *da_reservoir,
+                                           double *primary_flux_m, double *secondary_flux_m);
+
+#endif

include/giuh.h

@@ -1,8 +1,6 @@
 #ifndef _GIUH_H
-
 #define _GIUH_H
 
-
 extern double giuh_convolution_integral(double runoff_m, int num_giuh_ordinates, 
                                    double *giuh_ordinates, double *runoff_queue_m_per_timestep);
 

src/bmi_cfe.c

@@ -1,8 +1,8 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <stdbool.h>
-#include "../include/bmi.h"
-#include "../include/bmi_cfe.h"
+#include "bmi.h"
+#include "bmi_cfe.h"
 #include <time.h>
 #include <float.h>
 #ifndef WATER_SPECIFIC_WEIGHT

src/cfe.c

@@ -1,4 +1,5 @@
-#include "../include/cfe.h"
+#include "cfe.h"
+
 #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; })
 
@@ -325,81 +326,6 @@ return (outflow_m);
 }
 
 
-//##############################################################
-//########## SINGLE OUTLET EXPONENTIAL RESERVOIR ###############
-//##########                -or-                 ###############
-//##########    TWO OUTLET NONLINEAR RESERVOIR   ###############
-//################################################################
-// This function calculates the flux from a linear, or nonlinear 
-// conceptual reservoir with one or two outlets, or from an
-// exponential nonlinear conceptual reservoir with only one outlet.
-// In the non-exponential instance, each outlet can have its own
-// activation storage threshold.  Flow from the second outlet is 
-// turned off by setting the discharge coeff. to 0.0.
-//################################################################
-extern void conceptual_reservoir_flux_calc(struct conceptual_reservoir *da_reservoir,
-                                           double *primary_flux_m,double *secondary_flux_m)
-{
-//struct conceptual_reservoir  <<<<INCLUDED HERE FOR REFERENCE.>>>>
-//{
-// int    is_exponential;  // set this true TRUE to use the exponential form of the discharge equation
-// double storage_max_m;
-// double storage_m;
-// double coeff_primary;
-// double exponent_secondary;
-// double storage_threshold_primary_m;
-// double storage_threshold_secondary_m;
-// double coeff_secondary;
-// double exponent_secondary;
-// };
-// THIS FUNCTION CALCULATES THE FLUXES FROM A CONCEPTUAL NON-LINEAR (OR LINEAR) RESERVOIR WITH TWO OUTLETS
-// all fluxes calculated by this routine are instantaneous with units of the coefficient.
-
-//local variables
-double storage_above_threshold_m;
-
-if(da_reservoir->is_exponential==TRUE)  // single outlet reservoir like the NWM V1.2 exponential conceptual gw reservoir
-  {
-  // calculate the one flux and return.
-  *primary_flux_m=da_reservoir->coeff_primary*
-                    (exp(da_reservoir->exponent_primary*da_reservoir->storage_m/da_reservoir->storage_max_m)-1.0);
-
-  *secondary_flux_m=0.0;
-  /*printf("primary_flux_m %lf  \n",*primary_flux_m);
-  printf("da_reservoir->coeff_primary %lf \n",da_reservoir->coeff_primary);
-  printf("da_reservoir->exponent_primary %lf \n",da_reservoir->exponent_primary);
-  printf("da_reservoir->storage_m %lf \n",da_reservoir->storage_m);
-  printf("da_reservoir->storage_max_m %lf \n",da_reservoir->storage_max_m);*/
-  return;
-  }
-// code goes past here iff it is not a single outlet exponential deep groundwater reservoir of the NWM variety
-// The vertical outlet is assumed to be primary and satisfied first.
-
-*primary_flux_m=0.0;
-storage_above_threshold_m=da_reservoir->storage_m-da_reservoir->storage_threshold_primary_m;
-if(storage_above_threshold_m>0.0)
-  {
-  // flow is possible from the primary outlet
-  *primary_flux_m=da_reservoir->coeff_primary*
-                pow(storage_above_threshold_m/(da_reservoir->storage_max_m-da_reservoir->storage_threshold_primary_m),
-                    da_reservoir->exponent_primary);
-  if(*primary_flux_m > storage_above_threshold_m) 
-                    *primary_flux_m=storage_above_threshold_m;  // limit to max. available
-  }
-*secondary_flux_m=0.0;
-storage_above_threshold_m=da_reservoir->storage_m-da_reservoir->storage_threshold_secondary_m;
-if(storage_above_threshold_m>0.0)
-  {
-  // flow is possible from the secondary outlet
-  *secondary_flux_m=da_reservoir->coeff_secondary*
-                  pow(storage_above_threshold_m/(da_reservoir->storage_max_m-da_reservoir->storage_threshold_secondary_m),
-                      da_reservoir->exponent_secondary);
-  if(*secondary_flux_m > (storage_above_threshold_m-(*primary_flux_m))) 
-                    *secondary_flux_m=storage_above_threshold_m-(*primary_flux_m);  // limit to max. available
-  }
-return;
-}
-
 
 //##############################################################
 //#########   SCHAAKE RUNOFF PARTITIONING SCHEME   #############

src/conceptual_reservoir.c

@@ -0,0 +1,86 @@
+#ifndef _CONCEPTUAL_RESERVOIR_C
+#define _CONCEPTUAL_RESERVOIR_C
+
+#include "conceptual_reservoir.h"
+
+
+//##############################################################
+//########## SINGLE OUTLET EXPONENTIAL RESERVOIR ###############
+//##########                -or-                 ###############
+//##########    TWO OUTLET NONLINEAR RESERVOIR   ###############
+//################################################################
+// This function calculates the flux from a linear, or nonlinear 
+// conceptual reservoir with one or two outlets, or from an
+// exponential nonlinear conceptual reservoir with only one outlet.
+// In the non-exponential instance, each outlet can have its own
+// activation storage threshold.  Flow from the second outlet is 
+// turned off by setting the discharge coeff. to 0.0.
+//################################################################
+extern void conceptual_reservoir_flux_calc(struct conceptual_reservoir *da_reservoir,
+                                           double *primary_flux_m,double *secondary_flux_m)
+{
+  //struct conceptual_reservoir  <<<<INCLUDED HERE FOR REFERENCE.>>>>
+  //{
+  // int    is_exponential;  // set this true TRUE to use the exponential form of the discharge equation
+  // double storage_max_m;
+  // double storage_m;
+  // double coeff_primary;
+  // double exponent_secondary;
+  // double storage_threshold_primary_m;
+  // double storage_threshold_secondary_m;
+  // double coeff_secondary;
+  // double exponent_secondary;
+  // };
+  // THIS FUNCTION CALCULATES THE FLUXES FROM A CONCEPTUAL NON-LINEAR (OR LINEAR) RESERVOIR WITH TWO OUTLETS
+  // all fluxes calculated by this routine are instantaneous with units of the coefficient.
+
+  //local variables
+  double storage_above_threshold_m;
+
+  // *****************************************************************************
+  // ------------------ Conceptual Ground Water Reservoir -------------------------
+  // single outlet reservoir like the NWM V1.2 exponential conceptual gw reservoir
+  if (da_reservoir->is_exponential == TRUE) { 
+    // calculate the one flux and return.
+    double exp_term = exp(da_reservoir->exponent_primary * da_reservoir->storage_m / da_reservoir->storage_max_m);
+
+    *primary_flux_m = da_reservoir->coeff_primary * (exp_term - 1.0);
+    *secondary_flux_m = 0.0;
+
+    return;
+  }
+  // *****************************************************************************
+
+  // code goes past here iff it is not a single outlet exponential deep groundwater reservoir of the NWM variety
+  // The vertical outlet is assumed to be primary and satisfied first.
+
+  *primary_flux_m=0.0;
+  storage_above_threshold_m = da_reservoir->storage_m - da_reservoir->storage_threshold_primary_m;
+
+  if (storage_above_threshold_m > 0.0) {
+    // flow is possible from the primary outlet
+    *primary_flux_m = da_reservoir->coeff_primary *
+      pow(storage_above_threshold_m / (da_reservoir->storage_max_m - da_reservoir->storage_threshold_primary_m),
+	  da_reservoir->exponent_primary);
+
+    if(*primary_flux_m > storage_above_threshold_m) 
+      *primary_flux_m = storage_above_threshold_m;  // limit to max. available
+  }
+
+  *secondary_flux_m=0.0;
+  storage_above_threshold_m = da_reservoir->storage_m - da_reservoir->storage_threshold_secondary_m;
+
+  if (storage_above_threshold_m > 0.0) {
+    // flow is possible from the secondary outlet
+    *secondary_flux_m = da_reservoir->coeff_secondary *
+      pow(storage_above_threshold_m / (da_reservoir->storage_max_m - da_reservoir->storage_threshold_secondary_m),
+	  da_reservoir->exponent_secondary);
+
+    if (*secondary_flux_m > (storage_above_threshold_m-(*primary_flux_m))) 
+      *secondary_flux_m = storage_above_threshold_m-(*primary_flux_m);  // limit to max. available
+  }
+  return;
+}
+
+
+#endif

src/giuh.c

@@ -1,8 +1,7 @@
 #ifndef _GIUH_C
-
 #define _GIUH_C
 
-#include "../include/giuh.h"
+#include "giuh.h"
 
 
 //##############################################################

src/main.c

@@ -1,9 +1,9 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <stdbool.h>
-#include "../include/bmi.h"
-#include "../include/bmi_cfe.h"
-#include "../include/cfe.h"
+#include "bmi.h"
+#include "bmi_cfe.h"
+#include "cfe.h"
 
 /*
 This main program is a mock framwork.