diff --git a/data/topmod.run b/data/topmod.run index f730eac..1565470 100644 --- a/data/topmod.run +++ b/data/topmod.run @@ -4,4 +4,4 @@ data/inputs.dat data/subcat.dat data/params.dat topmod.out -hyd.out \ No newline at end of file +hyd.out diff --git a/include/topmodel.h b/include/topmodel.h index 4ba300c..fc8815b 100755 --- a/include/topmodel.h +++ b/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 ********************/ diff --git a/src/bmi_topmodel.c b/src/bmi_topmodel.c index d31d05d..a41bdf3 100755 --- a/src/bmi_topmodel.c +++ b/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; idata 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; diff --git a/src/main.c b/src/main.c index 954a46d..e7d15d0 100755 --- a/src/main.c +++ b/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 diff --git a/src/topmodel.c b/src/topmodel.c index 27e8c2c..0bef0c4 100755 --- a/src/topmodel.c +++ b/src/topmodel.c @@ -99,6 +99,14 @@ *****************************************************************/ #include "../include/topmodel.h" +void shift_Q(double* Q, int num_ordinates){ //NJF Really don't like the +1 size assumption here + //memmove puts num_ordinates*sizeof(double) bytes starting at Q[1] + //into Q[0], effectively shifting the array by one + memmove(&Q[0], &Q[1], num_ordinates*sizeof(*Q)); + //set the last value to 0.0 + Q[num_ordinates] = 0.0; +} + 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, @@ -135,7 +143,13 @@ extern void topmod(FILE *output_fptr, int nstep, int num_topodex_values, current_time_step, *sump, *sumae, *sumq, stand_alone added as function input parameters */ -double ex[31]; +double ex[num_topodex_values+1]; //+1 to maintin 1 based array indexing +//NJF TODO consider warning on all program limits here since this is essentially +//"the model" where those assumptions may not be valid... +if(num_topodex_values > WARN_TOPODEX_INCREMENTS){ + printf("WARNING: num_topodex_values, %d, is greater than %d\n", + num_topodex_values, WARN_TOPODEX_INCREMENTS); +} int ia,ib,in,irof,it,ir; double rex,cumf,max_contrib_area,ea,rint,acm,df; double acf,uz,sae,of; @@ -159,8 +173,15 @@ if(yes_print_output==TRUE && current_time_step==1) /* BMI Adaption: Set iteration to bmi's current_time_step (standalone) or 1 (framework) Counter++ is handled by bmi's update()*/ -if (stand_alone == TRUE)it=current_time_step; -else it=1; +if (stand_alone == TRUE) { + it=current_time_step; +} +else { + it=1; + //shift Q array to align current time step + shift_Q(Q, num_delay + num_time_delay_histo_ords); +} + *qof=0.0; *quz=0.0; @@ -300,7 +321,7 @@ if(irof==1) max_contrib_area=1.0; (*qb)=szq*exp(-(*sbar)/szm); (*sbar)+=(-(*quz)+(*qb)); *Qout=(*qb)+(*qof); -*sumq+=(*Qout); /* BMI Adaption: *sumq now as pointer var; incl in model struct */ + /* CHANNEL ROUTING CALCULATIONS */ /* allow for time delay to catchment outlet num_delay as well as */ @@ -312,9 +333,12 @@ for(ir=1;ir<=num_time_delay_histo_ords;ir++) { in=it+num_delay+ir-1; if(in>current_time_step) break; - Q[in]=(*Qout)*time_delay_histogram[ir]; + //Accumulate previous time dealyed flow with current + Q[in]+=(*Qout)*time_delay_histogram[ir]; } + //Add current time flow to mass balance variable + *sumq += Q[it]; /* BMI Adaption: replace nstep with current_time_step */ if(yes_print_output==TRUE && in<=current_time_step) { @@ -394,7 +418,11 @@ This change should be make by the AGU 2021 CONUS scale demonstration*/ d_alloc(r,(*nstep)); d_alloc(pe,(*nstep)); d_alloc(Qobs,(*nstep)); -d_alloc(Q,(*nstep)); +//NJF This is dangerous, there is no validation between nstep and +//number of historgram ordinates, which is used to iterate Q in later steps +//so this could easily overflow if nstep is smaller than historgram ords +d_alloc(Q,(*nstep)); //NJF TODO validate that this works correctly + //When used in "stand alone" mode d_alloc(contrib_area,(*nstep)); //--------------------------------------- @@ -414,8 +442,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 max_num_subcatchments, - int max_num_increments) + double **dist_from_outlet) { /************************************************************** @@ -426,19 +453,8 @@ extern void tread(FILE *subcat_fptr,FILE *output_fptr,char *subcat, double tarea,sumac; int j; - -if((*cum_dist_area_with_dist)==NULL) - { - d_alloc(cum_dist_area_with_dist,max_num_subcatchments+1); - d_alloc(dist_from_outlet,max_num_subcatchments+1); - } -if((*dist_area_lnaotb)==NULL) - { - d_alloc(dist_area_lnaotb,max_num_increments+1); - d_alloc(lnaotb,max_num_increments+1); - } - - +// NJF This won't work unless the first line of subcat file has already been consumed +// Should probably document this behavior for this function... fgets(subcat,256,subcat_fptr); /* do twice to read in line-feed */ fgets(subcat,256,subcat_fptr); @@ -448,30 +464,37 @@ if (yes_print_output == TRUE) } fscanf(subcat_fptr,"%d %lf",num_topodex_values,area); - -if((*num_topodex_values)>max_num_increments) - { - printf("Error, you have too many increments in your ln(a/tanB) file.\n"); - printf(" The maximum is: %d\n",max_num_increments); - printf("program stopped.\n"); - exit(-9); - } +//Setup the topoindex arrays +if(*num_topodex_values > WARN_TOPODEX_INCREMENTS){ + printf("WARNING: Number of ln(a/tanB) increments,%d, is > than %d\n", + *num_topodex_values, + WARN_TOPODEX_INCREMENTS); +} +if((*dist_area_lnaotb)==NULL) +{ + //Need one extra value (0.0) at the end of this array to do + //sum in topmod, e.g. dist_area_lnaotb[i] + dist_area_lnaotb[i+1] + //for each increment + d_alloc(dist_area_lnaotb, *num_topodex_values+1); +} +if((*lnaotb) == NULL){ + d_alloc(lnaotb, *num_topodex_values); //NJF why +1 in old code? +} +//NJF TODO if not NULL should probably free and then allocate? +//I don't see any reasonable expectation that these arrays should persist +//across calls to tread... /* num_topodex_values IS NUMBER OF A/TANB ORDINATES */ /* AREA IS SUBCATCHMENT AREA AS PROPORTION OF TOTAL CATCHMENT */ +tarea = 0; //Accumulate area as it is read for(j=1;j<=(*num_topodex_values);j++) { fscanf(subcat_fptr,"%lf %lf",&(*dist_area_lnaotb)[j],&(*lnaotb)[j]); + tarea += (*dist_area_lnaotb)[j]; } /* dist_area_lnaotb IS DISTRIBUTION OF AREA WITH LN(A/TANB) */ /* lnaotb IS LN(A/TANB) VALUE */ -tarea = (*dist_area_lnaotb)[1]; -for(j=2;j<=(*num_topodex_values);j++) - { - tarea+=(*dist_area_lnaotb)[j]; - } - /* CALCULATE AREAL INTEGRAL OF LN(A/TANB) * NB. a/tanB values should be ordered from high to low with lnaotb[1] * as an upper limit such that dist_area_lnaotb[1] should be zero, with dist_area_lnaotb(2) representing @@ -486,10 +509,28 @@ for(j=2;j<=(*num_topodex_values);j++) sumac+=(*dist_area_lnaotb)[j]; (*tl)+=(*dist_area_lnaotb)[j]*((*lnaotb)[j]+(*lnaotb)[j-1])/2.0; } +//init last value to additive identity (*dist_area_lnaotb)[(*num_topodex_values)+1]=0.0; /* READ CHANNEL NETWORK DATA */ fscanf(subcat_fptr,"%d",num_channels); + +if(*num_channels > WARN_NUM_SUBCATCHMENTS){ + printf("WARNING: Number of channels, %d, is greater than %d\n", + *num_channels, + WARN_NUM_SUBCATCHMENTS); +} + +if((*cum_dist_area_with_dist)==NULL){ + d_alloc(cum_dist_area_with_dist, *num_channels); //NJF why +1 in old code? +} +if((*dist_from_outlet) == NULL){ + d_alloc(dist_from_outlet, *num_channels); //NJF why +1 in old code? +} +//NJF TODO if not NULL should probably free and then allocate? +//I don't see any reasonable expectation that these arrays should persist +//calls to tread... + for(j=1;j<=(*num_channels);j++) { fscanf(subcat_fptr,"%lf %lf", @@ -508,16 +549,378 @@ if(yes_print_output==TRUE) return; } -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) + +/** + * Define functions to be used in init() and in Set_value in bmi_topmodel.c + * to enable calibratable parameters to be updated. + * BChoat 2023/08/29 + */ + + +/** + * Function to convert distance/area form to time delay histogram ordinates + * + * Converts parameters to m/time step DT + * Based on the calculation of values going into output tch, it seems TOPMODEL assumes + * one main channel with up to 9 areas of overland flow contributing to the + * one main channel. tch[1] is the main channel, and tch[>1] are the areas of + * overland flow (BChoat). + * + * There is one assumed velocity for the main channel and one assumed velocity applied + * to all overland flow areas. (BChoat) + * + * It is not clear why tch is hardcoded as an array of a fixed length. This is also done + * in the original fortran code on which Fred Ogden based this C-code. In Fortran + * it is coded as a 10-dimensional array and tch(10), and here as tch[11]. + * + * TCH(10) in Fortran gives you a ten element array indexed from 1..10. When this was ported, + * the 1-based indexing was directly ported as well, so the tch buffer was over allocated by one, + * giving an array indexed 0..10, but index 0 is ignored/unsued. (BChoat) + * + * NJF Refactored to dynamically allocate tch based on the number of channels provided + * + * @params[in] dist_from_outlet, pointer to array of length num_channels of type double, + * distance from outlet to point on channel with area known (i.e., to a channel; BChoat) + * @params[in] num_channels, integer, how many channels are present (m/hr) + * @params[in] CALIBRATABLE chv, pointer of type double and length one, the average channel velocity (m/hr) + * for the main channel (BChoat) + * @params[in] CALIBRATABLE rv, pointer of type double and length one, average internal overland flow routing velocity + * @params[in] dt, double of length one, timestep in hours + + * @params[out] tch, double* of length num_channels, holds histogram ordinates for each channel (BChoat) + * tch is used as input in subsequent functions + * Note, position 0 is ignored. It was + * written this way when translated from the original code. + */ + +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) +{ + //This function ASSUMES tch is overallocated by 1 and is indexable from (1, num_channels) + //validate invariants + if(num_channels < 1){ + printf("ERROR: convert_dist_to_histords, num channels < 1, must have at least one channel\n"); + exit(-1); //TODO return int error code and handle error externally + } + // declare local variables + double chvdt, rvdt; + int j; + + chvdt = *chv*dt; // distance water travels in one timestep within channel + rvdt = *rv*dt; //distance water travels as overland flow in one timestep + + tch[1]=dist_from_outlet[1]/chvdt; + for(j=2;j<=num_channels;j++) + { + + tch[j]=tch[1]+(dist_from_outlet[j]-dist_from_outlet[1])/rvdt; + } + + return; + +} + + + +/** + * Function to calculate the time delay histogram + * + * @params[in] num_channels, int, defined in subcat.dat file + * @params[in] area, double between 0 and 1 defining catchment area as ratio of + * entire catchment + * @params[in] tch, double* of length num_channels, holds histogram ordinates for each channel + * output from conver_dist_to_histords() + * Note, position 0 is ignored. It was + * written this way when translated from the original code. + * @params[in], cum_dist_area_with_dist, pointer of length num_channels-1 and type double, + * cumulative distribution of area with dist_from_outlet. + + * @params[out] num_time_delay_histo_ords, pointer of type int, number of time delay histogram ordinates + * @params[out] num_delay, pointer of type int, number of time steps lag (delay) in channel within + * catchment to outlet. + * @parms[out], time_delay_histogram, pointer of size and length max_time_delay_ordinates and type double, + * time lag of outflows due to channel routing + */ + +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) + +{ + // declare local variables + double time, sumar; //, a1, a2; + int j, ir; + + // casting tch[num_channels] to int truncates tch[num_channels] + // (e.g., 7.9 becomes 7) + //Determine how many ROUTING ORDINATES + //computed, essentially, by the `dist_from_outlet` of the FARTHEST channel segment + (*num_time_delay_histo_ords)=(int)tch[num_channels]; + + // this is here to round up. Since casting tch as int effectively rounds down, + // we add a value of 1 effectively rounding up. + if((double)(*num_time_delay_histo_ords) WARN_HISTOGRAM_ORDINATES){ + printf("WARNING: number of time delay hisogram ordinates, %d, is greater than %d\n", + *num_time_delay_histo_ords, + WARN_HISTOGRAM_ORDINATES); + } + + if((*time_delay_histogram) == NULL){ + d_alloc(time_delay_histogram, *num_time_delay_histo_ords); + } //FIXME always free/realloc + //If not, the caller must ensure the correct size of time_delay_histogram prior to calling + + //NJF so we build histogram with ordinates between 1 and "distance" between first and last channel + for(ir=1;ir<=(*num_time_delay_histo_ords);ir++) + { + time=(double)(*num_delay)+(double)ir; + if(time>tch[num_channels]) + { + (*time_delay_histogram)[ir]=1.0; + } + else + { + for(j=2;j<=num_channels;j++) + { + if(time<=tch[j]) + { + (*time_delay_histogram)[ir]= + cum_dist_area_with_dist[j-1]+ + (cum_dist_area_with_dist[j]-cum_dist_area_with_dist[j-1])* + (time-tch[j-1])/(tch[j]-tch[j-1]); + break; /* exits this for loop */ + } + } + } + } + + sumar = 0; + //Convert ordinates to cummulative area, should sum to 1 + for(int i = *num_time_delay_histo_ords; i > 1; i--){ + (*time_delay_histogram)[i] -= (*time_delay_histogram)[i-1]*area; + (*time_delay_histogram)[i] *= area; + sumar += (*time_delay_histogram)[i]; + } + sumar += (*time_delay_histogram)[1]; + if(sumar < 0.99999 || sumar > 1.00001){ + printf("ERROR: Histogram oridnates do not sum to 1.\n"); + printf("Check that the correct number of values for cum_dist_area_with_dist and dist_from_outlet are provided.\n"); + printf("The number of values for each variable should be equal to the number of channels (i.e., num_channels).\n\n"); + exit(-1); //FIXME this fuction should probably return an error code + //and the error be handled elsewhere, not just an exit here... + } + + return; +} + + + +/** + * Function to (re)initialize discharge array + * @params[in] stand_alone, int, 0 for running with ngen and 1 for running in stand alone mode + * @params[in] num_delay, pointer of type int, number of time steps lag (delay) in channel within + * catchment to outlet. Output from calc_time_delay_histogram + * @params[in] Q0, pointer of type double, initial subsurface flow per unit area + * @params[in] area, pointer of type double between 0 and 1 defining catchment area as ratio of + * entire catchment + * @params[in] num_time_delay_histo_ords, pointer of type int, number of time delay histogram ordinates, + * output from calc_time_delay_histogram + * @parms[in], time_delay_histogram, double pointer of type double, time lag of outflows due to channel routing, + * output from calc_time_delay_histogram + * + * @params[out] Q, pointer of type double and length num_delay+num_time_delay_histo_ords, simulated discharge + */ +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) +{ + // if not in stand alone mode, then reallocate Q + if(stand_alone != TRUE) { + // if Q is already allocated, need to free and re-compute the required size + if(*Q != NULL){ + free(*Q); + *Q = NULL; + } + //*Q = calloc(*num_delay + *num_time_delay_histo_ords + 1, sizeof(double)); + d_alloc(Q, *num_delay + *num_time_delay_histo_ords); + } + + // declare local variables + double sum; + int i, in; + + sum=0.0; + + for(i=1;i<=(*num_delay);i++) + { + (*Q)[i]+=(*Q0)*area; + } + + for(i=1;i<=(*num_time_delay_histo_ords);i++) + { + sum+=(*time_delay_histogram)[i]; + in=(*num_delay)+i; + (*Q)[in]+=(*Q0)*(area-sum); + }; + + return; +} + + +/** + * Function to initialize unsaturated zone storage and deficit + * + * @params[in] num_topodex_values, int, number of topodex histogram values + * (i.e., number of A/TANB ordinates), defines the size of stor_unsat_zone, deficit_root_zone, and deficit_local + * @params[in] dt, double, timestep in hours + * @params[in] CALIBRATABLE sr0, pointer of type double, initial root zone storage deficit + * @params[in] CALIBRATABLE szm, pointer of type double, exponential scaling parameter for decline + * of transmissivity with increase in storage deficit + * @params[in] Q0, pointer of type double, initial subsurface flow per unit area + * @params[in] CALIBRATABLE t0, pointer of type double, downslope transmissivity when soil is just saturated + * to the surface. 'input as LN(t0) + * @params[in] tl, double, not well defined, but related to time lag + + * @params[out] stor_unsat_zone, double pointer of type double of size num_topodex_values, + * storage in the unsaturated zone + * @params[out] szq, pointer of type double, not well defined, but related to rate at which moisture + * is lost from the subsruface. + * @params[out] deficit_local, double pointer of type double of size num_topodex_values, local storage (or saturation) + * deficit + * @params[out] deficit_root_zone, double pointer of type double of size num_topodex_values, root zone storage deficit + * @params[out] sbar, pointer of type double, catchment average soil mositure deficit + * edited dynamically as model steps through time. + * @params[out] bal, pointer of type double, residual of water balance + */ + +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) +{ + // declare local variables + double t0dt; + int ia; + + if((*stor_unsat_zone) == NULL){ + d_alloc(stor_unsat_zone, num_topodex_values); + } + if((*deficit_root_zone) == NULL){ + d_alloc(deficit_root_zone, num_topodex_values); + } + if((*deficit_local) == NULL){ + d_alloc(deficit_local, num_topodex_values); + } + //FIXME realloc if any of these were NULL so they are guaranteed the correct size? Or carefully + //document the assumption and the requirement for caller to size these arrays to num_topodex_values + + t0dt=(*t0)+log(dt); /* was ALOG - specific log function in fortran*/ + + /* Calculate SZQ parameter */ + (*szq)=exp(t0dt-tl); + + for(ia=1;ia<=num_topodex_values;ia++) + { + (*stor_unsat_zone)[ia]=0.0; + (*deficit_root_zone)[ia]=(*sr0); + } + + (*sbar)=-(*szm)*log((*Q0)/(*szq)); + + // Initialise water balance. BAL is positive for storage + (*bal)=-(*sbar)-(*sr0); + + return; +} + + +/** + * Main initialize function. Calls convert_dist_to_histords(), calc_time_delay_histogram(), + * init_water_balance(), and init_discharge_array() + * + * @params[in] in_param_fptr, FILE pointer, file with parameters (e.g., params.dat) + * @params[in] output_fptr, FILE pointer, file to which output will be written (e.g., topmod-cat.out) + * @params[in] subcat, pointer of type char, name of subcatchment, read in from in_param_fptr file + * @params[in] stand_alone, int, 0 for running with ngen and 1 for running in stand alone mode + * @params[in] num_channels, int, defined in subcat.dat file + * @params[in] num_topodex_values, int, number of topodex histogram values + * (i.e., number of A/TANB ordinates) + * @params[in] yes_print_output, int (boolean), 1 = TRUE; write output files, 0 = FALSE; do not write + * @params[in] area, double between 0 and 1 defining catchment area as ratio of + * entire catchment + * @parms[in], time_delay_histogram, double pointer of type double, time lag of outflows due to channel routing, + * output from calc_time_delay_histogram + * @params[in], cum_dist_area_with_dist, pointer of type double, cumulative distribution of area with + * dist_from_outlet + * @params[in] dt, double, timestep in hours + * @params[in] tl, double, not well defined, but related to time lag + * @params[in] dist_from_outlet, pointer to array of length num_channels of type double, + * distance from outlet to point on channel with area known (i.e., to a channel; BChoat) + * @params[out] num_time_delay_histo_ords, pointer of type int, number of time delay histogram ordinates, + * output from calc_time_delay_histogram + * @params[out] num_delay, pointer of type int, number of time steps lag (delay) in channel within + * catchment to outlet. Output from calc_time_delay_histogram + + * The following 12 variables are read in from the params.dat file using fscanf + * ----------------- + * @params[out] CALIBRATABLE szm, pointer of type double, exponential scaling parameter for decline + * of transmissivity with increase in storage deficit + * @params[out] CALIBRATABLE t0, pointer of type double, downslope transmissivity when soil is just saturated + * to the surface. 'input as LN(t0) + * @params[out] CALIBRATABLE td, pointer of type double, unsaturated zone time delay for recharge to the + * saturated zone per unit of deficit, provided in params.dat + * Only read in within init(), not used. + * @params[out] CALIBRATABLE chv, double of length one, the average channel velocity (m/hr) + * for main channel (BChoat) + * @params[out] CALIBRATABLE rv, double of length one, average internal overland flow routing velocity + * @params[out] CALIBRATABLE srmax, pointer of type double, maximum root zone storage deficit provided in params.dat + * Only read in within init(), not used. + * @params[out] Q0, pointer of type double, initial subsurface flow per unit area + * @params[out] CALIBRATABLE sr0, pointer of type double, initial root zone storage deficit + * @params[out] infex, pointer of type int (boolean), 1 = TRUE; call subroutine to do infiltration excess calcs, + * Not typically appropiate in catchments where TOPMODEL is applicable (i.e., shallow highly + * permeable soils). 0 = FALSE (default), + * Only read in within init(), not used. + * @params[out] xk0, pointer of type double, surface soil hydraulic conductivity. + * Only read in within init(), not used. + * @params[out] hf, pointer of type double, wetting from suction for G&A soln. + * Only read in within init(), not used. + * @params[out] dth, pointer of type double, Water content change across the wetting front + * Only read in within init(), not used. + * ----------------- + + * @params[out] stor_unsat_zone, double pointer of type double of size num_topodex_values, + * storage in the unsaturated zone + * @params[out] deficit_local, double pointer of type double of size num_topodex_values, local storage (or saturation) + * deficit + * @params[out] deficit_root_zone, double pointer of type double of size num_topodex_valuesnum_topodex_values, + * root zone storage deficit + * @params[out] szq, pointer of type double, not well defined, but related to rate at which moisture + * is lost from the subsurface. + * @params[out] Q, pointer of type double of length num_delay, simulated discharge + * @params[out] sbar, pointer of type double, catchment average soil mositure deficit + * edited dynamically as model steps through time. + * @params[out] bal, pointer of type double, residual of water balance + * + */ +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) { /*************************************************************** @@ -525,95 +928,68 @@ extern void init(FILE *in_param_fptr,FILE *output_fptr,char *subcat, READS PARAMETER DATA ******************************************/ -double rv; /* internal overland flow routing velocity */ -double chv; /* average channel flow velocity */ -double tch[11],rvdt,chvdt,t0dt,time,a1,sumar; -double sum,ar2,a2; -int i,j,ia,in,ir; +double tch[num_channels+1]; //+1 to maintain 1 based indexing used by other routines +double sumar; +int ir; + -if((*stor_unsat_zone)==NULL) - { - d_alloc(stor_unsat_zone,max_atb_increments); - d_alloc(deficit_root_zone,max_atb_increments); - d_alloc(deficit_local,max_atb_increments); - } -if((*time_delay_histogram)==NULL) - { - d_alloc(time_delay_histogram,max_time_delay_ordinates); - } /* read in run parameters */ fgets(subcat,256,in_param_fptr); +printf("subcat: %s\n", subcat); + fscanf(in_param_fptr,"%lf %lf %lf %lf %lf %lf %lf %lf %d %lf %lf %lf", - szm,t0,td,&chv,&rv,srmax,Q0,sr0,infex,xk0,hf,dth); + szm,t0,td,chv,rv,srmax,Q0,sr0,infex,xk0,hf,dth); -/* Convert parameters to m/time step DT -* with exception of XK0 which must stay in m/h -* Q0 is already in m/time step - T0 is input as Ln(To) */ -rvdt=rv*dt; -chvdt=chv*dt; -t0dt=(*t0)+log(dt); /* was ALOG */ -/* Calculate SZQ parameter */ -(*szq)=exp(t0dt-tl); +printf("\n\nCalibratable parameters from params*.dat:\n"); -/* CONVERT DISTANCE/AREA FORM TO TIME DELAY HISTOGRAM ORDINATES */ +printf("\nET and recharge:\n"); +printf("srmax = %f\n", *srmax); +printf("td = %f\n", *td); -tch[1]=dist_from_outlet[1]/chvdt; -for(j=2;j<=num_channels;j++) - { - tch[j]=tch[1]+(dist_from_outlet[j]-dist_from_outlet[1])/rvdt; - } -(*num_time_delay_histo_ords)=(int)tch[num_channels]; +printf("\nDischarge:\n"); +printf("chv = %f\n", *chv); +printf("rv = %f\n", *rv); +printf("\nWater balance:\n"); +printf("szm = %f\n", *szm); +printf("sr0 = %f\n", *sr0); +printf("t0 = %f\n\n", *t0); -if((double)(*num_time_delay_histo_ords)tch[num_channels]) - { - (*time_delay_histogram)[ir]=1.0; - } - else - { - for(j=2;j<=num_channels;j++) - { - if(time<=tch[j]) - { - (*time_delay_histogram)[ir]= - cum_dist_area_with_dist[j-1]+ - (cum_dist_area_with_dist[j]-cum_dist_area_with_dist[j-1])* - (time-tch[j-1])/(tch[j]-tch[j-1]); - break; /* exits this for loop */ - } - } - } - } -a1=(*time_delay_histogram)[1]; -sumar=(*time_delay_histogram)[1]; -(*time_delay_histogram)[1]*=area; -if((*num_time_delay_histo_ords)>1) - { - for(ir=2;ir<=(*num_time_delay_histo_ords);ir++) - { - a2=(*time_delay_histogram)[ir]; - (*time_delay_histogram)[ir]=a2-a1; - a1=a2; - sumar+=(*time_delay_histogram)[ir]; - (*time_delay_histogram)[ir]*=area; - } - } if(yes_print_output==TRUE) { + + // calculate sum of hist ords for printing to output file + sumar = 0; + for(int i = *num_time_delay_histo_ords; i > 1; i--){ + sumar += (*time_delay_histogram)[i]; + } + sumar += (*time_delay_histogram)[1]; + fprintf(output_fptr,"SZQ = %12.5lf\n",(*szq)); fprintf(output_fptr,"Subcatchment routing data:\n"); fprintf(output_fptr,"Maximum Routing Delay %12.5lf\n",tch[num_channels]); @@ -623,42 +999,12 @@ if(yes_print_output==TRUE) fprintf(output_fptr,"%12.5lf ",(*time_delay_histogram)[ir]); } fprintf(output_fptr,"\n"); - } - -/* INITIALISE deficit_root_zone AND Q0 VALUES HERE - * SR0 IS INITIAL ROOT ZONE STORAGE DEFICIT BELOW FIELD CAPACITY - * Q0 IS THE INITIAL DISCHARGE FOR THIS SUBCATCHMENT - * - * INITIALISE STORES */ -for(ia=1;ia<=num_topodex_values;ia++) - { - (*stor_unsat_zone)[ia]=0.0; - (*deficit_root_zone)[ia]=(*sr0); - } -(*sbar)=-(*szm)*log((*Q0)/(*szq)); - -/* Reinitialise discharge array */ -sum=0.0; -for(i=1;i<=(*num_delay);i++) - { - Q[i]+=(*Q0)*area; - } - -for(i=1;i<=(*num_time_delay_histo_ords);i++) - { - sum+=(*time_delay_histogram)[i]; - in=(*num_delay)+i; - Q[in]+=(*Q0)*(area-sum); - } - -/* Initialise water balance. BAL is positive for storage */ -(*bal)=-(*sbar)-(*sr0); -if(yes_print_output==TRUE) - { + fprintf(output_fptr,"Initial BAL %12.5f\n",(*bal)); fprintf(output_fptr,"Initial SBAR %12.5f\n",(*sbar)); fprintf(output_fptr,"Initial SR0 %12.5f\n",(*sr0)); } + return; } @@ -678,6 +1024,7 @@ extern void expinf(int irof, int it, int rint, double* df, double* cumf, // BMI Adaption: df and cumf are now passed as pointers + double sum,fc,func,cd,xke,e,tp,r2,f1; double dfunc,fx,add,fac,constant,f,f2,xkf,szf,dth; int i,j; diff --git a/test/Makefile b/test/Makefile index aded50b..c99a47e 100644 --- a/test/Makefile +++ b/test/Makefile @@ -1,10 +1,10 @@ # define the C compiler to use -CC = gcc +CC = gcc -g # define any compile-time flags # for newer compiler, the right hand side is not needed... # ... for instance, remove the flags if using gnu/10.1.0 on Cheyenne. -CFLAGS = +CFLAGS = -fsanitize=address INCLUDES = -I../include diff --git a/test/main_unit_test_bmi.c b/test/main_unit_test_bmi.c index 3f5c993..0ab1a6d 100755 --- a/test/main_unit_test_bmi.c +++ b/test/main_unit_test_bmi.c @@ -274,7 +274,7 @@ main(void){ { status = model->get_value_ptr(model, var_name, (void**)(&var)); if (status == BMI_FAILURE)return BMI_FAILURE; - printf(" get value ptr: %f\n",var); + printf(" get value ptr: %f\n",*var); } // Go ahead and test set_value_*() for last time step here if (n == test_nstep){ @@ -338,7 +338,7 @@ main(void){ { status = model->get_value_ptr(model, var_name, (void**)(&var)); if (status == BMI_FAILURE)return BMI_FAILURE; - printf(" get value ptr: %f\n",var); + printf(" get value ptr: %f\n",*var); } // Go ahead and test set_value_*() for last time step here if (n == test_nstep){ @@ -375,7 +375,13 @@ main(void){ } } } + for (i=0; ifinalize(model); if (status == BMI_FAILURE) return BMI_FAILURE; + free(model); printf("\n******************\nEND BMI UNIT TEST\n\n"); } return 0;