Adding Xinanjiang direct runoff option, changed volume trackers to st…

…ructure in CFE(#25)

Also fixes update until for BMI 2.0 to work appropriately.

Co-authored-by: madMatchstick <[email protected]>

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:

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

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
+verbosity=1
+surface_partitioning_scheme=Schaake

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

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
+verbosity=2
+surface_partitioning_scheme=Xinanjiang
+a_Xinanjiang_inflection_point_parameter=1
+b_Xinanjiang_shape_parameter=1
+x_Xinanjiang_shape_parameter=1

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
+verbosity=1
+surface_partitioning_scheme=Xinanjiang
+a_Xinanjiang_inflection_point_parameter=1
+b_Xinanjiang_shape_parameter=1
+x_Xinanjiang_shape_parameter=1

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

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

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

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(&copy, ",");
-    // 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(&copy, ",");
-    // Not sure what this is
-
-    // DLWRF_surface
-    value = strsep(&copy, ",");
-    aorc->incoming_longwave_W_per_m2 = strtof(value, NULL);
-    // DSWRF_surface
-    value = strsep(&copy, ",");
-    aorc->incoming_shortwave_W_per_m2 = strtof(value, NULL);
-    // PRES_surface
-    value = strsep(&copy, ",");
-    aorc->surface_pressure_Pa = strtof(value, NULL);
-    // SPFH_2maboveground
-    value = strsep(&copy, ",");
-    aorc->specific_humidity_2m_kg_per_kg = strtof(value, NULL);;
-    // TMP_2maboveground
-    value = strsep(&copy, ",");
-    aorc->air_temperature_2m_K = strtof(value, NULL);
-    // UGRD_10maboveground
-    value = strsep(&copy, ",");
-    aorc->u_wind_speed_10m_m_per_s = strtof(value, NULL);
-    // VGRD_10maboveground
-    value = strsep(&copy, ",");
-    aorc->v_wind_speed_10m_m_per_s = strtof(value, NULL);
-    // precip_rate
-    value = strsep(&copy, ",");
-    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 ############################*/

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;

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);
 

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;