diff --git a/CMakeLists.txt b/CMakeLists.txt index 6ceb4a3b..38840828 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -47,15 +47,15 @@ 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) @@ -63,6 +63,7 @@ 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) diff --git a/README.md b/README.md index 491fa19c..d89ca0d0 100644 --- a/README.md +++ b/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,7 +79,7 @@ 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` @@ -88,13 +87,13 @@ To compile and run CFE with locally read forcing data, run the following from th 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 diff --git a/include/cfe.h b/include/cfe.h index a71ca898..03f8809e 100644 --- a/include/cfe.h +++ b/include/cfe.h @@ -7,7 +7,7 @@ #include #include #include - +#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. @@ -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); diff --git a/include/conceptual_reservoir.h b/include/conceptual_reservoir.h new file mode 100644 index 00000000..88de7caa --- /dev/null +++ b/include/conceptual_reservoir.h @@ -0,0 +1,47 @@ +#ifndef _CONCEPTUAL_RESERVOIR_H +#define _CONCEPTUAL_RESERVOIR_H + +#include +#include +#include + +#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 diff --git a/include/giuh.h b/include/giuh.h index 99fb8726..82ce2f34 100644 --- a/include/giuh.h +++ b/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); diff --git a/src/bmi_cfe.c b/src/bmi_cfe.c index b85421ef..c000632d 100644 --- a/src/bmi_cfe.c +++ b/src/bmi_cfe.c @@ -1,8 +1,8 @@ #include #include #include -#include "../include/bmi.h" -#include "../include/bmi_cfe.h" +#include "bmi.h" +#include "bmi_cfe.h" #include #include #ifndef WATER_SPECIFIC_WEIGHT diff --git a/src/cfe.c b/src/cfe.c index 6079a19d..d1cee9dd 100644 --- a/src/cfe.c +++ b/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 <<<>>> -//{ -// 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 ############# diff --git a/src/conceptual_reservoir.c b/src/conceptual_reservoir.c new file mode 100644 index 00000000..496429ce --- /dev/null +++ b/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 <<<>>> + //{ + // 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 diff --git a/src/giuh.c b/src/giuh.c index f9cc0121..44aeb850 100644 --- a/src/giuh.c +++ b/src/giuh.c @@ -1,8 +1,7 @@ #ifndef _GIUH_C - #define _GIUH_C -#include "../include/giuh.h" +#include "giuh.h" //############################################################## diff --git a/src/main.c b/src/main.c index 6cfb70dc..105bd340 100644 --- a/src/main.c +++ b/src/main.c @@ -1,9 +1,9 @@ #include #include #include -#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. diff --git a/src/main_cfe_aorc_pet.c b/src/main_cfe_aorc_pet.c index daee76cc..49b02f03 100644 --- a/src/main_cfe_aorc_pet.c +++ b/src/main_cfe_aorc_pet.c @@ -1,9 +1,9 @@ #include #include #include -#include "../include/bmi.h" -#include "../include/cfe.h" -#include "../include/bmi_cfe.h" +#include "bmi.h" +#include "cfe.h" +#include "bmi_cfe.h" #include "../forcing_code/include/pet.h" #include "../forcing_code/include/bmi_pet.h" #include "../forcing_code/include/aorc.h" diff --git a/src/main_cfe_aorc_pet_rz_aet.cxx b/src/main_cfe_aorc_pet_rz_aet.cxx index 308e668b..2bd27bbb 100644 --- a/src/main_cfe_aorc_pet_rz_aet.cxx +++ b/src/main_cfe_aorc_pet_rz_aet.cxx @@ -3,9 +3,9 @@ #include -#include "../include/cfe.h" -#include "../include/bmi.h" -#include "../include/bmi_cfe.h" +#include "cfe.h" +#include "bmi.h" +#include "bmi_cfe.h" #include "../forcing_code/include/pet.h" #include "../forcing_code/include/bmi_pet.h" diff --git a/src/main_pass_forcings.c b/src/main_pass_forcings.c index 31bf7a3c..3187f795 100644 --- a/src/main_pass_forcings.c +++ b/src/main_pass_forcings.c @@ -1,9 +1,9 @@ #include #include #include -#include "../include/cfe.h" -#include "../include/bmi.h" -#include "../include/bmi_cfe.h" +#include "cfe.h" +#include "bmi.h" +#include "bmi_cfe.h" #include "../forcing_code/include/aorc.h" #include "../forcing_code/include/bmi_aorc.h"