diff --git a/README.md b/README.md index ac6426a3..de2700cb 100644 --- a/README.md +++ b/README.md @@ -6,6 +6,15 @@ There are actually multiple ways to run the CFE. 1. As written by the original author. This includes a full program to read and process atmospheric forcing data, print the model output and check for mass balance closure. This code can be run from the `original_author_code` directory. 2. Through the BMI commands. This method will be used by the Next Generation U.S. National Water Model (Nextgen NWM). The `make_and_run_bmi.sh` script in the main directory does that for you through a main run file in the `./src` directory. +## Direct runoff +The user has the option to pick a particular direct runoff (aka surface partitioning) method: +1. Schaake function (configuration: `surface_partitioning_scheme=Schaake`) +2. XinanJiang function (configuration: `surface_partitioning_scheme=Xinanjiang`). When using this runoff method the user must also include three parameters. +### If XinanJiang is choosen these parameters need to be included in the configuration file: +1. a_Xinanjiang_inflection_point_parameter +2. b_Xinanjiang_shape_parameter +3. x_Xinanjiang_shape_parameter + # Compiling this code to run examples with a "pseudo" or "mini" framework. ## Read local forcing file The BMI functionality was developed as a standalone module in C. To compile this code the developer used these steps: diff --git a/configs/cat_58_bmi_config_cfe.txt b/configs/cat_58_bmi_config_cfe.txt index 59ac8a40..3468803a 100644 --- a/configs/cat_58_bmi_config_cfe.txt +++ b/configs/cat_58_bmi_config_cfe.txt @@ -19,3 +19,4 @@ nash_storage=0.0,0.0 giuh_ordinates=0.06,0.51,0.28,0.12,0.03 num_timesteps=1 verbosity=1 +surface_partitioning_scheme=Schaake \ No newline at end of file diff --git a/configs/cat_58_bmi_config_cfe_pass.txt b/configs/cat_58_bmi_config_cfe_pass.txt index b515018e..6531feb5 100644 --- a/configs/cat_58_bmi_config_cfe_pass.txt +++ b/configs/cat_58_bmi_config_cfe_pass.txt @@ -18,4 +18,5 @@ K_lf=0.01 nash_storage=0.0,0.0 giuh_ordinates=0.06,0.51,0.28,0.12,0.03 num_timesteps=1 -verbosity=1 \ No newline at end of file +verbosity=1 +surface_partitioning_scheme=Schaake \ No newline at end of file diff --git a/configs/cat_87_bmi_config_aorc.txt b/configs/cat_87_bmi_config_aorc.txt index 265c3fc0..3a4ca72e 100644 --- a/configs/cat_87_bmi_config_aorc.txt +++ b/configs/cat_87_bmi_config_aorc.txt @@ -1,5 +1,5 @@ verbose=0 -forcing_file=./forcings/cat87_01Dec2015-.csv +forcing_file=./forcings/cat87_01Dec2015.csv run_unit_tests=0 wind_speed_measurement_height_m=10.0 humidity_measurement_height_m=2.0 diff --git a/configs/cat_87_bmi_config_cfe.txt b/configs/cat_87_bmi_config_cfe.txt index 1b5d1f0b..e2cb9973 100644 --- a/configs/cat_87_bmi_config_cfe.txt +++ b/configs/cat_87_bmi_config_cfe.txt @@ -18,4 +18,8 @@ K_lf=0.01 nash_storage=0.0,0.0 giuh_ordinates=0.06,0.51,0.28,0.12,0.03 num_timesteps=1 -verbosity=1 \ No newline at end of file +verbosity=2 +surface_partitioning_scheme=Xinanjiang +a_Xinanjiang_inflection_point_parameter=1 +b_Xinanjiang_shape_parameter=1 +x_Xinanjiang_shape_parameter=1 \ No newline at end of file diff --git a/configs/cat_87_bmi_config_cfe_pass.txt b/configs/cat_87_bmi_config_cfe_pass.txt index 1b5d1f0b..5a04f299 100644 --- a/configs/cat_87_bmi_config_cfe_pass.txt +++ b/configs/cat_87_bmi_config_cfe_pass.txt @@ -1,4 +1,4 @@ -forcing_file=./forcings/cat87_01Dec2015.csv +forcing_file=BMI soil_params.depth=2.0 soil_params.b=4.05 soil_params.mult=1000.0 @@ -18,4 +18,8 @@ K_lf=0.01 nash_storage=0.0,0.0 giuh_ordinates=0.06,0.51,0.28,0.12,0.03 num_timesteps=1 -verbosity=1 \ No newline at end of file +verbosity=1 +surface_partitioning_scheme=Xinanjiang +a_Xinanjiang_inflection_point_parameter=1 +b_Xinanjiang_shape_parameter=1 +x_Xinanjiang_shape_parameter=1 \ No newline at end of file diff --git a/configs/cat_89_bmi_config_cfe.txt b/configs/cat_89_bmi_config_cfe.txt index 970ec950..b3ef9cd3 100644 --- a/configs/cat_89_bmi_config_cfe.txt +++ b/configs/cat_89_bmi_config_cfe.txt @@ -19,3 +19,4 @@ nash_storage=0.0,0.0 giuh_ordinates=0.06,0.51,0.28,0.12,0.03 num_timesteps=1 verbosity=1 +surface_partitioning_scheme=Schaake \ No newline at end of file diff --git a/configs/cat_89_bmi_config_cfe_pass.txt b/configs/cat_89_bmi_config_cfe_pass.txt index 764f04ca..c13d8e3e 100644 --- a/configs/cat_89_bmi_config_cfe_pass.txt +++ b/configs/cat_89_bmi_config_cfe_pass.txt @@ -19,3 +19,4 @@ nash_storage=0.0,0.0 giuh_ordinates=0.06,0.51,0.28,0.12,0.03 num_timesteps=1 verbosity=1 +surface_partitioning_scheme=Schaake diff --git a/configs/cat_89_bmi_config_cfe_unit_test.txt b/configs/cat_89_bmi_config_cfe_unit_test.txt index 344646e6..cced8b49 100644 --- a/configs/cat_89_bmi_config_cfe_unit_test.txt +++ b/configs/cat_89_bmi_config_cfe_unit_test.txt @@ -19,3 +19,4 @@ nash_storage=0.0,0.0 giuh_ordinates=0.06,0.51,0.28,0.12,0.03 num_timesteps=1 verbosity=0 +surface_partitioning_scheme=Schaake diff --git a/forcing_code/src/aorc.c b/forcing_code/src/aorc.c index 8a78aed0..1713af81 100644 --- a/forcing_code/src/aorc.c +++ b/forcing_code/src/aorc.c @@ -71,7 +71,7 @@ extern int run_aorc(aorc_model* model) aorc_calculate_solar_radiation(model); if (model->bmi.verbose > 1) - printf("calculate the net radiation"); + printf("calculate the net radiation\n"); // NOTE don't call this function use_aerodynamic_method option is TRUE model->other_forcing.net_radiation_W_per_sq_m=aorc_calculate_net_radiation_W_per_sq_m(model); @@ -148,79 +148,67 @@ void aorc_unit_tests(aorc_model* model) /*####################################################################*/ /*########################### PARSE LINE #############################*/ /*####################################################################*/ -void parse_aorc_line_aorc(char *theString, long *year, long *month, long *day, long *hour, long *minute, double *second, - struct aorc_forcing_data *aorc) { - char str[20]; - long yr, mo, da, hr, mi; - double mm, julian, se; - double val; - int i, start, end, len; - int yes_pm, wordlen; +void parse_aorc_line_aorc(char *theString,long *year,long *month, long *day,long *hour,long *minute, double *second, + struct aorc_forcing_data *aorc) +{ + int start,end; + int wordlen; char theWord[150]; - - len = strlen(theString); - - char *copy, *copy_to_free, *value; - copy_to_free = copy = strdup(theString); - - // time - value = strsep(©, ","); - // TODO: handle this - // struct tm{ - // int tm_year; - // int tm_mon; - // int tm_mday; - // int tm_hour; - // int tm_min; - // int tm_sec; - // int tm_isdst; - // } t; - struct tm t; - time_t t_of_day; - - t.tm_year = (int)strtol(strsep(&value, "-"), NULL, 10) - 1900; - t.tm_mon = (int)strtol(strsep(&value, "-"), NULL, 10); - t.tm_mday = (int)strtol(strsep(&value, " "), NULL, 10); - t.tm_hour = (int)strtol(strsep(&value, ":"), NULL, 10); - t.tm_min = (int)strtol(strsep(&value, ":"), NULL, 10); - t.tm_sec = (int)strtol(value, NULL, 10); - t.tm_isdst = -1; // Is DST on? 1 = yes, 0 = no, -1 = unknown - aorc->time = mktime(&t); - - // APCP_surface - value = strsep(©, ","); - // Not sure what this is - - // DLWRF_surface - value = strsep(©, ","); - aorc->incoming_longwave_W_per_m2 = strtof(value, NULL); - // DSWRF_surface - value = strsep(©, ","); - aorc->incoming_shortwave_W_per_m2 = strtof(value, NULL); - // PRES_surface - value = strsep(©, ","); - aorc->surface_pressure_Pa = strtof(value, NULL); - // SPFH_2maboveground - value = strsep(©, ","); - aorc->specific_humidity_2m_kg_per_kg = strtof(value, NULL);; - // TMP_2maboveground - value = strsep(©, ","); - aorc->air_temperature_2m_K = strtof(value, NULL); - // UGRD_10maboveground - value = strsep(©, ","); - aorc->u_wind_speed_10m_m_per_s = strtof(value, NULL); - // VGRD_10maboveground - value = strsep(©, ","); - aorc->v_wind_speed_10m_m_per_s = strtof(value, NULL); - // precip_rate - value = strsep(©, ","); - aorc->precip_kg_per_m2 = strtof(value, NULL); - - // Go ahead and free the duplicate copy now - free(copy_to_free); - + + //len=strlen(theString); + + start=0; /* begin at the beginning of theString */ + get_word_aorc(theString,&start,&end,theWord,&wordlen); + *year=atol(theWord); + + get_word_aorc(theString,&start,&end,theWord,&wordlen); + *month=atol(theWord); + + get_word_aorc(theString,&start,&end,theWord,&wordlen); + *day=atol(theWord); + + get_word_aorc(theString,&start,&end,theWord,&wordlen); + *hour=atol(theWord); + + get_word_aorc(theString,&start,&end,theWord,&wordlen); + *minute=atol(theWord); + + get_word_aorc(theString,&start,&end,theWord,&wordlen); + *second=(double)atof(theWord); + + get_word_aorc(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_aorc(theString,&start,&end,theWord,&wordlen); + aorc->incoming_longwave_W_per_m2=(double)atof(theWord); + + get_word_aorc(theString,&start,&end,theWord,&wordlen); + aorc->incoming_shortwave_W_per_m2=(double)atof(theWord); + + get_word_aorc(theString,&start,&end,theWord,&wordlen); + aorc->surface_pressure_Pa=(double)atof(theWord); + + get_word_aorc(theString,&start,&end,theWord,&wordlen); + aorc->specific_humidity_2m_kg_per_kg=(double)atof(theWord); + + get_word_aorc(theString,&start,&end,theWord,&wordlen); + aorc->air_temperature_2m_K=(double)atof(theWord); + + get_word_aorc(theString,&start,&end,theWord,&wordlen); + aorc->u_wind_speed_10m_m_per_s=(double)atof(theWord); + + get_word_aorc(theString,&start,&end,theWord,&wordlen); + 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; -} + } /*####################################################################*/ /*############################## GET WORD ############################*/ diff --git a/forcing_code/src/bmi_aorc.c b/forcing_code/src/bmi_aorc.c index 2bfc3295..059cf875 100644 --- a/forcing_code/src/bmi_aorc.c +++ b/forcing_code/src/bmi_aorc.c @@ -72,7 +72,7 @@ Initialize (Bmi *self, const char *cfg_file) for (int i = 0; i < aorc->bmi.num_timesteps; i++) { fgets(line_str, max_forcing_line_length + 1, ffp); // read in a line of AORC data. parse_aorc_line_aorc(line_str, &year, &month, &day, &hour, &minute, &dsec, &forcings); - aorc->forcing_data_precip_kg_per_m2[i] = forcings.precip_kg_per_m2 * ((double)aorc->bmi.time_step_size); + aorc->forcing_data_precip_kg_per_m2[i] = forcings.precip_kg_per_m2; //jmframe * ((double)aorc->bmi.time_step_size); if (aorc->bmi.verbose >4) printf("precip %f \n", aorc->forcing_data_precip_kg_per_m2[i]); aorc->forcing_data_surface_pressure_Pa[i] = forcings.surface_pressure_Pa; diff --git a/include/bmi.h b/include/bmi.h index 061a17a0..e6bb6d12 100755 --- a/include/bmi.h +++ b/include/bmi.h @@ -67,9 +67,7 @@ typedef struct Bmi { int (*get_output_var_names)(struct Bmi *self, char **names); - /* Variable information */ - // int (*get_var_role)(struct Bmi *self, const char *name, char* role); - + /* Variable information */ int (*get_var_grid)(struct Bmi *self, const char *name, int *grid); int (*get_var_type)(struct Bmi *self, const char *name, char *type); @@ -99,15 +97,7 @@ typedef struct Bmi { int (*get_value_ptr)(struct Bmi *self, const char *name, void **dest_ptr); int (*get_value_at_indices)(struct Bmi *self, const char *name, void *dest, int *inds, int count); - - // New BMI functions to support serialization - int (*get_state_var_count)(struct Bmi *self, int *count); - int (*get_state_var_names)(struct Bmi *self, char ** names); - int (*get_state_var_types)(struct Bmi *self, char ** types); - int (*get_state_var_ptrs)(struct Bmi *self, void *ptr_list[]); - int (*get_state_var_sizes)(struct Bmi *self, unsigned int sizes[]); - int (*set_state_var)(struct Bmi *self, void *src, int index); - + /* Setters */ int (*set_value)(struct Bmi *self, const char *name, void *src); diff --git a/include/bmi_cfe.h b/include/bmi_cfe.h index f4962344..02fad65d 100644 --- a/include/bmi_cfe.h +++ b/include/bmi_cfe.h @@ -44,25 +44,6 @@ struct aorc_forcing_data_cfe } ; typedef struct aorc_forcing_data_cfe aorc_forcing_data_cfe; -struct vol_tracking_struct{ - double vol_sch_runoff; - double vol_sch_infilt; - double vol_to_soil; - double vol_to_gw; - double vol_soil_to_gw; - double vol_soil_to_lat_flow; - double volstart; - double volout; - double volin; - double vol_from_gw; - double vol_out_giuh; - double vol_in_nash; - double vol_out_nash; - double vol_in_gw_start; - double vol_soil_start; -}; -typedef struct vol_tracking_struct vol_tracking_struct; - struct cfe_state_struct { // ************************************* @@ -80,7 +61,10 @@ struct cfe_state_struct { struct conceptual_reservoir gw_reservoir; struct NWM_soil_parameters NWM_soil_params; struct evapotranspiration_structure et_struct; - struct vol_tracking_struct vol_struct; + struct massbal vol_struct; + + /* xinanjiang_dev */ + struct direct_runoff_parameters_structure direct_runoff_params_struct; // Epoch-based start time (BMI start time is considered 0.0) long epoch_start_time; @@ -93,7 +77,9 @@ struct cfe_state_struct { char* forcing_file; - double Schaake_adjusted_magic_constant_by_soil_type; + /* xinanjiang_dev + double Schaake_adjusted_magic_constant_by_soil_type; */ + int num_lateral_flow_nash_reservoirs; double K_lf; @@ -114,9 +100,11 @@ struct cfe_state_struct { double* nash_storage; double* runoff_queue_m_per_timestep; - // These are likely only single values, but should be allocated as pointers so the pointer can be returned -// double* flux_overland_m; NOT NEEDED, redundant with flux_Schaake_output_runoff_m - double* flux_Schaake_output_runoff_m; + /* xinanjiang_dev + changing the name to the more general "direct runoff" + double* flux_Schaake_output_runoff_m;*/ + double* flux_output_direct_runoff_m ; + double* flux_giuh_runoff_m; double* flux_nash_lateral_runoff_m; double* flux_from_deep_gw_to_chan_m; diff --git a/include/cfe.h b/include/cfe.h index 8126c41c..0b867649 100644 --- a/include/cfe.h +++ b/include/cfe.h @@ -106,11 +106,57 @@ struct evapotranspiration_structure { }; typedef struct evapotranspiration_structure evapotranspiration_structure; +struct massbal +{ + double volstart ; + double vol_runoff ; + double vol_infilt ; + double vol_out_giuh ; + double vol_end_giuh ; + double vol_to_gw ; + double vol_in_gw_start ; + double vol_in_gw_end ; + double vol_from_gw ; + double vol_in_nash ; + double vol_in_nash_end ; // note the nash cascade is empty at start of simulation. + double vol_out_nash ; + double vol_soil_start ; + double vol_to_soil ; + double vol_soil_to_lat_flow; + double vol_soil_to_gw ; // this should equal vol_to_gw + double vol_soil_end ; + double volin ; + double volout ; + double volend ; +}; +typedef struct massbal massbal; + +// define data types +//-------------------------- +typedef enum {Schaake, Xinanjiang} surface_water_partition_type; + +/* xinanjiang_dev*/ +struct direct_runoff_parameters_structure{ + surface_water_partition_type surface_partitioning_scheme; + double Schaake_adjusted_magic_constant_by_soil_type; + double a_Xinanjiang_inflection_point_parameter; + double b_Xinanjiang_shape_parameter; + double x_Xinanjiang_shape_parameter; +}; +typedef struct direct_runoff_parameters_structure direct_runoff_parameters_structure; + + // function prototypes // -------------------------------- extern void Schaake_partitioning_scheme(double dt, double magic_number, double deficit, double qinsur, double *runsrf, double *pddum); +// xinanjiang_dev: XinJiang function written by Rachel adapted by Jmframe and FLO, +extern void Xinanjiang_partitioning_scheme(double water_input_depth_m, double field_capacity_m, + double max_soil_moisture_storage_m, double column_total_soil_water_m, + struct direct_runoff_parameters_structure *parms, + double *surface_runoff_depth_m, double *infiltration_depth_m); + extern void conceptual_reservoir_flux_calc(struct conceptual_reservoir *da_reservoir, double *primary_flux_m, double *secondary_flux_m); @@ -131,37 +177,36 @@ extern void cfe( struct NWM_soil_parameters NWM_soil_params_struct, struct conceptual_reservoir *soil_reservoir_struct, double timestep_h, - double Schaake_adjusted_magic_constant_by_soil_type, + + /* xinanjiang_dev: since we are doing the option for Schaake and XinJiang, + instead of passing in the constants + pass in a structure with the constants for both subroutines. + //double Schaake_adjusted_magic_constant_by_soil_type,*/ + struct direct_runoff_parameters_structure direct_runoff_param_struct, + double timestep_rainfall_input_m, - double *Schaake_output_runoff_m_ptr, + + /* xinanjiang_dev: + double *Schaake_output_runoff_m_ptr,*/ + double *flux_output_direct_runoff_m, + double *infiltration_depth_m_ptr, -// double *flux_overland_m_ptr, // NOT NEEDED redundant with Schaake_output_runoff_m_fptr - double *vol_sch_runoff_ptr, - double *vol_sch_infilt_ptr, double *flux_perc_m_ptr, - double *vol_to_soil_ptr, double *flux_lat_m_ptr, double *gw_reservoir_storage_deficit_m_ptr, struct conceptual_reservoir *gw_reservoir_struct, - double *vol_to_gw_ptr, - double *vol_soil_to_gw_ptr, - double *vol_soil_to_lat_flow_ptr, - double *volout_ptr, double *flux_from_deep_gw_to_chan_m_ptr, - double *vol_from_gw_ptr, double *giuh_runoff_m_ptr, int num_giuh_ordinates, double *giuh_ordinates_arr, double *runoff_queue_m_per_timestep_arr, - double *vol_out_giuh_ptr, double *nash_lateral_runoff_m_ptr, int num_lateral_flow_nash_reservoirs, double K_nash, double *nash_storage_arr, - double *vol_in_nash_ptr, - double *vol_out_nash_ptr, struct evapotranspiration_structure *evap_struct, - double *Qout_m_ptr + double *Qout_m_ptr, + struct massbal *massbal_struct ); #endif //CFE_CFE_H diff --git a/include/serialize_state.h b/include/serialize_state.h deleted file mode 100644 index 95113093..00000000 --- a/include/serialize_state.h +++ /dev/null @@ -1,19 +0,0 @@ - -#ifndef SERIALIZE_STATE_H -#define SERIALIZE_STATE_H - -#if defined(__cplusplus) -extern "C" { -#endif - -int serialize(Bmi* model1, const char *ser_file); - -int deserialize_to_state(const char *ser_file, Bmi* model2, int print_obj); - -int compare_states( Bmi* model1, Bmi* model2); - -#if defined(__cplusplus) -} -#endif - -#endif \ No newline at end of file diff --git a/original_author_code/CFE_1.1.c b/original_author_code/CFE_1.1.c new file mode 100644 index 00000000..c35ae399 --- /dev/null +++ b/original_author_code/CFE_1.1.c @@ -0,0 +1,1427 @@ +#include +#include +#include +#include +#include + +#define TRUE 1 +#define FALSE 0 +#define DEBUG 1 +#define MAX_NUM_GIUH_ORDINATES 10 +#define MAX_NUM_NASH_CASCADE 3 +#define MAX_NUM_RAIN_DATA 720 + +// 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 +// conceptual model. The hypothesis is supported by a number of observations made during a 2017-2018 deep dive +// into the NWM code. Thesed are: +// +// 1. Rainfall/throughfall/melt partitioning in the NWM is based on a simple curve-number like approach that +// was developed by Schaake et al. (1996) and which is very similar to the Probability Distributed Moisture (PDM) +// 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 +// 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 +// assuming no-flow lower boundary condition. Since this form of Richards' equation cannot simulate heterogeneous +// soil layers, it can be replaced with a conceptual reservoir. +// 4. The lateral flow routing function in the NWM is purely conceptual. It is activated whenever the soil water +// content in one or more of the four Richards-equation discretizations reaches the wilting point water content. +// This activation threshold is physically unrealistic, because in most soils lateral subsurface flow is not +// active until pore water pressures become positive at some point in the soil profile. Furthermore, the lateral +// flow hydraulic conductivity is assumed to be the vertical hydraulic conductivity multiplied by a calibration +// factor "LKSATFAC" which is allowed to vary between 10 and 10,000 during calibration, resulting in an anisotropy +// ratio that varies over the same range, without correlation with physiographic characteristics or other support. +// +// This code implements these assumptions using pure conceptualizations. The formulation consists of the following: +// +// 1. Rainfall is partitioned into direct runoff and soil moisture using the Schaake function. +// 2. Rainfall that becomes direct runoff is routed to the catchment outlet using a geomorphological instantanteous +// unit hydrograph (GIUH) approach, eliminating the 250 m NWM routing grid, and the incorrect use of the Schaake +// function to simulate run-on infiltration. +// 3. Water partitioned by the Schaake function to be soil moisture is placed into a conceptual linear reservoir +// that consists of two outlets that apply a minimum storage activation threshold. This activation threshold +// is identical for both outlets, and is based on an integral solution of the storage in the soil assuming +// Clapp-Hornberger parameters equal to those used in the NWM to determine that storage corresponding to a +// soil water content 0.5 m above the soil column bottom that produces a soil suction head equal to -1/3 atm, +// which is a commonly applied assumption used to estimate the field capacity water content. +// The first outlet calculates vertical percolation of water to deep groundwater using the saturated hydraulic +// conductivity of the soil multiplied by the NWM "slope" parameter, which when 1.0 indicates free drainage and +// when 0.0 indicates a no-flow lower boundary condition. The second outlet is used to calculate the flux to +// the soil lateral flow path, using a conceptual LKSATFAC-like calibration parameter. +// 4. The lateral flow is routed to the catchment outlet using a Nash-cascade of reservoirs to produce a mass- +// conserving delayed response, and elminates the need for the 250 m lateral flow routing grid. +// 5. The groundwater contribution to base flow is modeled using either (a) an exponential nonlinear reservoir +// identical to the one in the NWM formulation, or (b) a nonlinear reservoir forumulation, which can also be +// made linear by assuming an exponent value equal to 1.0. +// +// This code was written entirely by Fred L. Ogden, May 22-24, 2020, in the service of the NOAA-NWS Office of Water +// Prediction, in Tuscaloosa, Alabama. + + +// 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 storage_max_m; // maximum storage in this reservoir +double storage_m; // state variable. +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; +}; + +struct NWM_soil_parameters +{ +// using same variable names as used in NWM. +double smcmax; // effective porosity [V/V] +double wltsmc; // wilting point soil moisture content [V/V] +double satdk; // saturated hydraulic conductivity [m s-1] +double satpsi; // saturated capillary head [m] +double bb; // beta exponent on Clapp-Hornberger (1978) soil water relations [-] +double mult; // the multiplier applied to satdk to route water rapidly downslope +double slop; // this factor (0-1) modifies the gradient of the hydraulic head at the soil bottom. 0=no-flow. +double D; // soil depth [m] +}; + +struct direct_runoff_parameters_structure{ + int method; + double Schaake_adjusted_magic_constant_by_soil_type; + double a_Xinanjiang_inflection_point_parameter; + double b_Xinanjiang_shape_parameter; + double x_Xinanjiang_shape_parameter; +}; + +//DATA STRUCTURE TO HOLD AORC FORCING DATA +struct aorc_forcing_data +{ +// struct NAME DESCRIPTION ORIGINAL AORC NAME +//____________________________________________________________________________________________________________________ +float precip_kg_per_m2; // Surface precipitation "kg/m^2" | APCP_surface +float incoming_longwave_W_per_m2 ; // Downward Long-Wave Rad. Flux at 0m height, W/m^2 | DLWRF_surface +float incoming_shortwave_W_per_m2; // Downward Short-Wave Radiation Flux at 0m height, W/m^2 | DSWRF_surface +float surface_pressure_Pa; // Surface atmospheric pressure, Pa | PRES_surface +float specific_humidity_2m_kg_per_kg; // Specific Humidity at 2m height, kg/kg | SPFH_2maboveground +float air_temperature_2m_K; // Air temparture at 2m height, K | TMP_2maboveground +float u_wind_speed_10m_m_per_s; // U-component of Wind at 10m height, m/s | UGRD_10maboveground +float v_wind_speed_10m_m_per_s; // V-component of Wind at 10m height, m/s | VGRD_10maboveground +float latitude; // degrees north of the equator. Negative south | latitude +float longitude; // degrees east of prime meridian. Negative west | longitude +long int time; //TODO: type? // seconds since 1970-01-01 00:00:00.0 0:00 | time +} ; + +struct massbal +{ +double volstart ; +double vol_runoff ; +double vol_infilt ; +double vol_out_giuh ; +double vol_end_giuh ; +double vol_to_gw ; +double vol_in_gw_start ; +double vol_in_gw_end ; +double vol_from_gw ; +double vol_in_nash ; +double vol_in_nash_end ; // note the nash cascade is empty at start of simulation. +double vol_out_nash ; +double vol_soil_start ; +double vol_to_soil ; +double vol_soil_to_lat_flow; +double vol_soil_to_gw ; // this should equal vol_to_gw +double vol_soil_end ; +double volin ; +double volout ; +double volend ; +}; + +// FLO NEW +// define data types +//-------------------------- + +typedef enum {Schaake, Xinanjiang} surface_water_partition_type; + + + +// function prototypes +// -------------------------------- +extern void Schaake_partitioning_scheme(double dt, double magic_number, double deficit, double qinsur, + double *runsrf, double *pddum); + +extern void Xinanjiang_partitioning_scheme(double water_input_depth_m, double field_capacity_m, + double max_soil_moisture_storage_m, double column_total_soil_water_m, + struct direct_runoff_parameters_structure *parms, + double *surface_runoff_depth_m, double *infiltration_depth_m); + +extern void conceptual_reservoir_flux_calc(struct conceptual_reservoir *da_reservoir, + double *primary_flux,double *secondary_flux); + +extern double convolution_integral(double runoff_m, int num_giuh_ordinates, + double *giuh_ordinates, double *runoff_queue_m_per_timestep); + +extern double nash_cascade(double flux_lat_m,int num_lateral_flow_nash_reservoirs, + double K_nash,double *nash_storage); + + +extern int is_fabs_less_than_epsilon(double a,double epsilon); // returns TRUE iff fabs(a)MAX_NUM_GIUH_ORDINATES) + { + fprintf(stderr,"Big problem, number of giuh ordinates greater than MAX_NUM_GIUH_ORDINATES. Stopped.\n"); + exit(0); + } +giuh_ordinates[0]=0.06; // note these sum to 1.0. If we have N ordinates, we need a queue sized N+1 to perform +giuh_ordinates[1]=0.51; // the convolution. +giuh_ordinates[2]=0.28; +giuh_ordinates[3]=0.12; +giuh_ordinates[4]=0.03; + +// initialize Schaake parameters +//-------------------------------- +// in noah-mp refkdt=3.0. +direct_runoff_params_struct.Schaake_adjusted_magic_constant_by_soil_type=refkdt*NWM_soil_params.satdk/2.0e-06; // 2.0e-06 is used in noah-mp + + +// initialize lateral flow function parameters +//--------------------------------------------- +assumed_near_channel_water_table_slope=0.01; // [L/L] +lateral_flow_threshold_storage_m=field_capacity_storage_threshold_m; // making them the same, but they don't have 2B + +// Equation 10 in parameter equivalence document. +lateral_flow_linear_reservoir_constant=2.0*assumed_near_channel_water_table_slope*NWM_soil_params.mult* + NWM_soil_params.satdk*NWM_soil_params.D*drainage_density_km_per_km2; // m/s +lateral_flow_linear_reservoir_constant*=3600.0; // convert to m/h + +num_lateral_flow_nash_reservoirs=2; +if(num_lateral_flow_nash_reservoirs>MAX_NUM_NASH_CASCADE) + { + fprintf(stdout,"Number of Nash Cascade linear reservoirs greater than MAX_NUM_NASH_CASCADE. Stopped.\n"); + return(-2); + } +K_nash=0.03; // lateral_flow_nash_cascade_reservoir_constant. TODO Calibrate. + +d_alloc(&nash_storage,MAX_NUM_NASH_CASCADE); +for(i=0;i 1.0 +gw_reservoir.storage_threshold_primary_m=0.0; // 0.0 means no threshold applied +gw_reservoir.storage_threshold_secondary_m=0.0; // 0.0 means no threshold applied +gw_reservoir.coeff_secondary=0.0; // 0.0 means that secondary outlet is not applied +gw_reservoir.exponent_secondary=1.0; // linear + +gw_reservoir.storage_m=gw_reservoir.storage_max_m*0.5; // INITIALIZE HALF FULL. +//------------------------------------------------------------------------------------------------------------- + +massbal_struct.volstart += gw_reservoir.storage_m; // initial mass balance checks in g.w. reservoir +massbal_struct.vol_in_gw_start = gw_reservoir.storage_m; + +// Initialize the soil conceptual reservoir data structure. Indented here to highlight different purposes +//------------------------------------------------------------------------------------------------------------- +// soil conceptual reservoir first, two outlets, two thresholds, linear (exponent=1.0). +soil_reservoir.is_exponential=FALSE; // set this true TRUE to use the exponential form of the discharge equation + // this should NEVER be set to true in the soil reservoir. +soil_reservoir.storage_max_m=NWM_soil_params.smcmax*NWM_soil_params.D; + // vertical percolation parameters------------------------------------------------ + soil_reservoir.coeff_primary=NWM_soil_params.satdk*NWM_soil_params.slop*3600.0; // m per h + soil_reservoir.exponent_primary=1.0; // 1.0=linear + soil_reservoir.storage_threshold_primary_m=field_capacity_storage_threshold_m; + // lateral flow parameters -------------------------------------------------------- + soil_reservoir.coeff_secondary=0.01; // 0.0 to deactiv. else =lateral_flow_linear_reservoir_constant; // m per h + soil_reservoir.exponent_secondary=1.0; // 1.0=linear + soil_reservoir.storage_threshold_secondary_m=lateral_flow_threshold_storage_m; + +soil_reservoir.storage_m=soil_reservoir.storage_max_m*0.667; // INITIALIZE SOIL STORAGE +//------------------------------------------------------------------------------------------------------------- + +massbal_struct.volstart += soil_reservoir.storage_m; // initial mass balance checks in soil reservoir +massbal_struct.vol_soil_start = soil_reservoir.storage_m; + + +d_alloc(&rain_rate,MAX_NUM_RAIN_DATA); + +//######################################################################################################## +// read in the aorc forcing data, just save the rain rate in mm/h in an array, ignore other vars ffor now. +//######################################################################################################## +if(yes_aorc==TRUE) // reading a csv file containing all the AORC meteorological/radiation and rainfall + { + if((in_fptr=fopen("cat87_01Dec2015.csv","r"))==NULL) + {printf("Can't open input file\n");exit(0);} + + num_rain_dat=720; // hard coded number of lines in the forcing input file. + fgets(theString,512,in_fptr); // read in the header. + for(i=0;isoil_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; + 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 + flux_overland_m+=diff; // bug found by Nels. This was missing and fixes it. + } + + massbal_struct.vol_to_soil += infiltration_depth_m; + soil_reservoir.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,&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.storage_max_m-gw_reservoir.storage_m; + + // limit amount transferred to deficit iff there is insuffienct avail. storage + 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 + massbal_struct.vol_infilt-=diff; // correct overprediction of infilt. + } + + massbal_struct.vol_to_gw +=flux_perc_m; + massbal_struct.vol_soil_to_gw +=flux_perc_m; + + gw_reservoir.storage_m += flux_perc_m; + soil_reservoir.storage_m -= flux_perc_m; + soil_reservoir.storage_m -= flux_lat_m; + massbal_struct.vol_soil_to_lat_flow += flux_lat_m; //TODO add this to nash cascade as input + massbal_struct.volout=massbal_struct.volout+flux_lat_m; + + conceptual_reservoir_flux_calc(&gw_reservoir,&primary_flux,&secondary_flux); + + flux_from_deep_gw_to_chan_m=primary_flux; // m/h <<<<<<<<<< BASE FLOW FLUX >>>>>>>>> + 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"); + + + // adjust state of deep groundwater conceptual nonlinear reservoir + //----------------------------------------------------------------- + + gw_reservoir.storage_m -= flux_from_deep_gw_to_chan_m; + + + // Solve the convolution integral ffor this time step + + giuh_runoff_m = convolution_integral(direct_output_runoff_m,num_giuh_ordinates, + giuh_ordinates,runoff_queue_m_per_timestep); + massbal_struct.vol_out_giuh+=giuh_runoff_m; + + massbal_struct.volout+=giuh_runoff_m; + massbal_struct.volout+=flux_from_deep_gw_to_chan_m; + + // Route lateral flow through the Nash cascade. + nash_lateral_runoff_m = nash_cascade(flux_lat_m,num_lateral_flow_nash_reservoirs, + K_nash,nash_storage); + massbal_struct.vol_in_nash += flux_lat_m; + massbal_struct.vol_out_nash += nash_lateral_runoff_m; + +#ifdef DEBUG + fprintf(out_debug_fptr,"%d %lf %lf\n",tstep,flux_lat_m,nash_lateral_runoff_m); +#endif + + Qout_m = giuh_runoff_m + nash_lateral_runoff_m + flux_from_deep_gw_to_chan_m; + + // <<<<<<<>>>>>>>> + // These fluxs are all in units of meters per time step. Multiply them by the "catchment_area_km2" variable + // to convert them into cubic meters per time step. + + // WRITE OUTPUTS IN mm/h ffor aiding in interpretation + fprintf(out_fptr,"%16.8lf %lf %lf %lf %lf %lf %lf\n",jdate_start+(double)tstep*1.0/24.0, + timestep_rainfall_input_m*1000.0,direct_output_runoff_m*1000.0, + giuh_runoff_m*1000.0,nash_lateral_runoff_m*1000.0, flux_from_deep_gw_to_chan_m*1000.0, + Qout_m*1000.0 ); + } // END OF TIME LOOP + + +// +// PERFORM MASS BALANCE CHECKS AND WRITE RESULTS TO stderr. +//---------------------------------------------------------- + +massbal_struct.volend= soil_reservoir.storage_m+gw_reservoir.storage_m; +massbal_struct.vol_in_gw_end = gw_reservoir.storage_m; + +#ifdef DEBUG + fclose(out_debug_fptr); +#endif + +// the GIUH queue might have water in it at the end of the simulation, so sum it up. +for(i=0;i0.0) fprintf(stderr,"global pct. err: %6.4e percent of inputs\n",global_residual/massbal_struct.volin*100.0); +else fprintf(stderr,"global pct. err: %6.4e percent of initial\n",global_residual/massbal_struct.volstart*100.0); +if(!is_fabs_less_than_epsilon(global_residual,1.0e-12)) + fprintf(stderr, "WARNING: GLOBAL MASS BALANCE CHECK FAILED\n"); + +partitioning_residual=massbal_struct.volin-massbal_struct.vol_runoff-massbal_struct.vol_infilt; +fprintf(stderr," SURFACE PARTITIONING MASS BALANCE\n"); +fprintf(stderr," surface runoff: %8.4lf m\n",massbal_struct.vol_runoff); +fprintf(stderr," infiltration: %8.4lf m\n",massbal_struct.vol_infilt); +fprintf(stderr,"partitioning residual: %6.4e m\n",partitioning_residual); // should equal 0.0 +if(!is_fabs_less_than_epsilon(partitioning_residual,1.0e-12)) + fprintf(stderr,"WARNING: SURFACE PARTITIONING MASS BALANCE CHECK FAILED\n"); + +giuh_residual=massbal_struct.vol_runoff-massbal_struct.vol_out_giuh-massbal_struct.vol_end_giuh; +fprintf(stderr," GIUH MASS BALANCE\n"); +fprintf(stderr," vol. into giuh: %8.4lf m\n",massbal_struct.vol_runoff); +fprintf(stderr," vol. out giuh: %8.4lf m\n",massbal_struct.vol_out_giuh); +fprintf(stderr," vol. end giuh q: %8.4lf m\n",massbal_struct.vol_end_giuh); +fprintf(stderr," giuh residual: %6.4e m\n",giuh_residual); // should equal zero +if(!is_fabs_less_than_epsilon(giuh_residual,1.0e-12)) + fprintf(stderr,"WARNING: GIUH MASS BALANCE CHECK FAILED\n"); + +soil_residual=massbal_struct.vol_soil_start+massbal_struct.vol_infilt-massbal_struct.vol_soil_to_lat_flow-massbal_struct.vol_soil_end-massbal_struct.vol_to_gw; +fprintf(stderr," SOIL WATER CONCEPTUAL RESERVOIR MASS BALANCE\n"); +fprintf(stderr," init soil vol: %8.4lf m\n",massbal_struct.vol_soil_start); +fprintf(stderr," vol. into soil: %8.4lf m\n",massbal_struct.vol_infilt); +fprintf(stderr,"vol.soil2latflow: %8.4lf m\n",massbal_struct.vol_soil_to_lat_flow); +fprintf(stderr," vol. soil to gw: %8.4lf m\n",massbal_struct.vol_soil_to_gw); +fprintf(stderr," final vol. soil: %8.4lf m\n",massbal_struct.vol_soil_end); +fprintf(stderr,"vol. soil resid.: %6.4e m\n",soil_residual); +if(!is_fabs_less_than_epsilon(soil_residual,1.0e-12)) + fprintf(stderr,"WARNING: SOIL CONCEPTUAL RESERVOIR MASS BALANCE CHECK FAILED\n"); + +nash_residual= massbal_struct.vol_in_nash - massbal_struct.vol_out_nash - massbal_struct.vol_in_nash_end; +fprintf(stderr," NASH CASCADE CONCEPTUAL RESERVOIR MASS BALANCE\n"); +fprintf(stderr," vol. to nash: %8.4lf m\n",massbal_struct.vol_in_nash); +fprintf(stderr," vol. from nash: %8.4lf m\n",massbal_struct.vol_out_nash); +fprintf(stderr," final vol. nash: %8.4lf m\n",massbal_struct.vol_in_nash_end); +fprintf(stderr,"nash casc resid.: %6.4e m\n",nash_residual); +if(!is_fabs_less_than_epsilon(nash_residual,1.0e-12)) + fprintf(stderr,"WARNING: NASH CASCADE CONCEPTUAL RESERVOIR MASS BALANCE CHECK FAILED\n"); + + +gw_residual=massbal_struct.vol_in_gw_start+massbal_struct.vol_to_gw-massbal_struct.vol_from_gw-massbal_struct.vol_in_gw_end; +fprintf(stderr," GROUNDWATER CONCEPTUAL RESERVOIR MASS BALANCE\n"); +fprintf(stderr,"init gw. storage: %8.4lf m\n",massbal_struct.vol_in_gw_start); +fprintf(stderr," vol to gw: %8.4lf m\n",massbal_struct.vol_to_gw); +fprintf(stderr," vol from gw: %8.4lf m\n",massbal_struct.vol_from_gw); +fprintf(stderr,"final gw.storage: %8.4lf m\n",massbal_struct.vol_in_gw_end); +fprintf(stderr," 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"); + +fclose(out_fptr); + +} //<<<<<<<<<<<< END OF MAIN PROGRAM + +//############################################################## +//################### NASH CASCADE ######################### +//############################################################## +extern double nash_cascade(double flux_lat_m,int num_lateral_flow_nash_reservoirs, + double K_nash,double *nash_storage) +{ +//############################################################## +// Solve for the flow through the Nash cascade to delay the +// arrival of the lateral flow into the channel +//############################################################## +// local vars +int i; +double outflow_m; +static double Q[MAX_NUM_NASH_CASCADE]; + +//Loop through reservoirs +for(i = 0; i < num_lateral_flow_nash_reservoirs; i++) + { + Q[i] = K_nash*nash_storage[i]; + nash_storage[i] -= Q[i]; + + if (i==0) nash_storage[i] += flux_lat_m; + else nash_storage[i] += Q[i-1]; + } + +/* Get Qout */ +outflow_m = Q[num_lateral_flow_nash_reservoirs-1]; + +//Return the flow output +return (outflow_m); + +} + +//############################################################## +//############### GIUH CONVOLUTION INTEGRAL ################## +//############################################################## +extern double convolution_integral(double runoff_m,int num_giuh_ordinates, + double *giuh_ordinates, double *runoff_queue_m_per_timestep) +{ +//############################################################## +// This function solves the convolution integral involving N +// GIUH ordinates. +//############################################################## +double runoff_m_now; +int N,i; + +N=num_giuh_ordinates; +runoff_queue_m_per_timestep[N]=0.0; + +for(i=0;i>>> +//{ +// 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(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=reservoir->coeff_primary* + (exp(reservoir->exponent_primary*reservoir->storage_m/reservoir->storage_max_m)-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=reservoir->storage_m-reservoir->storage_threshold_primary_m; +if(storage_above_threshold_m>0.0) + { + // flow is possible from the primary outlet + *primary_flux_m=reservoir->coeff_primary* + pow(storage_above_threshold_m/(reservoir->storage_max_m-reservoir->storage_threshold_primary_m), + 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=reservoir->storage_m-reservoir->storage_threshold_secondary_m; +if(storage_above_threshold_m>0.0) + { + // flow is possible from the secondary outlet + *secondary_flux_m=reservoir->coeff_secondary* + pow(storage_above_threshold_m/(reservoir->storage_max_m-reservoir->storage_threshold_secondary_m), + 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 ############# +//############################################################## +void Schaake_partitioning_scheme(double timestep_h, double Schaake_adjusted_magic_constant_by_soil_type, + double column_total_soil_moisture_deficit_m, + double water_input_depth_m,double *surface_runoff_depth_m,double *infiltration_depth_m) +{ + + +/*! =============================================================================== + This subtroutine takes water_input_depth_m and partitions it into surface_runoff_depth_m and + infiltration_depth_m using the scheme from Schaake et al. 1996. +! -------------------------------------------------------------------------------- +! ! modified by FLO April 2020 to eliminate reference to ice processes, +! ! and to de-obfuscate and use descriptive and dimensionally consistent variable names. +! -------------------------------------------------------------------------------- +! inputs + double timestep_h + double Schaake_adjusted_magic_constant_by_soil_type = C*Ks(soiltype)/Ks_ref, where C=3, and Ks_ref=2.0E-06 m/s + double column_total_soil_moisture_deficit_m + double water_input_depth_m amount of water input to soil surface this time step [m] + +! outputs + double surface_runoff_depth_m amount of water partitioned to surface water this time step [m] + + +--------------------------------------------------------------------------------*/ +int k; +double timestep_d,Schaake_parenthetical_term,Ic,Px,infilt_dep_m; + + + +if (0.0 >= column_total_soil_moisture_deficit_m) + { + *surface_runoff_depth_m=water_input_depth_m; + *infiltration_depth_m=0.0; + } +else + { + // partition time-step total applied water as per Schaake et al. 1996. + // change from dt in [s] to dt1 in [d] because kdt has units of [d^(-1)] + timestep_d = timestep_h /24.0; // timestep_d is the time step in days. + + // calculate the parenthetical part of Eqn. 34 from Schaake et al. Note the magic constant has units of [d^(-1)] + + Schaake_parenthetical_term = (1.0 - exp ( - Schaake_adjusted_magic_constant_by_soil_type * timestep_d)); + + // From Schaake et al. Eqn. 2., using the column total moisture deficit + // BUT the way it is used here, it is the cumulative soil moisture deficit in the entire soil profile. + // "Layer" info not used in this subroutine in noah-mp, except to sum up the total soil moisture storage. + // NOTE: when column_total_soil_moisture_deficit_m becomes zero, which occurs when the soil column is saturated, + // then Ic=0, where Ic in the Schaake paper is called the "spatially averaged infiltration capacity", + // and is defined in Eqn. 12. + + Ic = column_total_soil_moisture_deficit_m * Schaake_parenthetical_term; + + Px=water_input_depth_m; // Total water input to partitioning scheme this time step [m] + + // This is eqn 24 from Schaake et al. NOTE: this is 0 in the case of a saturated soil column, when Ic=0. + // Physically happens only if soil has no-flow lower b.c. + + *infiltration_depth_m = (Px * (Ic / (Px + Ic))); + + + if( 0.0 < (water_input_depth_m-(*infiltration_depth_m)) ) + { + *surface_runoff_depth_m = water_input_depth_m - (*infiltration_depth_m); + } + else *surface_runoff_depth_m=0.0; + *infiltration_depth_m = water_input_depth_m - (*surface_runoff_depth_m); + } + +return; +} + +//############################################################## +//######## XINANJIANG RUNOFF PARTITIONING SCHEME ########### +//############################################################## + +void Xinanjiang_partitioning_scheme(double water_input_depth_m, double field_capacity_m, + double max_soil_moisture_storage_m, double column_total_soil_water_m, + struct direct_runoff_parameters_structure *parms, + double *surface_runoff_depth_m, double *infiltration_depth_m) +{ + //------------------------------------------------------------------------ + // This module takes the water_input_depth_m and separates it into surface_runoff_depth_m + // and infiltration_depth_m by calculating the saturated area and runoff based on a scheme developed + // for the Xinanjiang model by Jaywardena and Zhou (2000). According to Knoben et al. + // (2019) "the model uses a variable contributing area to simulate runoff. [It] uses + // a double parabolic curve to simulate tension water capacities within the catchment, + // instead of the original single parabolic curve" which is also used as the standard + // VIC fomulation. This runoff scheme was selected for implementation into NWM v3.0. + // REFERENCES: + // 1. Jaywardena, A.W. and M.C. Zhou, 2000. A modified spatial soil moisture storage + // capacity distribution curve for the Xinanjiang model. Journal of Hydrology 227: 93-113 + // 2. Knoben, W.J.M. et al., 2019. Supplement of Modular Assessment of Rainfall-Runoff Models + // Toolbox (MARRMoT) v1.2: an open-source, extendable framework providing implementations + // of 46 conceptual hydrologic models as continuous state-space formulations. Supplement of + // Geosci. Model Dev. 12: 2463-2480. + //------------------------------------------------------------------------- + // Written by RLM May 2021 + // Adapted by JMFrame September 2021 for new version of CFE + //------------------------------------------------------------------------- + // Inputs + // double water_input_depth_m amount of water input to soil surface this time step [m] + // double field_capacity_m amount of water stored in soil reservoir when at field capacity [m] + // double max_soil_moisture_storage_m total storage of the soil moisture reservoir (porosity*soil thickness) [m] + // double column_total_soil_water_m current storage of the soil moisture reservoir [m] + // double a_inflection_point_parameter a parameter + // double b_shape_parameter b parameter + // double x_shape_parameter x parameter + // + // Outputs + // double surface_runoff_depth_m amount of water partitioned to surface water this time step [m] + // double infiltration_depth_m amount of water partitioned as infiltration (soil water input) this time step [m] + //------------------------------------------------------------------------- + + double tension_water_m, free_water_m, max_tension_water_m, max_free_water_m, pervious_runoff_m; + + //could move this if statement outside of both the Schaake and Xinanjiang subroutines edit FLO- moved to main(). + + // partition the total soil water in the column between free water and tension water + free_water_m = column_total_soil_water_m - field_capacity_m; + + if(0.0 < free_water_m) { //edit FLO + tension_water_m = field_capacity_m; + } else { + free_water_m = 0.0; + tension_water_m = column_total_soil_water_m; + } + + // estimate the maximum free water and tension water available in the soil column + max_free_water_m = max_soil_moisture_storage_m - field_capacity_m; + max_tension_water_m = field_capacity_m; + + // check that the free_water_m and tension_water_m do not exceed the maximum and if so, change to the max value + if(max_free_water_m < free_water_m) free_water_m = max_free_water_m; + if(max_tension_water_m < tension_water_m) tension_water_m = max_tension_water_m; + + // NOTE: the impervious surface runoff assumptions due to frozen soil used in NWM 3.0 have not been included. + // We are assuming an impervious area due to frozen soils equal to 0 (see eq. 309 from Knoben et al). + + // The total (pervious) runoff is first estimated before partitioning into surface and subsurface components. + // See Knoben et al eq 310 for total runoff and eqs 313-315 for partitioning between surface and subsurface + // components. + + // Calculate total estimated pervious runoff. + // NOTE: If the impervious surface runoff due to frozen soils is added, + // the pervious_runoff_m equation will need to be adjusted by the fraction of pervious area. + if ((tension_water_m/max_tension_water_m) <= (0.5 - parms->a_Xinanjiang_inflection_point_parameter)) { + pervious_runoff_m = water_input_depth_m * (pow((0.5 - parms->a_Xinanjiang_inflection_point_parameter), + (1.0 - parms->b_Xinanjiang_shape_parameter)) * + pow((1.0 - (tension_water_m/max_tension_water_m)), + parms->b_Xinanjiang_shape_parameter)); + + } else { + pervious_runoff_m = water_input_depth_m * (1.0 - pow((0.5 + parms->a_Xinanjiang_inflection_point_parameter), + (1.0 - parms->b_Xinanjiang_shape_parameter)) * + pow((1.0 - (tension_water_m/max_tension_water_m)), + (parms->b_Xinanjiang_shape_parameter))); + } + // Separate the surface water from the pervious runoff + // NOTE: If impervious runoff is added to this subroutine, impervious runoff should be added to + // the surface_runoff_depth_m. + *surface_runoff_depth_m = pervious_runoff_m * (1.0 - pow((1.0 - (free_water_m/max_free_water_m)),parms->x_Xinanjiang_shape_parameter)); + + // The surface runoff depth is bounded by a minimum of 0 and a maximum of the water input depth. + // Check that the estimated surface runoff is not less than 0.0 and if so, change the value to 0.0. + if(*surface_runoff_depth_m < 0.0) *surface_runoff_depth_m = 0.0; + // Check that the estimated surface runoff does not exceed the amount of water input to the soil surface. If it does, + // change the surface water runoff value to the water input depth. + if(*surface_runoff_depth_m > water_input_depth_m) *surface_runoff_depth_m = water_input_depth_m; + // Separate the infiltration from the total water input depth to the soil surface. + *infiltration_depth_m = water_input_depth_m - *surface_runoff_depth_m; + +return; +} + + + +extern int is_fabs_less_than_epsilon(double a,double epsilon) // returns true if fabs(a)precip_kg_per_m2=atof(theWord); + +get_word(theString,&start,&end,theWord,&wordlen); +aorc->incoming_longwave_W_per_m2=atof(theWord); + +get_word(theString,&start,&end,theWord,&wordlen); +aorc->incoming_shortwave_W_per_m2=atof(theWord); + +get_word(theString,&start,&end,theWord,&wordlen); +aorc->surface_pressure_Pa=atof(theWord); + +get_word(theString,&start,&end,theWord,&wordlen); +aorc->specific_humidity_2m_kg_per_kg=atof(theWord); + +get_word(theString,&start,&end,theWord,&wordlen); +aorc->air_temperature_2m_K=atof(theWord); + +get_word(theString,&start,&end,theWord,&wordlen); +aorc->u_wind_speed_10m_m_per_s=atof(theWord); + +get_word(theString,&start,&end,theWord,&wordlen); +aorc->v_wind_speed_10m_m_per_s=atof(theWord); + + +return; +} + +/*####################################################################*/ +/*############################## GET WORD ############################*/ +/*####################################################################*/ +void get_word(char *theString,int *start,int *end,char *theWord,int *wordlen) +{ +int i,lenny,j; +lenny=strlen(theString); + +while(theString[*start]==' ' || theString[*start]=='\t') + { + (*start)++; + }; + +j=0; +for(i=(*start);i 2) + m -= 3; + else { + m += 9; + --y; + } + c = y / 100L; + ya = y - (100L * c); + j = (146097L * c) / 4L + (1461L * ya) / 4L + (153L * m + 2L) / 5L + d + 1721119L; + if (seconds < 12 * 3600.0) { + j--; + seconds += 12.0 * 3600.0; + } + else { + seconds = seconds - 12.0 * 3600.0; + } + return (j + (seconds / 3600.0) / 24.0); +} + +/* Julian date converter. Takes a julian date (the number of days since +** some distant epoch or other), and returns an int pointer to static space. +** ip[0] = month; +** ip[1] = day of month; +** ip[2] = year (actual year, like 1977, not 77 unless it was 77 a.d.); +** ip[3] = day of week (0->Sunday to 6->Saturday) +** These are Gregorian. +** 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 +** +*/ +void calc_date(double jd, long *y, long *m, long *d, long *h, long *mi, + double *sec) +{ + static int ret[4]; + + long j; + double tmp; + double frac; + + j=(long)jd; + frac = jd - j; + + if (frac >= 0.5) { + frac = frac - 0.5; + j++; + } + else { + frac = frac + 0.5; + } + + ret[3] = (j + 1L) % 7L; + j -= 1721119L; + *y = (4L * j - 1L) / 146097L; + j = 4L * j - 1L - 146097L * *y; + *d = j / 4L; + j = (4L * *d + 3L) / 1461L; + *d = 4L * *d + 3L - 1461L * j; + *d = (*d + 4L) / 4L; + *m = (5L * *d - 3L) / 153L; + *d = 5L * *d - 3 - 153L * *m; + *d = (*d + 5L) / 5L; + *y = 100L * *y + j; + if (*m < 10) + *m += 3; + else { + *m -= 9; + *y=*y+1; /* Invalid use: *y++. Modified by Tony */ + } + + /* if (*m < 3) *y++; */ + /* incorrectly repeated the above if-else statement. Deleted by Tony.*/ + + tmp = 3600.0 * (frac * 24.0); + *h = (long) (tmp / 3600.0); + tmp = tmp - *h * 3600.0; + *mi = (long) (tmp / 60.0); + *sec = tmp - *mi * 60.0; +} + +int dayofweek(double j) +{ + j += 0.5; + return (int) (j + 1) % 7; +} + diff --git a/original_author_code/make_and_run.sh b/original_author_code/make_and_run.sh index 9eec47bc..94e9dabd 100755 --- a/original_author_code/make_and_run.sh +++ b/original_author_code/make_and_run.sh @@ -1,5 +1,7 @@ #!/bin/bash gcc -lm cfe_driver.c cfe.c -o run_cfe gcc -lm t-shirt_0.99f.c -o run_cfe_f +gcc -lm CFE_1.1.c -o run_cfe_1.1 ./run_cfe ./run_cfe_f +./run_cfe_1.1 diff --git a/run_bmi_unit_test.py b/run_bmi_unit_test.py deleted file mode 100644 index d7052f5e..00000000 --- a/run_bmi_unit_test.py +++ /dev/null @@ -1,419 +0,0 @@ -"""Run BMI Unit Testing. -Author: jgarrett -Date: 08/31/2021""" - -import os -import sys -import numpy as np - -#import torch -#from torch import nn -from pathlib import Path -from netCDF4 import Dataset -import bmi_lstm # This is the BMI LSTM that we will be running - - -# setup a "success counter" for number of passing and failing bmi functions -# keep track of function def fails (vs function call) -pass_count = 0 -fail_count = 0 -var_name_counter = 0 -fail_list = [] - -def bmi_except(fstring): - """Prints message and updates counter and list - - Parameters - ---------- - fstring : str - Name of failing BMI function - """ - - global fail_count, fail_list, var_name_counter - print("**BMI ERROR** in " + fstring) - if (var_name_counter == 0): - fail_count += 1 - fail_list.append(fstring) - -bmi=bmi_lstm.bmi_LSTM() - -print("\nBEGIN BMI UNIT TEST\n*******************\n"); - -# Define config path -cfg_file=Path('./bmi_config_files/01022500_hourly_all_attributes_forcings.yml') -#cfg_file=Path('./bmi_config_files/01022500_hourly_slope_mean_precip_temp.yml') - -if os.path.exists(cfg_file): - print(" configuration found: " + str(cfg_file)) -else: - print(" no configuration found, exiting...") - sys.exit() - -#------------------------------------------------------------------- -# initialize() -try: - bmi.initialize(cfg_file) - print(" initializing..."); - pass_count += 1 -except: - bmi_except('initialize()') - -#------------------------------------------------------------------- -#------------------------------------------------------------------- -# BMI: Model Information Functions -#------------------------------------------------------------------- -#------------------------------------------------------------------- -print("\nMODEL INFORMATION\n*****************") - -#------------------------------------------------------------------- -# get_component_name() -try: - print (" component name: " + bmi.get_component_name()) - pass_count += 1 -except: - bmi_except('get_component_name()') - -#------------------------------------------------------------------- -# get_input_item_count() -try: - print (" input item count: " + str(bmi.get_input_item_count())) - pass_count += 1 -except: - bmi_except('get_input_item_count()') - -#------------------------------------------------------------------- -# get_output_item_count() -try: - print (" output item count: " + str(bmi.get_output_item_count())) - pass_count += 1 -except: - bmi_except('get_output_item_count()') - -#------------------------------------------------------------------- -# get_input_var_names() -try: - # only print statement if names exist - test_get_input_var_names = bmi.get_input_var_names() - if len(test_get_input_var_names) > 0: - print (" input var names: ") - for var_in in test_get_input_var_names: - print (" " + var_in) - pass_count += 1 -except: - bmi_except('get_input_var_names()') - -#------------------------------------------------------------------- -# get_input_var_names() -try: - # only print statement if out var list not null - test_get_output_var_names = bmi.get_output_var_names() - if len(test_get_output_var_names) > 0: - print (" output var names: ") - for var_out in test_get_output_var_names: - print (" " + var_out) - pass_count += 1 -except: - bmi_except('get_output_item_count()') - - -#------------------------------------------------------------------- -#------------------------------------------------------------------- -# BMI: Variable Information Functions -#------------------------------------------------------------------- -#------------------------------------------------------------------- -print("\nVARIABLE INFORMATION\n********************") - -for var_name in (bmi.get_output_var_names() + bmi.get_input_var_names()): - print (" " + var_name + ":") - - #------------------------------------------------------------------- - # get_var_units() - try: - print (" units: " + bmi.get_var_units(var_name)) - if var_name_counter == 0: - pass_count += 1 - except: - bmi_except('get_var_units()') - - #------------------------------------------------------------------- - # get_var_itemsize() - try: - print (" itemsize: " + str(bmi.get_var_itemsize(var_name))) - if var_name_counter == 0: - pass_count += 1 - except: - bmi_except('get_var_itemsize()') - - #------------------------------------------------------------------- - # get_var_type() - try: - print (" type: " + str(bmi.get_var_type(var_name))) - if var_name_counter == 0: - pass_count += 1 - except: - bmi_except('get_var_type()') - - #------------------------------------------------------------------- - # get_var_nbytes() - try: - print (" nbytes: " + str(bmi.get_var_nbytes(var_name))) - if var_name_counter == 0: - pass_count += 1 - except: - bmi_except('get_var_nbytes()') - - #------------------------------------------------------------------- - # get_var_grid - try: - print (" grid id: " + str(bmi.get_var_grid(var_name))) - if var_name_counter == 0: - pass_count += 1 - except: - bmi_except('get_var_grid()') - - #------------------------------------------------------------------- - # get_var_location - try: - print (" location: " + bmi.get_var_location(var_name)) - if var_name_counter == 0: - pass_count += 1 - except: - bmi_except('get_var_location()') - - var_name_counter += 1 - -# reset back to zero -var_name_counter = 0 - -#------------------------------------------------------------------- -#------------------------------------------------------------------- -# BMI: Time Functions -#------------------------------------------------------------------- -#------------------------------------------------------------------- -print("\nTIME INFORMATION\n****************") - -#------------------------------------------------------------------- -# get_start_time() -try: - print (" start time: " + str(bmi.get_start_time())) - pass_count += 1 -except: - bmi_except('get_start_time()') - -#------------------------------------------------------------------- -# get_end_time() -try: - print (" end time: " + str(bmi.get_end_time())) - pass_count += 1 -except: - bmi_except('get_end_time()') - -#------------------------------------------------------------------- -# get_current_time() -try: - print (" current time: " + str(bmi.get_current_time())) - pass_count += 1 -except: - bmi_except('get_current_time()') - -#------------------------------------------------------------------- -# get_time_step() -try: - print (" time step: " + str(bmi.get_time_step())) - pass_count += 1 -except: - bmi_except('get_time_step()') - -#------------------------------------------------------------------- -# get_time_units() -try: - print (" time units: " + bmi.get_time_units()) - pass_count += 1 -except: - bmi_except('get_time_units()') - - -#------------------------------------------------------------------- -#------------------------------------------------------------------- -# BMI: Model Grid Functions -#------------------------------------------------------------------- -#------------------------------------------------------------------- -print("\nGRID INFORMATION\n****************") -grid_id = 0 # there is only 1 -print (" grid id: " + str(grid_id)) - -#------------------------------------------------------------------- -# get_grid_rank() -try: - print (" rank: " + str(bmi.get_grid_rank(grid_id))) - pass_count += 1 -except: - bmi_except('get_grid_rank()') - -#------------------------------------------------------------------- -# get_grid_size() -try: - print (" size: " + str(bmi.get_grid_size(grid_id))) - pass_count += 1 -except: - bmi_except('get_grid_size()') - -#------------------------------------------------------------------- -# get_grid_type() -try: - print (" type: " + bmi.get_grid_type(grid_id)) - pass_count += 1 -except: - bmi_except('get_grid_type()') - - -#------------------------------------------------------------------- -#------------------------------------------------------------------- -# BMI: Variable Getter and Setter Functions -#------------------------------------------------------------------- -#------------------------------------------------------------------- -print ("\nGET AND SET VALUES\n******************") - -for var_name in (bmi.get_output_var_names() + bmi.get_input_var_names()): - print (" " + var_name + ":" ) - - #------------------------------------------------------------------- - # set_value() - try: - this_set_value = -99.0 - bmi.set_value(var_name, this_set_value) - print (" set value: " + str(this_set_value)) - - if var_name_counter == 0: - pass_count += 1 - except: - bmi_except('set_value()') - - #------------------------------------------------------------------- - # get_value() - try: - that_get_value = bmi.get_value(var_name) - # check if set_value() passed then see if get/set values match - if 'set_value()' not in fail_list: - if that_get_value == this_set_value: - print (" get value: " + str(that_get_value) + " (values match)") - else: - print (" get value: " + str(that_get_value) + " (values DO NOT match)") - else: - print (" get value: " + str(that_get_value)) - if var_name_counter == 0: - pass_count += 1 - except: - bmi_except('get_value()') - - #------------------------------------------------------------------- - # get_value_ptr() - try: - that_get_value_ptr = bmi.get_value_ptr(var_name) - # check if set_value() passed then see if get/set values match - if 'set_value()' not in fail_list: - if that_get_value_ptr == this_set_value: - print (" get value ptr: " + str(that_get_value) + " (values match)") - else: - print (" get value ptr: " + str(that_get_value) + " (values DO NOT match)") - else: - print (" get value ptr: " + str(that_get_value_ptr)) - if var_name_counter == 0: - pass_count += 1 - except: - bmi_except('get_value_ptr()') - - #------------------------------------------------------------------- - # set_value_at_indices() - try: - this_set_value_at_indices = -11.0 - bmi.set_value_at_indices(var_name,[0], this_set_value_at_indices) - #print (" set value at indices: -9.0, and got value:", bmi.get_value(var_name)) - print (" set value at indices: " + str(this_set_value_at_indices)) - if var_name_counter == 0: - pass_count += 1 - except: - bmi_except('set_value_at_indices()') - - #------------------------------------------------------------------- - # get_value_at_indices() - try: - dest0 = np.empty(bmi.get_grid_size(0), dtype=float) - that_get_value_at_indices = bmi.get_value_at_indices(var_name, dest0, [0]) - # check if set_value_at_indices() passed then see if get/set values match - if 'set_value_at_indices()' not in fail_list: - if that_get_value_at_indices == this_set_value_at_indices: - print (" get value at indices: " + str(that_get_value_at_indices) + " (values match)") - else: - print (" get value at indices: " + str(that_get_value_at_indices) + " (values DO NOT match)") - # prob worth while to bounce against set_value() if get_value_at_indices() failed.. - elif 'set_value()' not in fail_list: - if that_get_value_at_indices == this_set_value: - print (" get value at indices: " + str(that_get_value_at_indices) + " (values match)") - else: - print (" get value at indices: " + str(that_get_value_at_indices) + " (values DO NOT match)") - else: - # JMFrame NOTE: converting a list/array to a string probably won't work - #print (" get value at indices: " + str(bmi.get_value_at_indices(var_name, dest0, [0]))) - - print (" get value at indices: ", that_get_value_at_indices) - - if var_name_counter == 0: - pass_count += 1 - except: - bmi_except('get_value_at_indices()') - - var_name_counter += 1 - -# set back to zero -var_name_counter = 0 - -#------------------------------------------------------------------- -#------------------------------------------------------------------- -# BMI: Control Functions -#------------------------------------------------------------------- -#------------------------------------------------------------------- -print ("\nCONTROL FUNCTIONS\n*****************") - -#------------------------------------------------------------------- -# update() -try: - bmi.update() - # go ahead and print time to show iteration - # wrap another try/except incase get_current_time() failed - try: - print (" updating... time " + str(bmi.get_current_time())); - except: - print (" updating..."); - pass_count += 1 -except: - bmi_except('update()') - -#------------------------------------------------------------------- -# update_until() -try: - bmi.update_until(100) - # go ahead and print time to show iteration - # wrap another try/except incase get_current_time() failed - try: - print (" updating until... time " + str(bmi.get_current_time())); - except: - print (" updating until..."); - pass_count += 1 -except: - bmi_except('update_until()') - -#------------------------------------------------------------------- -# finalize() -try: - bmi.finalize() - print (" finalizing...") - pass_count += 1 -except: - bmi_except('finalize()') - -# lastly - print test summary -print ("\n Total BMI function PASS: " + str(pass_count)) -print (" Total BMI function FAIL: " + str(fail_count)) -for ff in fail_list: - print (" " + ff) \ No newline at end of file diff --git a/src/bmi_cfe.c b/src/bmi_cfe.c index 0171e161..f891303d 100755 --- a/src/bmi_cfe.c +++ b/src/bmi_cfe.c @@ -13,7 +13,7 @@ #define INPUT_VAR_NAME_COUNT 2 #define OUTPUT_VAR_NAME_COUNT 6 -#define STATE_VAR_NAME_COUNT 85 // must match var_info array size +#define STATE_VAR_NAME_COUNT 89 // must match var_info array size //---------------------------------------------- // Put variable info into a struct to simplify @@ -99,7 +99,8 @@ Variable var_info[] = { { 53, "time_step_size", "int", 1 }, { 54, "is_forcing_from_bmi", "int", 1 }, { 55, "forcing_file", "string", 1 }, // strlen - { 56, "Schaake_adjusted_magic_constant_by_soil_type", "double", 1 }, + //{ 56, "Schaake_adjusted_magic_constant_by_soil_type", "double", 1 }, + { 56, "surface_partitioning_scheme", "int", 1 }, // from direct_runoff_params_struct { 57, "num_lateral_flow_nash_reservoirs", "int", 1 }, { 58, "K_lf", "double", 1 }, { 59, "K_nash", "double", 1 }, @@ -134,7 +135,16 @@ Variable var_info[] = { { 81, "flux_perc_m", "double*", 1 }, { 82, "flux_lat_m", "double*", 1 }, { 83, "flux_Qout_m", "double*", 1 }, - { 84, "verbosity", "int", 1 } + { 84, "verbosity", "int", 1 }, + //--------------------------------------- + // direct_runoff_params_struct vars + // xinanjiang or schaake flag [56] + //--------------------------------------- + { 85, "Schaake_adjusted_magic_constant_by_soil_type", "double", 1}, + { 86, "a_Xinanjiang_inflection_point_parameter", "double", 1}, + { 87, "b_Xinanjiang_shape_parameter", "double", 1}, + { 88, "x_Xinanjiang_shape_parameter", "double", 1} + //--------------------------------------- }; int i = 0; @@ -143,7 +153,9 @@ int j = 0; // Don't forget to update Get_value/Get_value_at_indices (and setter) implementation if these are adjusted static const char *output_var_names[OUTPUT_VAR_NAME_COUNT] = { "RAIN_RATE", - "SCHAAKE_OUTPUT_RUNOFF", + /* xinanjiang_dev + "SCHAAKE_OUTPUT_RUNOFF",*/ + "DIRECT_RUNOFF", "GIUH_RUNOFF", "NASH_LATERAL_RUNOFF", "DEEP_GW_TO_CHANNEL_FLUX", @@ -245,6 +257,9 @@ static int Get_end_time (Bmi *self, double * time) 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 + // BMI shouldn't always have a very large end time... if (cfe->is_forcing_from_bmi == TRUE){ // if BMI, set to FLT_MAX macro via float.h // See https://bmi.readthedocs.io/en/latest/#get-end-time @@ -370,6 +385,12 @@ int read_init_config_cfe(const char* config_file, cfe_state_struct* model, doubl int is_num_timesteps_set = FALSE; int is_verbosity_set = FALSE; + /* xinanjiang_dev*/ + int is_direct_runoff_method_set = FALSE; + int is_a_Xinanjiang_inflection_point_parameter_set = FALSE; + int is_b_Xinanjiang_shape_parameter_set = FALSE; + int is_x_Xinanjiang_shape_parameter_set = FALSE; + // Keep track these in particular, because the "true" storage value may be a ratio and need both storage and max int is_gw_max_set = FALSE; int is_gw_storage_set = FALSE; @@ -540,6 +561,29 @@ int read_init_config_cfe(const char* config_file, cfe_state_struct* model, doubl is_verbosity_set = TRUE; continue; } + /* xinanjiang_dev: Need the option to run either runoff method in the config file, + *////////////////////////////////////////////////////////////////////////////// + if (strcmp(param_key, "surface_partitioning_scheme") == 0) { + if (strcmp(param_value, "Schaake")==0 || strcmp(param_value, "schaake")==0 || strcmp(param_value,"1")==0 ) + model->direct_runoff_params_struct.surface_partitioning_scheme = Schaake; + if (strcmp(param_value, "Xinanjiang")==0 || strcmp(param_value, "xinanjiang")==0 || strcmp(param_value,"2")==0) + model->direct_runoff_params_struct.surface_partitioning_scheme = Xinanjiang; + is_direct_runoff_method_set = TRUE; + } + if (model->direct_runoff_params_struct.surface_partitioning_scheme == Xinanjiang) { //Check that logical statement is correct + if (strcmp(param_key, "a_Xinanjiang_inflection_point_parameter") == 0){ + model->direct_runoff_params_struct.a_Xinanjiang_inflection_point_parameter = strtod(param_value, NULL); + is_a_Xinanjiang_inflection_point_parameter_set = TRUE; + } + if (strcmp(param_key, "b_Xinanjiang_shape_parameter") == 0) { + model->direct_runoff_params_struct.b_Xinanjiang_shape_parameter = strtod(param_value, NULL); + is_b_Xinanjiang_shape_parameter_set = TRUE; + } + if (strcmp(param_key, "x_Xinanjiang_shape_parameter") == 0) { + model->direct_runoff_params_struct.x_Xinanjiang_shape_parameter = strtod(param_value, NULL); + is_x_Xinanjiang_shape_parameter_set = TRUE; + } + } } if (is_forcing_file_set == FALSE) { @@ -656,17 +700,38 @@ int read_init_config_cfe(const char* config_file, cfe_state_struct* model, doubl model->verbosity = 10; return BMI_FAILURE; } - +/* xinanjiang_dev*/ + if(model->direct_runoff_params_struct.surface_partitioning_scheme == Xinanjiang){ + if (is_a_Xinanjiang_inflection_point_parameter_set == FALSE) { #if CFE_DEGUG >= 1 - printf("All CFE config params present\n"); + printf("Config param 'a_Xinanjiang_inflection_point_parameter' not found in config file\n"); #endif + return BMI_FAILURE; + } + if (is_b_Xinanjiang_shape_parameter_set == FALSE) { +#if CFE_DEGUG >= 1 + printf("Config param 'b_Xinanjiang_shape_parameter' not found in config file\n"); +#endif + return BMI_FAILURE; + } + if (is_x_Xinanjiang_shape_parameter_set == FALSE) { +#if CFE_DEGUG >= 1 + printf("Config param 'x_Xinanjiang_shape_parameter' not found in config file\n"); +#endif + return BMI_FAILURE; + } + } - model->Schaake_adjusted_magic_constant_by_soil_type = refkdt * model->NWM_soil_params.satdk / 0.000002; - + if(model->direct_runoff_params_struct.surface_partitioning_scheme == Schaake){ + model->direct_runoff_params_struct.Schaake_adjusted_magic_constant_by_soil_type = refkdt * model->NWM_soil_params.satdk / 0.000002; #if CFE_DEGUG >= 1 printf("Schaake Magic Constant calculated\n"); #endif + } +#if CFE_DEGUG >= 1 + printf("All CFE config params present\n"); +#endif // Used for parsing strings representing arrays of values below char *copy, *value; @@ -767,8 +832,12 @@ static int Initialize (Bmi *self, const char *file) JMFRAME: Moved these up before the read forcing line, Since we need them even if we don't read forcings from file. ************************************************************************/ -// cfe_bmi_data_ptr->flux_overland_m = malloc(sizeof(double)); //NOT NEEDED, redundant with flux_Schaake_output_runoff_m - cfe_bmi_data_ptr->flux_Schaake_output_runoff_m = malloc(sizeof(double)); + + /* xinanjiang_dev + changing the name to the more general "direct runoff" + cfe_bmi_data_ptr->flux_Schaake_output_runoff_m = malloc(sizeof(double));*/ + cfe_bmi_data_ptr->flux_output_direct_runoff_m = malloc(sizeof(double)); + cfe_bmi_data_ptr->flux_Qout_m = malloc(sizeof(double)); cfe_bmi_data_ptr->flux_from_deep_gw_to_chan_m = malloc(sizeof(double)); cfe_bmi_data_ptr->flux_giuh_runoff_m = malloc(sizeof(double)); @@ -783,13 +852,13 @@ static int Initialize (Bmi *self, const char *file) 2. Get the forcings passed in through BMI *******************************************************/ if (strcmp(cfe_bmi_data_ptr->forcing_file, "BMI") == 0){ - cfe_bmi_data_ptr->is_forcing_from_bmi = 1; + cfe_bmi_data_ptr->is_forcing_from_bmi = TRUE; cfe_bmi_data_ptr->forcing_data_precip_kg_per_m2 = malloc(sizeof(double)); cfe_bmi_data_ptr->forcing_data_time = malloc(sizeof(long)); } else { - cfe_bmi_data_ptr->is_forcing_from_bmi = 0; + 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; @@ -895,8 +964,8 @@ static int Update (Bmi *self) cfe_state_struct* cfe_ptr = ((cfe_state_struct *) self->data); - // Two modes to get forcing data... 1) read from file, 2) pass with bmi - if (cfe_ptr->is_forcing_from_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 cfe_ptr->timestep_rainfall_input_m = cfe_ptr->aorc.precip_kg_per_m2 /1000; @@ -906,9 +975,9 @@ static int Update (Bmi *self) cfe_ptr->timestep_rainfall_input_m = cfe_ptr->forcing_data_precip_kg_per_m2[cfe_ptr->current_time_step] /1000; cfe_ptr->vol_struct.volin += cfe_ptr->timestep_rainfall_input_m; - // Delete me... printf("_______PRECIP IN CFE UPDATE FUNCTION________: %lf\n", cfe_ptr->aorc.precip_kg_per_m2); + run_cfe(cfe_ptr); - + // Advance the model time cfe_ptr->current_time_step += 1; @@ -918,88 +987,62 @@ static int Update (Bmi *self) static int Update_until (Bmi *self, double t) { - // Since this model's time units are seconds, it is assumed that the param is either a valid time in seconds, a - // relative number of time steps into the future, or invalid + // https://bmi.readthedocs.io/en/latest/#update-until + // "the time argument can be a non-integral multiple of time steps" - // Don't support negative parameter values - if (t < 0.0) - return BMI_FAILURE; + cfe_state_struct* cfe_ptr = ((cfe_state_struct *) self->data); + + double dt; + double now; - // Don't continue if current time is at or beyond end time (or we can't determine this) - double current_time, end_time; - int current_time_result = self->get_current_time(self, ¤t_time); - if (current_time_result == BMI_FAILURE) - return BMI_FAILURE; - int end_time_result = self->get_end_time(self, &end_time); - if (end_time_result == BMI_FAILURE || current_time >= end_time) + if(self->get_time_step (self, &dt) == BMI_FAILURE) return BMI_FAILURE; - // Handle easy case of t == current_time by just returning success - if (t == current_time) - return BMI_SUCCESS; + if(self->get_current_time(self, &now) == BMI_FAILURE) + return BMI_FAILURE; - cfe_state_struct* cfe_ptr = ((cfe_state_struct *) self->data); - - // First, determine if t is some future time that will be arrived at exactly after some number of future time steps - int is_exact_future_time = (t == end_time) ? TRUE : FALSE; - // Compare to time step endings unless obvious that t lines up (i.e., t == end_time) or doesn't (t <= current_time) - if (is_exact_future_time == FALSE && t > current_time) { - int future_time_step = cfe_ptr->current_time_step; - double future_time_step_time = current_time; - while (future_time_step < cfe_ptr->num_timesteps && future_time_step_time < end_time) { - future_time_step_time += cfe_ptr->time_step_size; - if (future_time_step_time == t) { - is_exact_future_time = TRUE; - break; - } - } + { + + int n; + double frac; + const double n_steps = (t - now) / dt; + for (n=0; n<(int)n_steps; n++) { + Update (self); } - // If it is an exact time, advance to that time step - if (is_exact_future_time == TRUE) { - while (current_time < t) { - - run_cfe(cfe_ptr); - - // Advance the model time - cfe_ptr->current_time_step += 1; + frac = n_steps - (int)n_steps; + if (frac > 0){ + printf("WARNING: CFE trying to update a fraction of a timestep\n"); - // Set the current rainfall input to the right place in the forcing. - // divide by 1000 to convert from mm/h to m w/ 1h timestep as per t-shirt_0.99f - cfe_ptr->timestep_rainfall_input_m = cfe_ptr->forcing_data_precip_kg_per_m2[cfe_ptr->current_time_step] /1000; - - self->get_current_time(self, ¤t_time); - - } - return BMI_SUCCESS; + // change timestep to remaining fraction & call update() + cfe_ptr->time_step_size = frac * dt; + Update (self); + // set back to original + cfe_ptr->time_step_size = dt; } - - // If t is not an exact time, it could be a number of time step forward to proceed - - // The model doesn't support partial time step value args (i.e., fractions) - int t_int = (int) t; - if ((t - ((double)t_int)) != 0) - return BMI_FAILURE; - - // Keep in mind the current_time_step hasn't been processed yet (hence, using <= for this test) - // E.g., if (unprocessed) current_time_step = 0, t = 2, num_timesteps = 2, this is valid a valid t (run 0, run 1) - if ((cfe_ptr->current_time_step + t_int) <= cfe_ptr->num_timesteps) { - for (i = 0; i < t_int; i++){ - - // Set the current rainfall input to the right place in the forcing. - // convert from mm/h to m w/ 1h timestep as per t-shirt_0.99f - cfe_ptr->timestep_rainfall_input_m = cfe_ptr->forcing_data_precip_kg_per_m2[i] /1000; - - run_cfe(cfe_ptr); - } - - return BMI_SUCCESS; + } - // If we arrive here, t wasn't an exact time at end of a time step or a valid relative time step jump, so invalid. - return BMI_FAILURE; + return BMI_SUCCESS; } +// cfe_state_struct* cfe_ptr = ((cfe_state_struct *) self->data); +// +// int t_int = (int) t; +// if ((t - ((double)t_int)) != 0) +// return BMI_FAILURE; +// +// for (int j = 0; j < t_int; j++){ +// +// self->update(self); +// if (cfe_ptr->verbosity > 1) +// print_cfe_flux_at_timestep(cfe_ptr); +// +// } +// return BMI_SUCCESS; +//} + + static int Finalize (Bmi *self) { // Function assumes everything that is needed is retrieved from the model before Finalize is called. @@ -1016,13 +1059,15 @@ static int Finalize (Bmi *self) free(model->nash_storage); if( model->runoff_queue_m_per_timestep != NULL ) free(model->runoff_queue_m_per_timestep); - - // if( model->flux_overland_m != NULL ) //NOT NEEDED redundant with flux_Schaake_runoff_m - // free(model->flux_overland_m); 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); + free(model->flux_Schaake_output_runoff_m);*/ + if( model->flux_output_direct_runoff_m != NULL ) + free(model->flux_output_direct_runoff_m); + if( model->flux_from_deep_gw_to_chan_m != NULL ) free(model->flux_from_deep_gw_to_chan_m); if( model->flux_giuh_runoff_m != NULL ) @@ -1047,13 +1092,15 @@ static int Get_adjusted_index_for_variable(const char *name) // associated names array, plus either: // 0 if it is in the output variable array or // OUTPUT_VAR_NAME_COUNT if it is in the input variable array - for (i = 0; i < OUTPUT_VAR_NAME_COUNT; i++) + for (i = 0; i < OUTPUT_VAR_NAME_COUNT; i++){ if (strcmp(name, output_var_names[i]) == 0) return i; + } - for (i = 0; i < INPUT_VAR_NAME_COUNT; i++) + for (i = 0; i < INPUT_VAR_NAME_COUNT; i++){ if (strcmp(name, input_var_names[i]) == 0) return i + OUTPUT_VAR_NAME_COUNT; + } return -1; } @@ -1226,9 +1273,15 @@ static int Get_value_ptr (Bmi *self, const char *name, void **dest) return BMI_SUCCESS; } + /* xinanjiang_dev + changing the name to the more general "direct runoff" if (strcmp (name, "SCHAAKE_OUTPUT_RUNOFF") == 0) { *dest = (void*) ((cfe_state_struct *)(self->data))->flux_Schaake_output_runoff_m; return BMI_SUCCESS; + }*/ + if (strcmp (name, "DIRECT_RUNOFF") == 0) { + *dest = (void*) ((cfe_state_struct *)(self->data))->flux_output_direct_runoff_m; + return BMI_SUCCESS; } if (strcmp (name, "GIUH_RUNOFF") == 0) { @@ -1616,8 +1669,8 @@ static int Get_state_var_ptrs (Bmi *self, void *ptr_list[]) //------------------------------ // Vars in vol_tracking_struct //------------------------------ - ptr_list[35] = &(state->vol_struct.vol_sch_runoff ); - ptr_list[36] = &(state->vol_struct.vol_sch_infilt ); + ptr_list[35] = &(state->vol_struct.vol_runoff ); + ptr_list[36] = &(state->vol_struct.vol_infilt ); ptr_list[37] = &(state->vol_struct.vol_to_soil ); ptr_list[38] = &(state->vol_struct.vol_to_gw ); ptr_list[39] = &(state->vol_struct.vol_soil_to_gw ); @@ -1640,8 +1693,9 @@ static int Get_state_var_ptrs (Bmi *self, void *ptr_list[]) ptr_list[53] = &(state->time_step_size ); ptr_list[54] = &(state->is_forcing_from_bmi ); ptr_list[55] = state->forcing_file; - // ####### ptr_list[55] = &(state->forcing_file ); - ptr_list[56] = &(state->Schaake_adjusted_magic_constant_by_soil_type ); + // ####### ptr_list[55] = &(state->forcing_file ); + ptr_list[56] = &(state->direct_runoff_params_struct.surface_partitioning_scheme ); + // ptr_list[56] = &(state->Schaake_adjusted_magic_constant_by_soil_type ); ptr_list[57] = &(state->num_lateral_flow_nash_reservoirs); ptr_list[58] = &(state->K_lf); ptr_list[59] = &(state->K_nash); @@ -1672,7 +1726,7 @@ static int Get_state_var_ptrs (Bmi *self, void *ptr_list[]) ptr_list[74] = state->giuh_ordinates; ptr_list[75] = state->nash_storage; ptr_list[76] = state->runoff_queue_m_per_timestep; - ptr_list[77] = state->flux_Schaake_output_runoff_m; + ptr_list[77] = state->flux_output_direct_runoff_m; ptr_list[78] = state->flux_giuh_runoff_m; ptr_list[79] = state->flux_nash_lateral_runoff_m; ptr_list[80] = state->flux_from_deep_gw_to_chan_m; @@ -1680,6 +1734,14 @@ static int Get_state_var_ptrs (Bmi *self, void *ptr_list[]) ptr_list[82] = state->flux_lat_m; ptr_list[83] = state->flux_Qout_m; ptr_list[84] = &(state->verbosity); + //--------------------------------------- + // direct_runoff_params_struct vars + // xinanjiang or schaake flag [56] + //--------------------------------------- + ptr_list[85] = &(state->direct_runoff_params_struct.Schaake_adjusted_magic_constant_by_soil_type ); + ptr_list[86] = &(state->direct_runoff_params_struct.a_Xinanjiang_inflection_point_parameter ); + ptr_list[87] = &(state->direct_runoff_params_struct.b_Xinanjiang_shape_parameter ); + ptr_list[88] = &(state->direct_runoff_params_struct.x_Xinanjiang_shape_parameter ); //------------------------------------------------------------- return BMI_SUCCESS; } @@ -1706,8 +1768,11 @@ static int Get_state_var_ptrs (Bmi *self, void *ptr_list[]) // } // //----------------------------------------------------------------------- -static int Set_state_var (Bmi *self, void *src, int index) +/*static int Set_state_var (Bmi *self, void *src, int index) { + // NOTE: 11.05.2021: this and other serialization functions (*_state_var_*) + // are in current working development outside of master branch + //---------------------------------------------------- // Set the value (or values) for a state variable // using its position index within the state struct. @@ -1830,9 +1895,9 @@ static int Set_state_var (Bmi *self, void *src, int index) // vol_struct vars //---------------------------------------------------------------- else if (index == 35){ - state->vol_struct.vol_sch_runoff = *(double *)src; } + state->vol_struct.vol_runoff = *(double *)src; } else if (index == 36){ - state->vol_struct.vol_sch_infilt = *(double *)src; } + state->vol_struct.vol_infilt = *(double *)src; } else if (index == 37){ state->vol_struct.vol_to_soil = *(double *)src; } else if (index == 38){ @@ -1876,8 +1941,8 @@ static int Set_state_var (Bmi *self, void *src, int index) // forcing_file is a string memcpy(state->forcing_file, src, size); } // state->forcing_file = (char *)src; } // Doesn't work - else if (index == 56){ - state->Schaake_adjusted_magic_constant_by_soil_type = *(double *)src; } + else if (index == 56){ + state->direct_runoff_params_struct.surface_partitioning_scheme = *(int *)src; } else if (index == 57){ state->num_lateral_flow_nash_reservoirs = *(int *)src; } else if (index == 58){ @@ -1943,7 +2008,7 @@ static int Set_state_var (Bmi *self, void *src, int index) state->runoff_queue_m_per_timestep[i] = *( ((double *)src) + i); } } else if (index == 77){ for (i=0; iflux_Schaake_output_runoff_m[i] = *( ((double *)src) + i); } } + state->flux_output_direct_runoff_m[i] = *( ((double *)src) + i); } } else if (index == 78){ for (i=0; iflux_giuh_runoff_m[i] = *( ((double *)src) + i); } } @@ -1965,9 +2030,20 @@ static int Set_state_var (Bmi *self, void *src, int index) else if (index == 84){ // verbosity is not a pointer state->verbosity = *(int *)src; } - + //-------------------------------------------------------------------------- + // direct_runoff_params_struc vars (includes xinanjiang AND schaake) + //-------------------------------------------------------------------------- + else if (index == 85){ + state->direct_runoff_params_struct.Schaake_adjusted_magic_constant_by_soil_type = *(double *)src; } + else if (index == 86){ + state->direct_runoff_params_struct.a_Xinanjiang_inflection_point_parameter = *(double *)src; } + else if (index == 87){ + state->direct_runoff_params_struct.b_Xinanjiang_shape_parameter = *(double *)src; } + else if (index == 88){ + state->direct_runoff_params_struct.x_Xinanjiang_shape_parameter = *(double *)src; } + return BMI_SUCCESS; -} +}*/ /* Grid information */ static int Get_grid_rank (Bmi *self, int grid, int * rank) @@ -2149,7 +2225,12 @@ cfe_state_struct *new_bmi_cfe(void) data->nash_storage = NULL; data->runoff_queue_m_per_timestep = NULL; data->flux_Qout_m = NULL; - data->flux_Schaake_output_runoff_m = NULL; + + /* xinanjiang_dev + changing the name to the more general "direct runoff" + data->flux_Schaake_output_runoff_m = NULL;*/ + data->flux_output_direct_runoff_m = NULL; + data->flux_from_deep_gw_to_chan_m = NULL; data->flux_giuh_runoff_m = NULL; data->flux_lat_m = NULL; @@ -2194,14 +2275,6 @@ Bmi* register_bmi_cfe(Bmi *model) { model->set_value = Set_value; model->set_value_at_indices = Set_value_at_indices; - // New BMI extensions to support serialization - model->get_state_var_count = Get_state_var_count; - model->get_state_var_names = Get_state_var_names; - model->get_state_var_types = Get_state_var_types; - model->get_state_var_ptrs = Get_state_var_ptrs; - model->get_state_var_sizes = Get_state_var_sizes; - model->set_state_var = Set_state_var; - model->get_grid_size = Get_grid_size; model->get_grid_rank = Get_grid_rank; model->get_grid_type = Get_grid_type; @@ -2233,37 +2306,35 @@ extern void run_cfe(cfe_state_struct* cfe_ptr){ cfe_ptr->NWM_soil_params, // Set by config file &cfe_ptr->soil_reservoir, // Set in "init_soil_reservoir" function cfe_ptr->timestep_h, // Set in initialize - cfe_ptr->Schaake_adjusted_magic_constant_by_soil_type, // Set by config file + + /* xinanjiang_dev + changing the name to the more general "direct runoff" + cfe_ptr->Schaake_adjusted_magic_constant_by_soil_type, // Set by config file*/ + cfe_ptr->direct_runoff_params_struct, // Set by config file, includes parameters for Schaake and/or XinanJiang*/ + cfe_ptr->timestep_rainfall_input_m, // Set by bmi (set value) or read from file. - cfe_ptr->flux_Schaake_output_runoff_m, // Set by cfe function + + /* xinanjiang_dev + cfe_ptr->flux_Schaake_output_runoff_m, // Set by cfe function*/ + cfe_ptr->flux_output_direct_runoff_m, + &cfe_ptr->infiltration_depth_m, // Set by Schaake partitioning scheme -// cfe_ptr->flux_overland_m, // Set by CFE function, after Schaake not needed, redundant with flux_Schaake_runoff_m - &cfe_ptr->vol_struct.vol_sch_runoff, // Set by set_volume_trackers_to_zero - &cfe_ptr->vol_struct.vol_sch_infilt, // Set by set_volume_trackers_to_zero cfe_ptr->flux_perc_m, // Set to zero in definition. - &cfe_ptr->vol_struct.vol_to_soil, // Set by set_volume_trackers_to_zero cfe_ptr->flux_lat_m, // Set by CFE function after soil_resevroir calc &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->vol_struct.vol_to_gw, // Set by set_volume_trackers_to_zero - &cfe_ptr->vol_struct.vol_soil_to_gw, // Set by set_volume_trackers_to_zero - &cfe_ptr->vol_struct.vol_soil_to_lat_flow, // Set by set_volume_trackers_to_zero - &cfe_ptr->vol_struct.volout, // Set by set_volume_trackers_to_zero cfe_ptr->flux_from_deep_gw_to_chan_m, // Set by CFE function after gw_reservoir calc - &cfe_ptr->vol_struct.vol_from_gw, // Set by set_volume_trackers_to_zero cfe_ptr->flux_giuh_runoff_m, // Set in CFE by convolution_integral 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 - &cfe_ptr->vol_struct.vol_out_giuh, // Set by set_volume_trackers_to_zero cfe_ptr->flux_nash_lateral_runoff_m, // Set in CFE from nash_cascade function cfe_ptr->num_lateral_flow_nash_reservoirs, // Set from config file cfe_ptr->K_nash, // Set from config file cfe_ptr->nash_storage, // Set from config file - &cfe_ptr->vol_struct.vol_in_nash, // Set by set_volume_trackers_to_zero - &cfe_ptr->vol_struct.vol_out_nash, // Set by set_volume_trackers_to_zero &cfe_ptr->et_struct, // Set to zero with initalize. Set by BMI (set_value) during run - cfe_ptr->flux_Qout_m // Set by CFE function + cfe_ptr->flux_Qout_m, // Set by CFE function + &cfe_ptr->vol_struct ); } @@ -2337,8 +2408,8 @@ extern double init_reservoir_storage(int is_ratio, double amount, double max_amo } extern void initialize_volume_trackers(cfe_state_struct* cfe_ptr){ - cfe_ptr->vol_struct.vol_sch_runoff = 0; - cfe_ptr->vol_struct.vol_sch_infilt = 0; + cfe_ptr->vol_struct.vol_runoff = 0; + cfe_ptr->vol_struct.vol_infilt = 0; cfe_ptr->vol_struct.vol_to_soil = 0; cfe_ptr->vol_struct.vol_to_gw = 0; cfe_ptr->vol_struct.vol_soil_to_gw = 0; @@ -2370,7 +2441,11 @@ extern void print_cfe_flux_at_timestep(cfe_state_struct* cfe_ptr){ printf("%d %lf %lf %lf %lf %lf %lf\n", cfe_ptr->current_time_step, cfe_ptr->timestep_rainfall_input_m*1000.0, - *cfe_ptr->flux_Schaake_output_runoff_m*1000.0, + + /* xinanjiang_dev + *cfe_ptr->flux_Schaake_output_runoff_m*1000.0,*/ + *cfe_ptr->flux_output_direct_runoff_m*1000.0, + *cfe_ptr->flux_giuh_runoff_m*1000.0, *cfe_ptr->flux_nash_lateral_runoff_m*1000.0, *cfe_ptr->flux_from_deep_gw_to_chan_m*1000.0, @@ -2396,7 +2471,11 @@ extern void mass_balance_check(cfe_state_struct* cfe_ptr){ vol_soil_end=cfe_ptr->soil_reservoir.storage_m; double global_residual; - double schaake_residual; + + /* xinanjiang_dev + double schaake_residual;*/ + double direct_residual; + double giuh_residual; double soil_residual; double nash_residual; @@ -2415,28 +2494,47 @@ extern void mass_balance_check(cfe_state_struct* cfe_ptr){ if(!is_fabs_less_than_epsilon(global_residual,1.0e-12)) printf("WARNING: GLOBAL MASS BALANCE CHECK FAILED\n"); + /* xinanjiang_dev schaake_residual = cfe_ptr->vol_struct.volin - cfe_ptr->vol_struct.vol_sch_runoff - cfe_ptr->vol_struct.vol_sch_infilt; printf(" SCHAAKE MASS BALANCE\n"); printf(" surface runoff: %8.4lf m\n",cfe_ptr->vol_struct.vol_sch_runoff); printf(" infiltration: %8.4lf m\n",cfe_ptr->vol_struct.vol_sch_infilt); printf("schaake residual: %6.4e m\n",schaake_residual); // should equal 0.0 if(!is_fabs_less_than_epsilon(schaake_residual,1.0e-12)) - printf("WARNING: SCHAAKE PARTITIONING MASS BALANCE CHECK FAILED\n"); + printf("WARNING: SCHAAKE PARTITIONING MASS BALANCE CHECK FAILED\n");*/ + direct_residual = cfe_ptr->vol_struct.volin - cfe_ptr->vol_struct.vol_runoff - cfe_ptr->vol_struct.vol_infilt; + printf(" DIRECT RUNOFF MASS BALANCE\n"); + printf(" surface runoff: %8.4lf m\n",cfe_ptr->vol_struct.vol_runoff); + printf(" infiltration: %8.4lf m\n",cfe_ptr->vol_struct.vol_infilt); + printf("direct 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"); - giuh_residual = cfe_ptr->vol_struct.vol_sch_runoff - cfe_ptr->vol_struct.vol_out_giuh - vol_end_giuh; + /* xinanjiang_dev + giuh_residual = cfe_ptr->vol_struct.vol_out_giuh - cfe_ptr->vol_struct.vol_sch_runoff - vol_end_giuh; */ + giuh_residual = cfe_ptr->vol_struct.vol_runoff - cfe_ptr->vol_struct.vol_out_giuh - vol_end_giuh; printf(" GIUH MASS BALANCE\n"); - printf(" vol. into giuh: %8.4lf m\n",cfe_ptr->vol_struct.vol_sch_runoff); - printf(" vol. out giuh: %8.4lf m\n",cfe_ptr->vol_struct.vol_out_giuh); - printf(" vol. end giuh q: %8.4lf m\n",vol_end_giuh); + + /* xinanjiang_dev + printf(" vol. into giuh: %8.4lf m\n",cfe_ptr->vol_struct.vol_sch_runoff); */ + printf(" vol. into giuh: %8.4lf m\n",cfe_ptr->vol_struct.vol_runoff); + fprintf(stderr," vol. out giuh: %8.4lf m\n",cfe_ptr->vol_struct.vol_out_giuh); + fprintf(stderr," vol. end giuh q: %8.4lf m\n",cfe_ptr->vol_struct.vol_end_giuh); printf(" giuh 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"); - - soil_residual=cfe_ptr->vol_struct.vol_soil_start + cfe_ptr->vol_struct.vol_sch_infilt - + + /* xinanjiang_dev + soil_residual=cfe_ptr->vol_struct.vol_soil_start + cfe_ptr->vol_struct.vol_sch_infilt - */ + soil_residual=cfe_ptr->vol_struct.vol_soil_start + cfe_ptr->vol_struct.vol_infilt - cfe_ptr->vol_struct.vol_soil_to_lat_flow - vol_soil_end - cfe_ptr->vol_struct.vol_to_gw; + printf(" SOIL WATER CONCEPTUAL RESERVOIR MASS BALANCE\n"); printf(" init soil vol: %8.4lf m\n",cfe_ptr->vol_struct.vol_soil_start); - printf(" vol. into soil: %8.4lf m\n",cfe_ptr->vol_struct.vol_sch_infilt); + + /* xinanjiang_dev + printf(" vol. into soil: %8.4lf m\n",cfe_ptr->vol_struct.vol_sch_infilt); */ + printf(" vol. into soil: %8.4lf m\n",cfe_ptr->vol_struct.vol_infilt); printf("vol.soil2latflow: %8.4lf m\n",cfe_ptr->vol_struct.vol_soil_to_lat_flow); printf(" vol. soil to gw: %8.4lf m\n",cfe_ptr->vol_struct.vol_soil_to_gw); printf(" final vol. soil: %8.4lf m\n",vol_soil_end); diff --git a/src/cfe.c b/src/cfe.c index 69d7217e..010a5ca1 100644 --- a/src/cfe.c +++ b/src/cfe.c @@ -10,37 +10,36 @@ extern void cfe( struct NWM_soil_parameters NWM_soil_params_struct, struct conceptual_reservoir *soil_reservoir_struct, double timestep_h, - double Schaake_adjusted_magic_constant_by_soil_type, + + /* xinanjiang_dev: since we are doing the option for Schaake and XinJiang, + instead of passing in the constants + pass in a structure with the constants for both subroutines. + //double Schaake_adjusted_magic_constant_by_soil_type,*/ + struct direct_runoff_parameters_structure direct_runoff_params_struct, + double timestep_rainfall_input_m, - double *Schaake_output_runoff_m_ptr, + + /* xinanjiang_dev: rename to the general "direct runoff" + double *Schaake_output_runoff_m_ptr,*/ + double *flux_output_direct_runoff_m, + double *infiltration_depth_m_ptr, - // double *flux_overland_m_ptr, FIXME NOT NEEDED, redundant with Schaake_outout_runoff_m_ptr - double *vol_sch_runoff_ptr, - double *vol_sch_infilt_ptr, double *flux_perc_m_ptr, - double *vol_to_soil_ptr, double *flux_lat_m_ptr, double *gw_reservoir_storage_deficit_m_ptr, struct conceptual_reservoir *gw_reservoir_struct, - double *vol_to_gw_ptr, - double *vol_soil_to_gw_ptr, - double *vol_soil_to_lat_flow_ptr, - double *volout_ptr, double *flux_from_deep_gw_to_chan_m_ptr, - double *vol_from_gw_ptr, double *giuh_runoff_m_ptr, int num_giuh_ordinates, double *giuh_ordinates_arr, double *runoff_queue_m_per_timestep_arr, - double *vol_out_giuh_ptr, double *nash_lateral_runoff_m_ptr, int num_lateral_flow_nash_reservoirs, double K_nash, double *nash_storage_arr, - double *vol_in_nash_ptr, - double *vol_out_nash_ptr, struct evapotranspiration_structure *evap_struct, - double *Qout_m_ptr + double *Qout_m_ptr, + struct massbal *massbal_struct ){ // ####################################################################### // CFE STATE SPACE FUNCTION // ####################################################################### @@ -48,9 +47,12 @@ extern void cfe( // #### Reason: so we don't have to re-write domain science code to de-reference a whole bunch of pointers // #### Note: all of thes variables are storages in [m] or fluxes in [m/timestep] double soil_reservoir_storage_deficit_m = *soil_reservoir_storage_deficit_m_ptr; // storage [m] - double Schaake_output_runoff_m = *Schaake_output_runoff_m_ptr; // Schaake partitioned runoff this timestep [m] + + /* xinanjiang_dev: rename to the general "direct runoff" + double Schaake_output_runoff_m = *Schaake_output_runoff_m_ptr; // Schaake partitioned runoff this timestep [m]*/ + double direct_output_runoff_m = *flux_output_direct_runoff_m; // Schaake partitioned runoff this timestep [m]*/ + double infiltration_depth_m = *infiltration_depth_m_ptr; // Schaake partitioned infiltration this timestep [m] - // double flux_overland_m = *flux_overland_m_ptr; FIXME NOT NEEDED, redundant with Schaake_output_runoff_m double flux_perc_m = *flux_perc_m_ptr; // water moved from soil reservoir to gw reservoir this timestep [m] 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] @@ -58,17 +60,6 @@ extern void cfe( double giuh_runoff_m = *giuh_runoff_m_ptr; // water leaving GIUH 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 (GIUH+Nash+GW) [m] - double vol_sch_runoff = *vol_sch_runoff_ptr; // mass balance change in storage [m] this timestep or flux [m/timestep] - double vol_sch_infilt = *vol_sch_infilt_ptr; // mass balance change in storage [m] this timestep or flux [m/timestep] - double vol_to_soil = *vol_to_soil_ptr; // mass balance change in storage [m] this timestep or flux [m/timestep] - double vol_to_gw = *vol_to_gw_ptr; // mass balance change in storage [m] this timestep or flux [m/timestep] - double vol_soil_to_gw = *vol_soil_to_gw_ptr; // mass balance change in storage [m] this timestep or flux [m/timestep] - double vol_soil_to_lat_flow = *vol_soil_to_lat_flow_ptr; // mass balance change in storage [m] this timestep or flux [m/timestep] - double vol_from_gw = *vol_from_gw_ptr; // mass balance change in storage [m] this timestep or flux [m/timestep] - double vol_out_giuh = *vol_out_giuh_ptr; // mass balance change in storage [m] this timestep or flux [m/timestep] - double vol_in_nash = *vol_in_nash_ptr; // mass balance change in storage [m] this timestep or flux [m/timestep] - double vol_out_nash = *vol_out_nash_ptr; // mass balance change in storage [m] this timestep or flux [m/timestep] - double volout = *volout_ptr; // mass balance change in storage [m] this timestep or flux [m/timestep] // LOCAL VARIABLES, the values of which are not important to describe the model state. They are like notes on scrap paper. @@ -83,50 +74,70 @@ extern void cfe( //################################################## soil_reservoir_storage_deficit_m=(NWM_soil_params_struct.smcmax*NWM_soil_params_struct.D-soil_reservoir_struct->storage_m); - - if(timestep_rainfall_input_m >0.0) // Call Schaake function only if rainfall input this time step is nonzero. + + // NEW FLO + if(0.0 < timestep_rainfall_input_m) { - Schaake_partitioning_scheme(timestep_h,Schaake_adjusted_magic_constant_by_soil_type,soil_reservoir_storage_deficit_m, - timestep_rainfall_input_m,&Schaake_output_runoff_m,&infiltration_depth_m); + if (direct_runoff_params_struct.surface_partitioning_scheme == Schaake) + { + Schaake_partitioning_scheme(timestep_h,direct_runoff_params_struct.Schaake_adjusted_magic_constant_by_soil_type,soil_reservoir_storage_deficit_m, + timestep_rainfall_input_m,&direct_output_runoff_m,&infiltration_depth_m); + } + else if (direct_runoff_params_struct.surface_partitioning_scheme == Xinanjiang) + { + 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, + &direct_runoff_params_struct, + &direct_output_runoff_m, &infiltration_depth_m); + } + else + { + fprintf(stderr,"Problem, must specify one of Schaake of Xinanjiang partitioning scheme.\n"); + fprintf(stderr,"Program terminating.\n"); + exit(-1); // note -1 is arbitrary #############BOMB################ NEW FLO + } } - else // No need to call the Schaake function. + else // NEW FLO no need to call Schaake or Xinanjiang if it's not raining. { - Schaake_output_runoff_m=0.0; - infiltration_depth_m=0.0; + direct_output_runoff_m = 0.0; + infiltration_depth_m = 0.0; } - + // check to make sure that there is storage available in soil to hold the water that does not runoff //-------------------------------------------------------------------------------------------------- if(soil_reservoir_storage_deficit_mstorage_m=soil_reservoir_struct->storage_max_m; } #ifdef DEBUG + /* xinanjiang_dev printf("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); -#endif - - vol_sch_runoff += Schaake_output_runoff_m; - vol_sch_infilt += infiltration_depth_m; + 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", + timestep_rainfall_input_m,direct_output_runoff_m*1000.0,infiltration_depth_m*1000.0, + timestep_rainfall_input_m-direct_output_runoff_m-infiltration_depth_m); +#endif + + massbal_struct->vol_runoff += direct_output_runoff_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 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - // REDUNDANT soil_reservoir_storage_deficit_m=soil_reservoir.storage_max_m-soil_reservoir.storage_m; <- commented out by FLO based on comments from Bartel 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; - vol_sch_runoff+=diff; // send excess water back to GIUH runoff - vol_sch_infilt-=diff; // correct overprediction of infilt. - Schaake_output_runoff_m+=diff; // bug found by Nels. This was missing and fixes it. + 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 + direct_output_runoff_m+=diff; // bug found by Nels. This was missing and fixes it. } - 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. @@ -148,23 +159,23 @@ extern void cfe( { diff=flux_perc_m-gw_reservoir_storage_deficit_m; flux_perc_m=gw_reservoir_storage_deficit_m; - vol_sch_runoff+=diff; // send excess water back to GIUH runoff - vol_sch_infilt-=diff; // correct overprediction of infilt. + massbal_struct->vol_runoff+=diff; // send excess water back to GIUH runoff + massbal_struct->vol_infilt-=diff; // correct overprediction of infilt. } - vol_to_gw +=flux_perc_m; - vol_soil_to_gw +=flux_perc_m; + 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; - vol_soil_to_lat_flow += flux_lat_m; //TODO add this to nash cascade as input - volout=volout+flux_lat_m; + massbal_struct->vol_soil_to_lat_flow += flux_lat_m; //TODO add this to nash cascade as input + massbal_struct->volout=massbal_struct->volout+flux_lat_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 >>>>>>>>> - vol_from_gw+=flux_from_deep_gw_to_chan_m; + 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"); @@ -178,18 +189,20 @@ extern void cfe( // Solve the convolution integral ffor this time step - giuh_runoff_m = convolution_integral(Schaake_output_runoff_m,num_giuh_ordinates, + /* xinanjiang_dev + giuh_runoff_m = convolution_integral(Schaake_output_runoff_m,num_giuh_ordinates, */ + giuh_runoff_m = convolution_integral(direct_output_runoff_m,num_giuh_ordinates, giuh_ordinates_arr,runoff_queue_m_per_timestep_arr); - vol_out_giuh+=giuh_runoff_m; + massbal_struct->vol_out_giuh+=giuh_runoff_m; - volout+=giuh_runoff_m; - volout+=flux_from_deep_gw_to_chan_m; + massbal_struct->volout+=giuh_runoff_m; + massbal_struct->volout+=flux_from_deep_gw_to_chan_m; // Route lateral flow through the Nash cascade. nash_lateral_runoff_m = nash_cascade(flux_lat_m,num_lateral_flow_nash_reservoirs, K_nash,nash_storage_arr); - vol_in_nash += flux_lat_m; - vol_out_nash += nash_lateral_runoff_m; + massbal_struct->vol_in_nash += flux_lat_m; + massbal_struct->vol_out_nash += nash_lateral_runoff_m; #ifdef DEBUG fprintf(out_debug_fptr,"%d %lf %lf\n",tstep,flux_lat_m,nash_lateral_runoff_m); @@ -199,26 +212,18 @@ extern void cfe( // #### COPY BACK STATE VALUES BY POINTER REFERENCE SO VISIBLE TO FRAMEWORK #### *soil_reservoir_storage_deficit_m_ptr = soil_reservoir_storage_deficit_m; - *Schaake_output_runoff_m_ptr = Schaake_output_runoff_m; + + /* xinanjiang_dev + *Schaake_output_runoff_m_ptr = Schaake_output_runoff_m; */ + *flux_output_direct_runoff_m = direct_output_runoff_m; + *infiltration_depth_m_ptr = infiltration_depth_m; -// *flux_overland_m_ptr = Schaake_output_runoff_m; NOT NEEDED, They are the same thing. - *vol_sch_runoff_ptr = vol_sch_runoff; - *vol_sch_infilt_ptr = vol_sch_infilt; *flux_perc_m_ptr = flux_perc_m; - *vol_to_soil_ptr = vol_to_soil; *flux_lat_m_ptr = flux_lat_m; *gw_reservoir_storage_deficit_m_ptr = gw_reservoir_storage_deficit_m; - *vol_to_gw_ptr = vol_to_gw; - *vol_soil_to_gw_ptr = vol_soil_to_gw; - *vol_soil_to_lat_flow_ptr = vol_soil_to_lat_flow; - *volout_ptr = volout; *flux_from_deep_gw_to_chan_m_ptr = flux_from_deep_gw_to_chan_m; - *vol_from_gw_ptr = vol_from_gw; *giuh_runoff_m_ptr = giuh_runoff_m; - *vol_out_giuh_ptr = vol_out_giuh; *nash_lateral_runoff_m_ptr = nash_lateral_runoff_m; - *vol_in_nash_ptr = vol_in_nash; - *vol_out_nash_ptr = vol_out_nash; *Qout_m_ptr = Qout_m; @@ -444,6 +449,215 @@ else return; } + +//############################################################## +//######## XINANJIANG RUNOFF PARTITIONING SCHEME ########### +//############################################################## + +void Xinanjiang_partitioning_scheme(double water_input_depth_m, double field_capacity_m, + double max_soil_moisture_storage_m, double column_total_soil_water_m, + struct direct_runoff_parameters_structure *parms, + double *surface_runoff_depth_m, double *infiltration_depth_m) +{ + //------------------------------------------------------------------------ + // This module takes the water_input_depth_m and separates it into surface_runoff_depth_m + // and infiltration_depth_m by calculating the saturated area and runoff based on a scheme developed + // for the Xinanjiang model by Jaywardena and Zhou (2000). According to Knoben et al. + // (2019) "the model uses a variable contributing area to simulate runoff. [It] uses + // a double parabolic curve to simulate tension water capacities within the catchment, + // instead of the original single parabolic curve" which is also used as the standard + // VIC fomulation. This runoff scheme was selected for implementation into NWM v3.0. + // REFERENCES: + // 1. Jaywardena, A.W. and M.C. Zhou, 2000. A modified spatial soil moisture storage + // capacity distribution curve for the Xinanjiang model. Journal of Hydrology 227: 93-113 + // 2. Knoben, W.J.M. et al., 2019. Supplement of Modular Assessment of Rainfall-Runoff Models + // Toolbox (MARRMoT) v1.2: an open-source, extendable framework providing implementations + // of 46 conceptual hydrologic models as continuous state-space formulations. Supplement of + // Geosci. Model Dev. 12: 2463-2480. + //------------------------------------------------------------------------- + // Written by RLM May 2021 + // Adapted by JMFrame September 2021 for new version of CFE + //------------------------------------------------------------------------- + // Inputs + // double water_input_depth_m amount of water input to soil surface this time step [m] + // double field_capacity_m amount of water stored in soil reservoir when at field capacity [m] + // double max_soil_moisture_storage_m total storage of the soil moisture reservoir (porosity*soil thickness) [m] + // double column_total_soil_water_m current storage of the soil moisture reservoir [m] + // double a_inflection_point_parameter a parameter + // double b_shape_parameter b parameter + // double x_shape_parameter x parameter + // + // Outputs + // double surface_runoff_depth_m amount of water partitioned to surface water this time step [m] + // double infiltration_depth_m amount of water partitioned as infiltration (soil water input) this time step [m] + //------------------------------------------------------------------------- + + double tension_water_m, free_water_m, max_tension_water_m, max_free_water_m, pervious_runoff_m; + + //could move this if statement outside of both the Schaake and Xinanjiang subroutines edit FLO- moved to main(). + + // partition the total soil water in the column between free water and tension water + free_water_m = column_total_soil_water_m - field_capacity_m; + + if(0.0 < free_water_m) { //edit FLO + tension_water_m = field_capacity_m; + } else { + free_water_m = 0.0; + tension_water_m = column_total_soil_water_m; + } + + // estimate the maximum free water and tension water available in the soil column + max_free_water_m = max_soil_moisture_storage_m - field_capacity_m; + max_tension_water_m = field_capacity_m; + + // check that the free_water_m and tension_water_m do not exceed the maximum and if so, change to the max value + if(max_free_water_m < free_water_m) free_water_m = max_free_water_m; + if(max_tension_water_m < tension_water_m) tension_water_m = max_tension_water_m; + + // NOTE: the impervious surface runoff assumptions due to frozen soil used in NWM 3.0 have not been included. + // We are assuming an impervious area due to frozen soils equal to 0 (see eq. 309 from Knoben et al). + + // The total (pervious) runoff is first estimated before partitioning into surface and subsurface components. + // See Knoben et al eq 310 for total runoff and eqs 313-315 for partitioning between surface and subsurface + // components. + + // Calculate total estimated pervious runoff. + // NOTE: If the impervious surface runoff due to frozen soils is added, + // the pervious_runoff_m equation will need to be adjusted by the fraction of pervious area. + if ((tension_water_m/max_tension_water_m) <= (0.5 - parms->a_Xinanjiang_inflection_point_parameter)) { + pervious_runoff_m = water_input_depth_m * (pow((0.5 - parms->a_Xinanjiang_inflection_point_parameter), + (1.0 - parms->b_Xinanjiang_shape_parameter)) * + pow((1.0 - (tension_water_m/max_tension_water_m)), + parms->b_Xinanjiang_shape_parameter)); + + } else { + pervious_runoff_m = water_input_depth_m * (1.0 - pow((0.5 + parms->a_Xinanjiang_inflection_point_parameter), + (1.0 - parms->b_Xinanjiang_shape_parameter)) * + pow((1.0 - (tension_water_m/max_tension_water_m)), + (parms->b_Xinanjiang_shape_parameter))); + } + // Separate the surface water from the pervious runoff + // NOTE: If impervious runoff is added to this subroutine, impervious runoff should be added to + // the surface_runoff_depth_m. + *surface_runoff_depth_m = pervious_runoff_m * (1.0 - pow((1.0 - (free_water_m/max_free_water_m)),parms->x_Xinanjiang_shape_parameter)); + + // The surface runoff depth is bounded by a minimum of 0 and a maximum of the water input depth. + // Check that the estimated surface runoff is not less than 0.0 and if so, change the value to 0.0. + if(*surface_runoff_depth_m < 0.0) *surface_runoff_depth_m = 0.0; + // Check that the estimated surface runoff does not exceed the amount of water input to the soil surface. If it does, + // change the surface water runoff value to the water input depth. + if(*surface_runoff_depth_m > water_input_depth_m) *surface_runoff_depth_m = water_input_depth_m; + // Separate the infiltration from the total water input depth to the soil surface. + *infiltration_depth_m = water_input_depth_m - *surface_runoff_depth_m; + +return; +} +//############################################################## +//######## XINANJIANG RUNOFF PARTITIONING SCHEME ########### +/*############################################################## + +void Xinanjiang_partitioning_scheme(double water_input_depth_m, double field_capacity_m, + double max_soil_moisture_storage_m, double column_total_soil_water_m, + struct direct_runoff_parameters_structure *parms, + double *surface_runoff_depth_m, double *infiltration_depth_m) +{ + //------------------------------------------------------------------------ + // This module takes the water_input_depth_m and separates it into surface_runoff_depth_m + // and infiltration_depth_m by calculating the saturated area and runoff based on a scheme developed + // for the Xinanjiang model by Jaywardena and Zhou (2000). According to Knoben et al. + // (2019) "the model uses a variable contributing area to simulate runoff. [It] uses + // a double parabolic curve to simulate tension water capacities within the catchment, + // instead of the original single parabolic curve" which is also used as the standard + // VIC fomulation. This runoff scheme was selected for implementation into NWM v3.0. + // REFERENCES: + // 1. Jaywardena, A.W. and M.C. Zhou, 2000. A modified spatial soil moisture storage + // capacity distribution curve for the Xinanjiang model. Journal of Hydrology 227: 93-113 + // 2. Knoben, W.J.M. et al., 2019. Supplement of Modular Assessment of Rainfall-Runoff Models + // Toolbox (MARRMoT) v1.2: an open-source, extendable framework providing implementations + // of 46 conceptual hydrologic models as continuous state-space formulations. Supplement of + // Geosci. Model Dev. 12: 2463-2480. + //------------------------------------------------------------------------- + // Written by RLM May 2021 + // Adapted by JMFrame September 2021 for new version of CFE + //------------------------------------------------------------------------- + // Inputs + // double water_input_depth_m amount of water input to soil surface this time step [m] + // double field_capacity_m + // double max_soil_moisture_storage_m + // double column_total_soil_water_m + // double a_inflection_point_parameter + // double b_shape_parameter + // double x_shape_parameter + // + // Outputs + // double surface_runoff_depth_m amount of water partitioned to surface water this time step [m] + // double infiltration_depth_m amount of water partitioned as infiltration (soil water input) this time step [m] + //------------------------------------------------------------------------- + + double tension_water_m, free_water_m, max_tension_water_m, max_free_water_m, pervious_runoff_m; + + if(0.0 < water_input_depth_m) { //could move this if statement outside of both the schaake and xinanjiang subroutines + + // partition the total soil water in the column between free water and tension water + free_water_m = column_total_soil_water_m - field_capacity_m; + + if(free_water_m > 0) { + tension_water_m = field_capacity_m; + } else { + free_water_m = 0.0; + tension_water_m = column_total_soil_water_m; + } + + // estimate the maximum free water and tension water available in the soil column + max_free_water_m = max_soil_moisture_storage_m - field_capacity_m; + max_tension_water_m = field_capacity_m; + + // check that the free_water_m and tension_water_m do not exceed the maximum and if so, change to the max value + if(max_free_water_m < free_water_m) free_water_m = max_free_water_m; + if(max_tension_water_m < tension_water_m) tension_water_m = max_tension_water_m; + + // NOTE: the impervious surface runoff assumptions due to frozen soil used in NWM 3.0 have not been included. + // We are assuming an impervious area due to frozen soils equal to 0 (see eq. 309 from Knoben et al). + + // The total (pervious) runoff is first estimated before partitioning into surface and subsurface components. + // See Knoben et al eq 310 for total runoff and eqs 313-315 for partitioning between surface and subsurface + // components. + + // Calculate total estimated pervious runoff. + // NOTE: If the impervious surface runoff due to frozen soils is added, + // the pervious_runoff_m equation will need to be adjusted by the fraction of pervious area. + if ((tension_water_m/max_tension_water_m) <= (0.5 - parms->a_Xinanjiang_inflection_point_parameter)) { + pervious_runoff_m = water_input_depth_m * (pow((0.5 - parms->a_Xinanjiang_inflection_point_parameter), + (1.0 - parms->b_Xinanjiang_shape_parameter)) * + pow((1.0 - (tension_water_m/max_tension_water_m)), + parms->b_Xinanjiang_shape_parameter)); + + } else { + pervious_runoff_m = water_input_depth_m * (1.0 - pow((0.5 + parms->a_Xinanjiang_inflection_point_parameter), + (1.0 - parms->b_Xinanjiang_shape_parameter)) * + pow((1.0 - (tension_water_m/max_tension_water_m)), + (parms->b_Xinanjiang_shape_parameter))); + } + // Separate the surface water from the pervious runoff + // NOTE: If impervious runoff is added to this subroutine, impervious runoff should be added to + // the surface_runoff_depth_m. + *surface_runoff_depth_m = pervious_runoff_m * (1.0 - pow((1.0 - (free_water_m/max_free_water_m)),parms->x_Xinanjiang_shape_parameter)); + // The surface runoff depth is bounded by a minimum of 0 and a maximum of the water input depth. + // Check that the estimated surface runoff is not less than 0.0 and if so, change the value to 0.0. + if(*surface_runoff_depth_m < 0.0) *surface_runoff_depth_m = 0.0; + // Check that the estimated surface runoff does not exceed the amount of water input to the soil surface. If it does, + // change the surface water runoff value to the water input depth. + if(*surface_runoff_depth_m > water_input_depth_m) *surface_runoff_depth_m = water_input_depth_m; + // Separate the infiltration from the total water input depth to the soil surface. + *infiltration_depth_m = water_input_depth_m - *surface_runoff_depth_m; + + } else { + *surface_runoff_depth_m = 0.0; + *infiltration_depth_m = 0.0; + } + return; +}*/ + //############################################################## //#################### ET FROM RAINFALL #################### //############################################################## diff --git a/src/main.c b/src/main.c index f895c11e..8335cc9b 100644 --- a/src/main.c +++ b/src/main.c @@ -12,6 +12,16 @@ This is not part of BMI, but acts as the driver that calls the model. int main(int argc, const char *argv[]) { + + //////////////////////////////////////////////////////////////// + ////////////// USING UPDATE ////////////////////////////// + //////////////////////////////////////////////////////////////// + ////////////// USING UPDATE ////////////////////////////// + //////////////////////////////////////////////////////////////// + ////////////// USING UPDATE ////////////////////////////// + //////////////////////////////////////////////////////////////// + printf("Running CFE using Update\n"); + printf("_____________________________\n"); // A configuration file is required for running this model through BMI if(argc<=1){ @@ -39,7 +49,7 @@ int if (cfe_main_data->verbosity > 0) print_cfe_flux_header(); - for (int i = 0; i < 720; i++){ + for (int i = 0; i < 45; i++){ cfe_bmi_model->update(cfe_bmi_model); @@ -54,5 +64,58 @@ int printf("Finalizing CFE model\n"); cfe_bmi_model->finalize(cfe_bmi_model); + //////////////////////////////////////////////////////////////// + ////////////// USING UPDATE UNTIL //////////////////////// + //////////////////////////////////////////////////////////////// + ////////////// USING UPDATE UNTIL //////////////////////// + //////////////////////////////////////////////////////////////// + ////////////// USING UPDATE UNTIL //////////////////////// + //////////////////////////////////////////////////////////////// + printf("Running CFE using UPDATE UNTIL\n"); + printf("_____________________________\n"); + + printf("allocating memory to store entire BMI structure for CFE\n"); + // allocating memory to store the entire BMI structure for CFE + Bmi *cfe_bmi_model2 = (Bmi *) malloc(sizeof(Bmi)); + + printf("Registering BMI CFE model\n"); + register_bmi_cfe(cfe_bmi_model2); + + printf("Initializeing BMI CFE model\n"); + const char *cfg_file2 = argv[1]; + cfe_bmi_model2->initialize(cfe_bmi_model2, cfg_file2); + + printf("Get the information from the configuration here in Main\n"); + // Get the information from the configuration here in Main + cfe_state_struct *cfe_main_data2; + cfe_main_data2 = (cfe_state_struct *) cfe_bmi_model2->data; + + /************************************************************************* + This will be used to advance the model with update_until + **************************************************************************/ + double model_time_step_size; + cfe_bmi_model2->get_time_step(cfe_bmi_model2, &model_time_step_size); + printf("The model time step size is: %lf\n", model_time_step_size); + + printf("Calling update_until()\n"); + if (cfe_main_data2->verbosity > 0) + print_cfe_flux_header(); + + // Run the model with the update until function + cfe_bmi_model2->update_until(cfe_bmi_model2, 45 * model_time_step_size); + + if (cfe_main_data->verbosity > 0) + print_cfe_flux_at_timestep(cfe_main_data); + + // Run the Mass Balance check + //mass_balance_check(cfe_main_data); + + printf("Finalizing CFE model\n"); + cfe_bmi_model2->finalize(cfe_bmi_model2); + + + ///////////////////////////////////////////////////////// + ///////////////////////////////////////////////////////// + return 0; } diff --git a/src/main_cfe_aorc_pet.c b/src/main_cfe_aorc_pet.c index 19134e21..daee76cc 100644 --- a/src/main_cfe_aorc_pet.c +++ b/src/main_cfe_aorc_pet.c @@ -159,6 +159,13 @@ int aorc1 = (aorc_model *) aorc_bmi_model->data; printf("forcing file for the AORC module %s\n", aorc1->forcing_file); + /************************************************************************* + This will be used to advance the model with update_until + **************************************************************************/ + double model_time_step_size; + cfe_bmi_model->get_time_step(cfe_bmi_model, &model_time_step_size); + printf("The model time step size is: %lf\n", model_time_step_size); + /************************************************************************ This is the basic process for getting the three things to talk through BMI 1. Update the AORC forcing data @@ -211,7 +218,7 @@ int printf("PET value from CFE is %8.9lf\n", cfe1->et_struct.potential_et_m_per_s); } - cfe_bmi_model->update(cfe_bmi_model); //Update model 2 + cfe_bmi_model->update_until(cfe_bmi_model, (i+1)*model_time_step_size); // Update model 2 if (cfe1->verbosity > 0) print_cfe_flux_at_timestep(cfe1); diff --git a/src/main_pass_forcings.c b/src/main_pass_forcings.c index 5ad83c61..31bf7a3c 100644 --- a/src/main_pass_forcings.c +++ b/src/main_pass_forcings.c @@ -84,21 +84,28 @@ int aorc_model *aorc; aorc = (aorc_model *) aorc_bmi_model->data; + /************************************************************************* + This will be used to advance the model with update_until + **************************************************************************/ + double model_time_step_size; + cfe_bmi_model->get_time_step(cfe_bmi_model, &model_time_step_size); + printf("The model time step size is: %lf\n", model_time_step_size); + /************************************************************************ This is the basic process for getting two things to talk through BMI 1. Update the AORC forcing data 2. Getting forcing from AORC and setting forcing for CFE 3. Update the CFE model. ************************************************************************/ - + /************************************************************************ Now loop through time and call the models with the intermediate get/set ************************************************************************/ printf("looping through and calling updata\n"); if (cfe_model_data->verbosity > 0) print_cfe_flux_header(); - for (int i = 0; i < 720; i++){ - + for (int i = 0; i < 45; i++){ + aorc_bmi_model->update(aorc_bmi_model); // Update model 1 pass_forcing_from_aorc_to_cfe(cfe_bmi_model, aorc_bmi_model); // Get and Set values @@ -109,7 +116,7 @@ int printf("precip value from CFE is %lf\n", cfe_model_data->aorc.precip_kg_per_m2); } - cfe_bmi_model->update(cfe_bmi_model); // Update model 2 + cfe_bmi_model->update_until(cfe_bmi_model, (i+1)*model_time_step_size); // Update model 2 if (cfe_model_data->verbosity > 0) print_cfe_flux_at_timestep(cfe_model_data); diff --git a/test/main_unit_test_bmi.c b/test/main_unit_test_bmi.c index b061b2f0..66df6cc2 100644 --- a/test/main_unit_test_bmi.c +++ b/test/main_unit_test_bmi.c @@ -386,11 +386,12 @@ main(int argc, const char *argv[]){ free(names_out); free(names_in); // Test BMI: CONTROL FUNCTION update_until() + { int added_nstep=5; int total_nstep= added_nstep + test_nstep; printf("\n updating until... new total timesteps in test loop: %i\n", total_nstep); - status = model->update_until(model,total_nstep); + status = model->update_until(model,total_nstep*dt); if (status == BMI_FAILURE) return BMI_FAILURE; // confirm updated current time model->get_current_time(model, &now); diff --git a/test_serialize/README.md b/test_serialize/README.md deleted file mode 100644 index 2ad04446..00000000 --- a/test_serialize/README.md +++ /dev/null @@ -1,67 +0,0 @@ -# CFE Serialization and Deserialization Test - -Several new functions have been added to the Basic Model Interface (BMI) -in order to support the serialization and deserialization of all of a -model's state variables. State variables are variables that describe -the current state of the model at any given time and that persist between -BMI function calls. This excludes variables that are created only for -temporary use within any given function. A model's state variables will -typically not be defined until after the model's BMI initialize() function -has been called. When a model's BMI update() function is called, the values -of many of the state variables will be modified (in-place) as the model -advances in time. If, after one or more BMI update() calls, all of a -model's state variables are serialized and saved (i.e. to a buffer or a -file), it becomes -possible to deserialize them and to then set them as the state of a new -instance of the same model. The new instance can then continue running -from the point in time that the original model's variables were serialized. - -Serialization is a complex and error-prone process that model developers -should not have to be concerned with. The strategy employed here therefore -follows a "separation of concerns" philosophy where the model developer -only implements new BMI functions that allow the model coupling -framework (e.g. nextGen) to get and set the values of all state variables -and also to retrieve metadata about them. A separate, general-purpose -framework utility can then be written to perform the tasks of serialization -and deserialization that calls and relies on these new BMI functions. -The current implementation of this general utility uses the -[msgpack-c library](https://github.com/msgpack/msgpack-c) -and can be found in -```cfe/test_serialize/serialize_state.c.``` -That file implements three functions: -```serialize(), deserialize_to_state(), and compare_states().``` -Msgpack serializes to binary and is both fast and compact. - -**Note:** You must install the msgpack-c library to run the test. -See instructions at: -[msgpack-c README](https://github.com/msgpack/msgpack-c/blob/c_master/README.md). - -The new BMI functions that support serialization of models written in C -are as follows. These new functions are defined in -```cfe/include/bmi.h```. -Their implementations are straight-forward but model specific. -They have been implemented for the CFE model in: -```cfe/src/bmi_ccfe.c``` - -``` -get_state_var_count(): Return the total number of state variables -get_state_var_names(): Return string array of state variable internal names -get_state_var_types(): Return string array of state variable C types -get_state_var_ptrs(): Return pointer array to all state variables -get_state_var_sizes(): Return unsigned int array of sizes (number of elements) -set_state_var(): Set the value of a state variable given its index -``` - -The code for the actual test can be found in: -```cfe/test_serialize/cfe_serialize_test.c```, -along with a small shell script called -```make_and_run_ser_test.sh``` for making and running the test. - -**Note:** In ```serialize_state.c```, there are constants called -```BUFFER_SIZE``` and ```UNPACKED_BUFFER_SIZE``` that must be big -enough to accommodate all state variables. These are set dynamically -using the size of the file to which the serialized state was written. - -**Note:** To print and inspect the values of all state variables after -they have been deserialized (usually from a file), you should change the -parameter ```print_obj``` from 0 to 1 in ```cfe_serialize_test.c```. diff --git a/test_serialize/cfe_serialize_test.c b/test_serialize/cfe_serialize_test.c deleted file mode 100644 index 223bd840..00000000 --- a/test_serialize/cfe_serialize_test.c +++ /dev/null @@ -1,234 +0,0 @@ -#include -#include -#include // for access() -#include "../include/bmi.h" -#include "../include/bmi_cfe.h" -#include "../include/serialize_state.h" - -/* -This main program creates two instances of a BMI-enabled version of -the CFE model. It's purpose is to test a new set of BMI functions that -are intended to support model serialization for models written in C. - -For Model 1, bmi.initialize() is called, and bmi.update() is called -several times to advance the model state. Then serialize_state() is -called to retrieve the entire model state. - -For Model 2, bmi.initialize() may be called, and then the function -deserialize_to_state() is called to put it in the same state as Model 1. - -For both Models, bmi.update() is then called several more times and -results are compared with the compare_states() function. - -In C, we cannot infer the type or size of an array, so we need BMI -functions to get: name, type, size and ptr. -Recall that array name is a pointer to first element. - -Arrays are stored "flattened", so the model will take care of reading -values into 1D, 2D, 3D arrays. We should only need the total number -of array elements. - -Later on, ser_file should be an argument to main(), not hard-coded. -*/ - -//------------------------------------------------------------------------ -int print_some(void *ptr_list[]){ - - //------------------------------------------------ - // NOTE! Typecast ptr first, then add offset, - // into array, then dereference the ptr. - // Careful with order of operations. - //------------------------------------------------ - puts("Printing some selected variables ..."); - - printf("ptr_list[42] = volout = %f\n", *(double *)ptr_list[42]); - printf("ptr_list[43] = volin = %f\n", *(double *)ptr_list[43]); - printf("ptr_list[51] = num_time_steps = %d\n", *(int *)ptr_list[51]); - printf("ptr_list[52] = current_time_step = %d\n", *(int *)ptr_list[52]); - printf("ptr_list[53] = time_step_size = %d\n", *(int *)ptr_list[53]); - //---------------------------------------------- - // The following 3 methods all work, but don't - // use "&" in get_state_var_ptrs(). - //---------------------------------------------- - printf("ptr_list[55] = forcing_file = %s\n", (char *)ptr_list[55]); - //--------------------------------------------------------------------- - // printf("ptr_list[55] = forcing_file = %s", ptr_list[55]); - //--------------------------------------------------------------------- -// char *p; -// // p =(char *)ptr_list[55]; -// p = ptr_list[55]; -// printf("ptr_list[55] = forcing_file = "); -// while (*p != '\0'){ -// printf("%c", *p++); } - //--------------------------------------------------------------- - printf("ptr_list[71] = aorc.time = %ld\n", *(long *)ptr_list[71]); - //--------------------------------------------------------------- - puts(""); // newline is added - - return 0; -} - -//------------------------------------------------------------------------ -// int main(void) -int main(int argc, const char *argv[]) -{ - // Check for config file arg - const char *cfg_file; - if(argc<=1){ - printf("WARNING: Missing config file argument.\n"); - printf(" Using default.\n\n"); - cfg_file = "./configs/cat_89_bmi_config_cfe.txt"; - } else { - cfg_file = argv[1]; - } - if( access( cfg_file, F_OK ) != 0 ) { - printf("ERROR: cfg_file not found.\n"); - puts( cfg_file); - puts(""); - exit(1); - } - //-------------------------------------------------------------- - const char *ser_file = "./test_serialize/model_state.ser"; // make arg later - int n_steps1 = 10; // n_steps for Model1 before serializing - int n_steps2 = 50; // n_steps for models after deserializing - int verbose = 1; - int print_obj = 1; // Set to 1 to print values after deserializing - int n_state_vars; - - //-------------------------------------------------------------- - if (verbose){ - puts(""); - puts("Allocating memory for BMI CFE instances 1 & 2 ..."); - } - Bmi *model1 = (Bmi *) malloc(sizeof(Bmi)); - Bmi *model2 = (Bmi *) malloc(sizeof(Bmi)); - - //-------------------------------------------------------------- - if (verbose){ puts("Registering CFE models 1 & 2 ..."); } - - register_bmi_cfe(model1); - register_bmi_cfe(model2); - - //-------------------------------------------------------------- - if (verbose){ puts("Initializing CFE models 1 & 2 ..."); } - - model1->initialize(model1, cfg_file); - model2->initialize(model2, cfg_file); - - //-------------------------------------------------------------- - if (verbose){ - puts("Updating CFE model 1 ..."); - printf(" n_steps1 = %i \n", n_steps1); - puts(""); - } - - for (int i=1; i<=n_steps1; i++){ - model1->update(model1); - } - - //-------------------------------------------------------------- - if (verbose){ - puts("Calling BMI.get_state_var_ptrs() on CFE model 1 ..."); - } - - model1->get_state_var_count(model1, &n_state_vars); - - //--------------------------------------------- - // For testing: All 3 print "8" on my MacPro - //--------------------------------------------- - //printf("Size of void* = %lu\n", sizeof(void*)); - //printf("Size of int* = %lu\n", sizeof(int*)); - //printf("Size of double* = %lu\n", sizeof(double*)); - //printf("\n"); - - //--------------------------------------------------------------- - // See: https://stackoverflow.com/questions/7798383/ - // array-of-pointers-to-multiple-types-c/7799543 - //--------------------------------------------------------------- - void *ptr_list[ n_state_vars ]; - model1->get_state_var_ptrs(model1, ptr_list); - - if (verbose){ print_some( ptr_list ); } - - //-------------------------------------------------------------- - if (verbose){ puts("Calling serialize() on CFE model 1 ..."); } - - //------------------------------------------------ - // Serialize Model1 state and save to: ser_file - //------------------------------------------------ - serialize( model1, ser_file ); - - //-------------------------------------------------------------- - if (verbose){ - puts("Calling deserialize_to_state() on CFE model 2 ..."); - } - - //----------------------------------------------- - // Deserialize Model1 state saved in "ser_file" - // and set it as the new state of Model2 - //----------------------------------------------- - deserialize_to_state( ser_file, model2, print_obj ); - - //-------------------------------------------------------------- - if (verbose){ - puts("Updating BMI CFE model 2 ..."); - printf("n_steps2 = %i \n", n_steps2); - puts(""); - } - - for (int i=1; i<=n_steps2; i++){ - model2->update(model2); - } - - //-------------------------------------------------------------- - if (verbose){ - puts("Updating BMI CFE model 1 ..."); - printf("n_steps2 = %i \n", n_steps2); - puts(""); - } - - for (int i=1; i<=n_steps2; i++){ - model1->update(model1); - } - - //-------------------------------------------------------------- - if (verbose){ - puts("Calling BMI.get_state_var_ptrs() on CFE model 1 ..."); - } - - model1->get_state_var_ptrs(model1, ptr_list); - - if (verbose){ print_some( ptr_list ); } - - //-------------------------------------------------------------- - if (verbose){ - puts("Calling BMI.get_state_var_ptrs() on CFE model 2 ..."); - } - - model2->get_state_var_ptrs(model2, ptr_list); - - if (verbose){ print_some( ptr_list ); } - - //-------------------------------------------------------------- - if (verbose){ puts("Comparing CFE model 1 & 2 state vars ..."); } - - compare_states( model1, model2 ); - - //-------------------------------------------------------------- - if (verbose){ puts("Finalizing BMI CFE models 1 & 2 ..."); } - - model1->finalize(model1); - model2->finalize(model2); - - if (verbose){ - puts("Finished with serialization test.\n"); - } - return 0; -} - -//------------------------------------------------------------------------ - - - - - diff --git a/test_serialize/make_and_run_ser_test.sh b/test_serialize/make_and_run_ser_test.sh deleted file mode 100644 index d1099dcf..00000000 --- a/test_serialize/make_and_run_ser_test.sh +++ /dev/null @@ -1,10 +0,0 @@ -#!/bin/bash -FILE=./run_ser_test -if test -f "$FILE"; then - # If compilation fails, don't want to use old one - rm "$FILE" -fi - -gcc ./cfe_serialize_test.c ./serialize_state.c ../src/bmi_cfe.c ../src/cfe.c -o run_ser_test -lm -lmsgpackc -cd .. -./test_serialize/run_ser_test diff --git a/test_serialize/serialize_state.c b/test_serialize/serialize_state.c deleted file mode 100644 index f40eed41..00000000 --- a/test_serialize/serialize_state.c +++ /dev/null @@ -1,738 +0,0 @@ -//----------------------------------------- -// Notes: fbuffer.h is not included in: -// /usr/local/include/msgpack.h -//----------------------------------------- -#include -#include -#include // for strchr() -#include // to get filesize -#include -#include "msgpack/fbuffer.h" // See note above -#include "../include/bmi.h" - -//------------------------------------------------ -// Now computing these from filesize of ser_file -//------------------------------------------------ -// static const int BUFFER_SIZE = 65536; -// static const int UNPACKED_BUFFER_SIZE = 131072; -// static const int BUFFER_SIZE = 1024; -// static const int UNPACKED_BUFFER_SIZE = 2048; - -/* -These functions are meant to provide a "framework utility" -that can serialize, save to file, and then deserialize all -of the state variables of a model written in C. - -The function "serialize()" uses new BMI functions to retrieve -the model's state variables called: get_state_var_ptrs(), -get_state_var_names(), get_state_var_types(), and -get_state_var_sizes(). - -The function "deserialize_to_state()" deserializes the saved model -(from a file) and then uses the new BMI function set_state_var() -to set all state variable values into a new instance of the model. -This new instance could be running on another node. - -The function "compare_states()" compares all of the state -variables of two instances of the given model, to confirm that -they are identical after calling "deserialize_to_state()". - -These new functions could be used for the problem of -load balancing or for recovery after hardware failures. - -See this documentation for msgpack: -https://github.com/msgpack/msgpack-c/wiki/v2_0_c_overview - -See this msgpack example: -https://blog.gypsyengineer.com/en/security/msgpack-fuzzing.html -*/ - -//----------------------------------------------------------------------- -int get_file_size(const char *ser_file, unsigned long int *file_size){ - - //------------------------------------- - // Get the file size, set BUFFER_SIZE - //------------------------------------- - struct stat st; - stat(ser_file, &st); - *file_size = st.st_size; - return 0; -} - -//----------------------------------------------------------------------- -// int get_buffer_size(Bmi* model1, unsigned int buffer_size){ -// -// //--------------------------------------------------------- -// // We don't need this now. Can just use get_file_size(). -// //--------------------------------------------------------- -// int n_state_vars; -// model1->get_state_var_count(model1, &n_state_vars); -// -// unsigned int size, sizes[ n_state_vars ], type_size; -// unsigned int unpacked_buffer_size = 0; -// int verbose = 1; -// int i; -// -// char *type; -// char **types = NULL; -// types = (char**) malloc (sizeof(char *) * n_state_vars); -// for (i=0; iget_state_var_types(model1, types); -// -// if (verbose){ puts("Calling BMI.get_state_var_sizes()..."); } -// model1->get_state_var_sizes(model1, sizes); -// -// for (i=0; iget_state_var_count(model1, &n_state_vars); - - FILE *fp = fopen(ser_file, "w+"); - int i, j, n_bytes; - void *ptr_list[ n_state_vars ]; - unsigned int size, sizes[ n_state_vars ]; - int verbose = 1; - int i_val; - long li_val; - float f_val; - double d_val; - - char **names = NULL; - names = (char**) malloc (sizeof(char *) * n_state_vars); - for (i=0; iget_state_var_names(model1, names); - - if (verbose){ puts("Calling BMI.get_state_var_types()..."); } - model1->get_state_var_types(model1, types); - - if (verbose){ puts("Calling BMI.get_state_var_sizes()..."); } - model1->get_state_var_sizes(model1, sizes); - - if (verbose){ puts("Calling BMI.get_state_var_ptrs()..."); } - model1->get_state_var_ptrs(model1, ptr_list); - - //-------------------------------------------- - // Prepare to write serialized state to file - // msgpack_fbuffer_write needs fbuffer.h - //-------------------------------------------- - msgpack_packer pk; - msgpack_packer_init(&pk, fp, msgpack_fbuffer_write); - - //---------------------------------------------- - // Prepare to write serialized state to buffer - //---------------------------------------------- - // msgpack_sbuffer* buffer = msgpack_sbuffer_new(); - // msgpack_packer* pk = msgpack_packer_new(buffer, msgpack_sbuffer_write); - - //------------------------------------------------------ - // Note: strings are serialized with 2 commands, like: - //------------------------------------------------------ - // msgpack_pack_str(&pk, 7); - // msgpack_pack_str_body(&pk, "example", 7); - - //-------------------------------------- - // Note: booleans are serialized like: - //-------------------------------------- - // msgpack_pack_true(&pk); - // msgpack_pack_false(&pk); - - //---------------------------------------------- - // Dereference the pointers & serialize values - // Note that we aren't using the names here. - //---------------------------------------------- - if (verbose){ - puts("Serializing the model state variables..."); - printf(" Number of state vars = %d\n", n_state_vars); - } - for (i=0; iget_state_var_count(model2, &n_state_vars); - - unsigned int sizes[ n_state_vars ], size; - model2->get_state_var_sizes(model2, sizes); - int i_val, *i_arr, j; - long li_val, *li_arr; - float f_val, *f_arr; - double d_val, *d_arr; - void *ptr; - - char *type, *sval; - char **types = NULL; - types = (char**) malloc (sizeof(char *) * n_state_vars); - for (int j=0; jget_state_var_types(model2, types); - - //------------------------------------- - // Get the file size, set buffer_size - //------------------------------------- - unsigned long int file_size, buffer_size, unpacked_buffer_size; - get_file_size( ser_file, &file_size ); - buffer_size = file_size; - unpacked_buffer_size = 2 * buffer_size; - char inbuffer[buffer_size]; - char unpacked_buffer[unpacked_buffer_size]; -// if (verbose){ -// printf("Buffer_size = %lu\n", buffer_size); -// printf("Unpacked buffer_size = %lu\n", unpacked_buffer_size); -// } - // char inbuffer[BUFFER_SIZE]; - // char unpacked_buffer[UNPACKED_BUFFER_SIZE]; - - FILE *fp = fopen(ser_file, "rb"); - int i = 0; - size_t off = 0; - size_t len = 0; - - msgpack_unpacked unpacked; - msgpack_unpack_return ret; - msgpack_unpacked_init(&unpacked); - if (verbose){ puts("Deserializing source model state vars..."); } - - //------------------------------------------------------------ - // In online ref 1, buffer is passed as an argument (buf). - // In online ref 2, buffer is read from ser_file (inbuffer). - //------------------------------------------------------------ - // In online ref 1, len is passed as an argument. - // In online ref 2, len=read is return value from fread. - //------------------------------------------------------------ - // Online ref 1: "result" = Online ref 2: "unpacked". - // Online ref 1: "buf" = Online ref 2: "inbuffer" ???? - // Online ref 1: "len" = Online ref 2: "read" - //------------------------------------------------------------ - // len = fread(inbuffer, sizeof(char), BUFFER_SIZE, fp); - len = fread(inbuffer, sizeof(char), buffer_size, fp); - ret = msgpack_unpack_next(&unpacked, inbuffer, len, &off); - - while (ret == MSGPACK_UNPACK_SUCCESS) { - msgpack_object obj = unpacked.data; - - size = sizes[i]; - type = types[i]; - // printf("type = %s\n", type); - if (strchr(type, '*') != NULL){ - type[strlen(type)-1] = '\0'; // remove last char - } - - if (size == 1){ - //-------------------------------------------- - // Note: This does not work, even though we - // typecast in set_state_var(). - //-------------------------------------------- - // ptr = &obj; - // model2->set_state_var(model2, ptr, i ); - //-------------------------------------------- - if (strcmp(type, "int") == 0){ - i_val = (int)obj.via.i64; - model2->set_state_var(model2, &i_val, i ); - } else if (strcmp(type, "long") == 0){ - li_val = (long)obj.via.i64; - model2->set_state_var(model2, &li_val, i ); - } else if (strcmp(type, "float") == 0){ - f_val = (float)obj.via.f64; - model2->set_state_var(model2, &f_val, i ); - } else if (strcmp(type, "double") == 0){ - d_val = (double)obj.via.f64; - model2->set_state_var(model2, &d_val, i ); - } else if (strcmp(type, "string") == 0){ - model2->set_state_var(model2, &obj, i ); - //------------------------------------------- - // Next 2 lines don't work - // sval = (char*)(obj.via.str); - // model2->set_state_var(model2, sval, i ); - //------------------------------------------- - // Next 2 lines don't work either - //sval = obj.via.str; - //model2->set_state_var(model2, sval, i ); - } else{ - model2->set_state_var(model2, &obj, i ); - } - } else{ - //-------------------------------------------------- - // Copy values from msgpack array object into an - // array of the correct type, then pass pointer to - // that array to set_state_var(). (This works.) - //-------------------------------------------------- - // Couldn't figure out a clean way to access the - // msgpack array elements with pointer math. This - // seems to be because it allows mixed types. - // Options: via.u64, via.i64, via.f64, via.str. - //-------------------------------------------------- - if (strcmp(type, "string") == 0){ - // Works for title & subcat strings in TOPMODEL. - ptr = obj.via.array.ptr; - model2->set_state_var(model2, ptr, i ); - } else if (strcmp(type, "int") == 0){ - i_arr = (int*) malloc(size * sizeof( int )); - for (int j=0; jset_state_var(model2, i_arr, i ); - free(i_arr); - } else if (strcmp(type, "long") == 0){ - li_arr = (long*) malloc(size * sizeof( long )); //########## - for (int j=0; jset_state_var(model2, li_arr, i ); - free(li_arr); - } else if (strcmp(type, "float") == 0){ - f_arr = (float*) malloc(size * sizeof( float )); - for (int j=0; jset_state_var(model2, f_arr, i ); - free(f_arr); - } else if (strcmp(type, "double") == 0){ - d_arr = (double*) malloc(size * sizeof( double )); - for (int j=0; jset_state_var(model2, d_arr, i ); - free(d_arr); - } - - // For testing - //printf("### obj.via.array.ptr[0].via.f64 = %f\n", obj.via.array.ptr[0].via.f64); - //printf("### obj.via.array.ptr[1].via.f64 = %f\n", obj.via.array.ptr[1].via.f64); - //printf("### obj.via.array.ptr[2].via.f64 = %f\n", obj.via.array.ptr[2].via.f64); - //printf("\n"); - } - - //------------------------------------------------------------ - // Failed attempts to access a single pointer to an array in - // the msgpack_array structure, to pass to set_state_var(). - // Note that typecasting occurs in set_state_var(). - //------------------------------------------------------------ -// if (strcmp(type, "double") == 0){ -// ptr = obj.via.array.ptr; -// model2->set_state_var(model2, ptr, i ); } -// //----------------------------------------------- -// if (strcmp(type, "double") == 0){ -// ptr = &obj; -// model2->set_state_var(model2, ptr, i ); } -// //----------------------------------------------- -// if (strcmp(type, "double") == 0){ -// ptr = obj.via.array.ptr[0]; -// model2->set_state_var(model2, ptr, i ); } -// //----------------------------------------------- -// if (strcmp(type, "double") == 0){ -// d_val = obj.via.array.ptr[0].via.f64; -// ptr = &d_val; -// model2->set_state_var(model2, ptr, i ); } -// //----------------------------------------------- - - if (print_obj){ - printf("Object no %d:\n", i); - msgpack_object_print(stdout, obj); - printf("\n"); - - //------------------------------ - // This just prints obj again - //------------------------------ - // msgpack_object_print_buffer(unpacked_buffer, UNPACKED_BUFFER_SIZE, obj); - // printf("%s\n", unpacked_buffer); - } - //----------------------------------------------------- - i++; - ret = msgpack_unpack_next(&unpacked, inbuffer, len, &off); - } // end of while loop - - //-------------------------------------------------------- - // Did we unpack the expected number of state variables? - //-------------------------------------------------------- - if (i < n_state_vars){ - printf("WARNING: Expected %d state variables \n", n_state_vars); - printf(" But unpacked only %d vars. \n", i); - printf(" BUFFER_SIZE may be too small."); - printf(""); - } else{ - printf("Unpacked %d state variables.\n", i); - printf(""); - } - - //if (ret == MSGPACK_UNPACK_CONTINUE) { - // printf("Every msgpack_object in the buffer was processed.\n"); - //} - //else if (ret == MSGPACK_UNPACK_PARSE_ERROR) { - // printf("The data in the buffer has an invalid format.\n"); - //} - - msgpack_unpacked_destroy(&unpacked); - fclose(fp); // Try moving this up to just after fread. - - //-------------------------------- - // Free memory for string arrays - //-------------------------------- - if (verbose){ puts("Freeing memory...");} - for (i=0; iget_state_var_count(model1, &n_state_vars); - - //---------------------------------------- - // Get the state variable internal names - //---------------------------------------- - char **names = NULL; - names = (char**) malloc (sizeof(char *) * n_state_vars); - for (i=0; iget_state_var_names(model1, names); - - //------------------------------------ - // Get the state variable data types - //------------------------------------ - char *type; - char **types = NULL; - types = (char**) malloc (sizeof(char *) * n_state_vars); - for (i=0; iget_state_var_types(model1, types); - - //------------------------------- - // Get the state variable sizes - //------------------------------- - unsigned int size, sizes[ n_state_vars ]; - model1->get_state_var_sizes(model1, sizes); - - //------------------------------------- - // Get pointers to Model 1 state vars - //------------------------------------- - void *ptr_list1[ n_state_vars ]; - model1->get_state_var_ptrs(model1, ptr_list1); - - //------------------------------------- - // Get pointers to Model 2 state vars - //------------------------------------- - void *ptr_list2[ n_state_vars ]; - model2->get_state_var_ptrs(model2, ptr_list2); - - //-------------------------------- - // Loop over all state variables - //-------------------------------- - for (int i=0; i