Fix dynamic alloc

data/topmod.run

@@ -4,4 +4,4 @@ data/inputs.dat
 data/subcat.dat
 data/params.dat
 topmod.out
-hyd.out
+hyd.out

include/topmodel.h

@@ -17,24 +17,51 @@
 // Note: All errors causing program to exit will print console message
 #define TOPMODEL_DEBUG 0
 
+//The following "limits" from the original beven code
+//are now treated as warnings
+#define WARN_NUM_SUBCATCHMENTS 10
+#define WARN_HISTOGRAM_ORDINATES 20
+#define WARN_TOPODEX_INCREMENTS 30
+
 /*** Function/subroutine prototypes ***/
-extern void init(FILE *in_param_fptr,FILE *output_fptr,char *subcat, 
-                int num_channels,int num_topodex_values,
-                int yes_print_output, double area,double **time_delay_histogram,
-                double *cum_dist_area_with_dist,double dt, double *szm, double *t0,
-                double tl, double *dist_from_outlet,double *td, double *srmax, 
-                double *Q0,double *sr0, int *infex, double *xk0, double *hf, 
-                double *dth,double **stor_unsat_zone,double **deficit_local,
-                double **deficit_root_zone,double *szq, double *Q,double *sbar,
-                int max_atb_increments, int max_time_delay_ordinates,
-                double *bal,int *num_time_delay_histo_ords, int *num_delay);
+extern void convert_dist_to_histords(const double * const dist_from_outlet, const int num_channels,
+					const double * const chv, const double * const rv, const double dt, double* const tch);
+
+extern void calc_time_delay_histogram(int num_channels, double area,
+		double* tch, double *cum_dist_area_with_dist,
+		int *num_time_delay_histo_ords, int *num_delay,	double **time_delay_histogram);
+
+extern void init_water_balance(
+				int num_topodex_values, double dt, double *sr0, 
+				double *szm, double *Q0, double *t0, double tl,
+				double **stor_unsat_zone, double *szq,
+				double **deficit_local, double **deficit_root_zone, 
+				double *sbar, double *bal);
+
+extern void init_discharge_array(int stand_alone, int *num_delay, double *Q0, double area, 
+			int *num_time_delay_histo_ords, double **time_delay_histogram,
+                        double **Q);
+
+extern void init(FILE *in_param_fptr, FILE *output_fptr, char *subcat, int stand_alone,
+	      int num_channels, int num_topodex_values, int yes_print_output,
+	      double area, double **time_delay_histogram,
+	      double *cum_dist_area_with_dist, double dt, 
+        double tl, double *dist_from_outlet, 
+	      int *num_time_delay_histo_ords,int *num_delay,
+	      double *szm, double *t0, double *chv, double *rv, double *td, double *srmax, 
+        double *Q0,double *sr0, int *infex, double *xk0, double *hf, double *dth,
+	      double **stor_unsat_zone, double **deficit_local,
+        double **deficit_root_zone,double *szq,double **Q,
+        double *sbar, double *bal);
+
+
 
 extern void inputs(FILE *input_fptr, int *nstep, double *dt, double **rain,
                 double **pe, double **Qobs, double **Q, double **contrib_area);
 
 extern void topmod(FILE *output_fptr, int nstep, int num_topodex_values,
                 int yes_print_output,int infex, double dt, double szm,
-                double *stor_unsat_zone, double *deficit_root_zone,
+	        double *stor_unsat_zone, double *deficit_root_zone,
                 double *deficit_local, double *pe, double *rain,double xk0,double hf, 
                 double *dist_area_lnaotb, double tl, double *lnaotb, double td,
                 double srmax, double *contrib_area, double szq, double *Qout, 
@@ -48,7 +75,7 @@ extern void tread(FILE *subcat_fptr,FILE *output_fptr,char *subcat,
                 int *num_topodex_values,int *num_channels,double *area,
                 double **dist_area_lnaotb, double **lnaotb, int yes_print_output,
                 double **cum_dist_area_with_dist, double *tl, 
-                double **dist_from_outlet, int maxsubcatch, int maxincr);
+                double **dist_from_outlet);
 
 extern void expinf(int irof, int it, int rint, double *df, double *cumf,
                 double dt,double xk0, double szm, double hf);                  
@@ -106,6 +133,8 @@ struct TopModel_Struct{
   double hf;    /* wetting front suction for G&A soln.  */
   double dth;   /* water content change across the wetting front */
   double area;  /* catchment area */
+  double chv;   /* average channel flow velocity */
+  double rv;    /* internal overland flow routing velocity */
 
   /************ Other variables of note ************/
   int num_delay;          /* number of time steps lag (delay) in channel within catchment to outlet */
@@ -133,9 +162,7 @@ struct TopModel_Struct{
 
   /************** Other internal vars **************/
   int num_sub_catchments;
-  int max_atb_increments;
-  int max_num_subcatchments;
-  int max_time_delay_ordinates;
+
   // int isc   //subcat loop no longer handled by topmod()
 
   /******************** Returns ********************/

src/bmi_topmodel.c

@@ -2,11 +2,12 @@
 #include "../include/bmi.h" 
 #include "../include/bmi_topmodel.h"
 
+
 /* BMI Adaption: Max i/o file name length changed from 30 to 256 */
 #define MAX_FILENAME_LENGTH 256
 #define OUTPUT_VAR_NAME_COUNT 14
 #define INPUT_VAR_NAME_COUNT 2
-#define PARAM_VAR_NAME_COUNT 6
+#define PARAM_VAR_NAME_COUNT 8
 
 static const char *output_var_names[OUTPUT_VAR_NAME_COUNT] = {
         "Qout",
@@ -149,7 +150,9 @@ static const char *param_var_names[PARAM_VAR_NAME_COUNT] = {
     "td",    // unsaturated zone time delay per unit storage deficit (h)
     "srmax", // maximum root zone storage deficit (m)
     "sr0",   // initial root zone storage deficit below field capacity (m)
-    "xk0"    // surface soil hydraulic conductivity (m/h)
+    "xk0",   // surface soil hydraulic conductivity (m/h)
+    "chv",   // average channel velcoity
+    "rv"     // internal overland flow routing velocity
 };
 
 static const char *param_var_types[PARAM_VAR_NAME_COUNT] = {
@@ -158,6 +161,8 @@ static const char *param_var_types[PARAM_VAR_NAME_COUNT] = {
     "double",
     "double",
     "double",
+    "double",
+    "double",
     "double"
 };
 
@@ -279,33 +284,32 @@ int init_config(const char* config_file, topmodel_model* model)
         d_alloc(&model->rain,model->nstep);
         d_alloc(&model->pe,model->nstep);
         d_alloc(&model->Qobs,model->nstep);   //TODO: Consider removing this all together when framework
+        //FIXME this is no longer useful, as the follow `init` will ultimatey
+        //reallocate Q to based on the number of computed histogram ordinates
         d_alloc(&model->Q,model->nstep);
         d_alloc(&model->contrib_area,model->nstep);
 
         (model->rain)[1]=0.0;
         (model->pe)[1]=0.0;
         (model->Qobs)[1]=0.0;
+        //FIXME not really useful, as `init` will reallocate and properly initialize
         (model->Q)[1]=0.0;
         (model->contrib_area)[1]=0.0;
     }
 
-    // Set up maxes for subcat and params read-in functions
-    model-> max_atb_increments=30;
-    model-> max_num_subcatchments=10;
-    model-> max_time_delay_ordinates=20;
-
     tread(model->subcat_fptr,model->output_fptr,model->subcat,&model->num_topodex_values,&model->num_channels,
         &model->area,&model->dist_area_lnaotb,&model->lnaotb,model->yes_print_output,
-        &model->cum_dist_area_with_dist,&model->tl,&model->dist_from_outlet,
-        model->max_num_subcatchments,model->max_atb_increments);
-    fclose(model->subcat_fptr);
+        &model->cum_dist_area_with_dist,&model->tl,&model->dist_from_outlet);
 
-    init(model->params_fptr,model->output_fptr,model->subcat,model->num_channels,model->num_topodex_values,
+    fclose(model->subcat_fptr);
+
+    init(model->params_fptr,model->output_fptr,model->subcat,model->stand_alone, model->num_channels, model->num_topodex_values,
         model->yes_print_output,model->area,&model->time_delay_histogram,model->cum_dist_area_with_dist,
-        model->dt,&model->szm,&model->t0,model->tl,model->dist_from_outlet,&model->td, &model->srmax,&model->Q0,&model->sr0,&model->infex,&model->xk0,&model->hf,
-        &model->dth,&model->stor_unsat_zone,&model->deficit_local,&model->deficit_root_zone,
-        &model->szq,model->Q,&model->sbar,model->max_atb_increments,model->max_time_delay_ordinates,
-        &model->bal,&model->num_time_delay_histo_ords,&model->num_delay);
+        model->dt,model->tl,model->dist_from_outlet,&model->num_time_delay_histo_ords,&model->num_delay,
+	&model->szm,&model->t0,&model->chv,&model->rv,&model->td, &model->srmax,
+	&model->Q0,&model->sr0,&model->infex,&model->xk0,&model->hf,&model->dth,
+	&model->stor_unsat_zone,&model->deficit_local,&model->deficit_root_zone,
+        &model->szq,&model->Q,&model->sbar, &model->bal);
     fclose(model->params_fptr);
 
 
@@ -430,7 +434,7 @@ static int Update (Bmi *self)
         &topmodel->sump,&topmodel->sumae,&topmodel->sumq,&topmodel->sumrz,&topmodel->sumuz,
         &topmodel->quz, &topmodel->qb, &topmodel->qof, &topmodel->p, &topmodel->ep );
 
-    if ((topmodel->stand_alone == FALSE) & (topmodel->yes_print_output == TRUE)){
+    if ((topmodel->stand_alone == FALSE) && (topmodel->yes_print_output == TRUE)){
         fprintf(topmodel->out_hyd_fptr,"%d %lf %lf\n",topmodel->current_time_step,topmodel->Qobs[1],topmodel->Q[1]);
     }      
 
@@ -468,6 +472,7 @@ static int Update_until (Bmi *self, double t)
 
 static int Finalize (Bmi *self)
 {
+
     if (self){
         topmodel_model* model = (topmodel_model *)(self->data);
 
@@ -489,9 +494,9 @@ static int Finalize (Bmi *self)
             // framework-controlled mode.
             //-----------------------------------------------------------
             if (model->stand_alone == TRUE){
-                results(model->output_fptr,model->out_hyd_fptr,model->nstep, 
+	        results(model->output_fptr,model->out_hyd_fptr,model->nstep, 
                 model->Qobs, model->Q, model->yes_print_output);                 
-            }
+	    }
         }    
 
         if( model->Q != NULL )
@@ -840,6 +845,22 @@ static int Get_value_ptr (Bmi *self, const char *name, void **dest)
         *dest = (void*)&topmodel-> t0;
         return BMI_SUCCESS;
     }
+    // chv (parameter)
+    if (strcmp (name, "chv") == 0) {
+        topmodel_model *topmodel;
+        topmodel = (topmodel_model *) self->data;
+        *dest = (void*)&topmodel-> chv;
+        return BMI_SUCCESS;
+    }
+    // rv (parameter)
+    if (strcmp (name, "rv") == 0) {
+        topmodel_model *topmodel;
+        topmodel = (topmodel_model *) self->data;
+        *dest = (void*)&topmodel-> rv;
+        return BMI_SUCCESS;
+    }
+
+
 
 
     // STANDALONE Note: 
@@ -918,9 +939,143 @@ static int Set_value (Bmi *self, const char *name, void *array)
 
     memcpy (dest, array, nbytes);
 
+
+    ///////////////////////////////////
+    //Handle calibratable parameters///
+    ///////////////////////////////////
+
+
+    // CHECK IF/WHICH CALIBRATABLE PARAMETERS TO UPDATE
+
+    // DISCHARGE RELATED PARAMETERS
+
+    // define array holding calibratable parameter names
+    char *calibQParams[] = {"szm", "t0", "chv", "rv", "sr0"};
+
+    // get number of strings to use for loop
+    int numQParams = sizeof(calibQParams) / sizeof(calibQParams[0]);
+
+    // check if name is == any calibratable parameter
+    // create holder to check if name == any calibQParam
+    int nameIsCalibQParam = 0;
+    for (int i = 0; i<numQParams; i++){
+ 	if (strcmp(name, calibQParams[i]) == 0) {
+	    nameIsCalibQParam = 1;
+	    break;
+        }
+    }
+
+
+    // WATER BALANCE RELATED PARAMETERS
+
+    // define array holding calibratable parameter names
+    char *calibWBParams[] = {"szm", "t0", "sr0"};
+
+    // get number of strings to use for loop
+    int numWBParams = sizeof(calibWBParams) / sizeof(calibWBParams[0]);
+
+    // check if name is == any calibratable parameter
+    // create holder to check if name == any calibWBParam
+    int nameIsCalibWBParam = 0;
+    for (int i = 0; i<numWBParams; i++){
+ 	if (strcmp(name, calibWBParams[i]) == 0) {
+	    nameIsCalibWBParam = 1;
+	    break;
+        }
+    }
+
+    // if any calibratable variables were provided in realization file, then 
+    // print an update with the values being used.
+    // srmax and td do not need to be updated, but are printed to confirm their
+    // presence in the realization.json file.
+    if (nameIsCalibQParam || nameIsCalibWBParam || \
+		    strcmp(name, "srmax") == 0 || strcmp(name, "td") == 0) {
+
+        // instantiate topmodel as a pointer topmodel of type topmodel_model    
+        topmodel_model *topmodel;
+        // assign self->data to topmodel pointer
+        topmodel = (topmodel_model *) self->data;
+
+        printf("\n\n\nAT LEAST ONE OF THE FOLLOWING CALIBRATABLE PARAMETERS "
+			"WAS PROVIDED IN THE REALIZATION.JSON FILE!\n");
+
+	// print updated calibratable parameters
+        printf("\n\nCalibratable parameters related to ET and recharge:\n");
+	printf("srmax = %f\n", topmodel->srmax);
+        printf("td = %f\n", topmodel->td);
+
+	printf("\nCalibratable parameters related to discharge:\n");
+        printf("chv = %f\n", topmodel->chv);
+        printf("rv = %f\n", topmodel->rv);
+
+
+	printf("\nCalibratable parameters related to water balance:\n");
+        printf("szm = %f\n", topmodel->szm); 
+	printf("sr0 = %f\n", topmodel->sr0);
+        printf("t0 = %f\n\n\n\n", topmodel->t0);
+
+    }
+
+    // UPDATE APPROPRIATE PARAMETERS
+
+    // DISCHARGE RELATED PARAMETERS
+
+    // if name is a calibratable parameter, then update outputs
+    if (nameIsCalibQParam) {
+
+         // instantiate topmodel as a pointer topmodel of type topmodel_model    
+        topmodel_model *topmodel;
+        // assign self->data to topmodel pointer
+        topmodel = (topmodel_model *) self->data;
+
+        // declare variables
+        double* tch; //FIXME put this in topmod struct and reuse it
+        d_alloc(&tch, topmodel->num_channels);
+
+
+     	// convert from distance/area to histogram ordinate form
+        convert_dist_to_histords(topmodel->dist_from_outlet, topmodel->num_channels,
+				&topmodel->chv, &topmodel->rv, topmodel->dt, tch);
+
+        // calculate the time_delay_histogram
+        calc_time_delay_histogram(topmodel->num_channels, 
+	  			  topmodel->area, tch, 
+ 				  topmodel->cum_dist_area_with_dist, 
+				  &topmodel->num_time_delay_histo_ords,
+                  &topmodel->num_delay, &topmodel->time_delay_histogram);
+        free(tch);
+        // Reinitialise discharge array
+        init_discharge_array(topmodel->stand_alone, &topmodel->num_delay, &topmodel->Q0, topmodel->area, 
+				&topmodel->num_time_delay_histo_ords, &topmodel->time_delay_histogram, 
+				&topmodel->Q);	
+    }
+
+
+    // WATER BALANCE RELATED PARAMETERS
+
+    // if name is a calibratable parameter, then update outputs
+    if (nameIsCalibWBParam) {
+
+        // instantiate topmodel as a pointer topmodel of type topmodel_model    
+        topmodel_model *topmodel;
+        // assign self->data to topmodel pointer
+        topmodel = (topmodel_model *) self->data;
+
+
+        // Initialize water balance and unsatrutaed storage and deficits
+        init_water_balance(topmodel->num_topodex_values, 
+					topmodel->dt, &topmodel->sr0, &topmodel->szm, 
+					&topmodel->Q0, &topmodel->t0, topmodel->tl,
+					&topmodel->stor_unsat_zone, &topmodel->szq, 
+					&topmodel->deficit_local, &topmodel->deficit_root_zone, 
+					&topmodel->sbar, &topmodel->bal);            
+    }
+
     return BMI_SUCCESS;
 }
 
+
+
 static int Set_value_at_indices (Bmi *self, const char *name, int * inds, int len, void *src)
 {
     void * to = NULL;

src/main.c

@@ -48,6 +48,13 @@ int main(void)
     // Gather number of steps from input file
     // when in standalone mode.
     n = topmodel->nstep;
+    //NJF in standalone mode, one MUST allocate Q for all time steps
+    //This should probably re-considered and additional validation done
+    //that the algorithms work as intended
+    //especially with the addition of `shift_Q` in the topmodel function
+    //but for now, this prevents a fault
+    if(topmodel->Q != NULL) free(topmodel->Q);
+    d_alloc(&topmodel->Q, n);
   }
   else{
     // Otherwise define loop here