Skip to content

Fix aet, replaces #142 #143

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
ajkhattak merged 4 commits into from
Mar 13, 2025
Merged

Fix aet, replaces #142 #143

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
4ce073b
chore: small formatting refactor
hellkite500 Mar 4, 2025
a88e820
feat: add method to compute et_from_retention_depth
hellkite500 Mar 4, 2025
3ee91c9
fix: take AET from surface retention depth
hellkite500 Mar 4, 2025
9270e6e
fix: try to prevent segfault when computing AET from retention depth
hellkite500 Mar 4, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 5 additions & 0 deletions include/cfe.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -91,6 +91,7 @@ struct evapotranspiration_structure {
double potential_et_m_per_timestep;
double reduced_potential_et_m_per_timestep;
double actual_et_from_rain_m_per_timestep;
double actual_et_from_retention_depth_m_per_timestep;
double actual_et_from_soil_m_per_timestep;
double actual_et_m_per_timestep;
};
Expand Down Expand Up @@ -122,6 +123,7 @@ struct massbal
double volin ;
double volout ;
double volend ;
double vol_et_from_retention_depth;
};
typedef struct massbal massbal;

Expand Down Expand Up @@ -158,6 +160,9 @@ extern void Xinanjiang_partitioning_scheme(double water_input_depth_m, double fi

extern void et_from_rainfall(double *timestep_rainfall_input_m, struct evapotranspiration_structure *et_struct);

extern void et_from_retention_depth(struct nash_cascade_parameters *nash_surface_params,
struct evapotranspiration_structure *et_struct);

extern void et_from_soil(struct conceptual_reservoir *soil_res, struct evapotranspiration_structure *et_struct,
struct NWM_soil_parameters *soil_parms);

Expand Down
14 changes: 9 additions & 5 deletions src/bmi_cfe.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1121,6 +1121,8 @@ int read_init_config_cfe(const char* config_file, cfe_state_struct* model)
// Finally, free the original string memory
free(giuh_originates_string_val);

//NJF Be explicit that nash_storage in this case is NULL
model->nash_surface_params.nash_storage = NULL;
}
else if(model->surface_runoff_scheme == NASH_CASCADE) {
if (is_N_nash_surface_set == FALSE) {
Expand Down Expand Up @@ -3227,6 +3229,7 @@ extern void initialize_volume_trackers(cfe_state_struct* cfe_ptr) {
cfe_ptr->vol_struct.vol_soil_start = cfe_ptr->soil_reservoir.storage_m;
cfe_ptr->vol_struct.vol_et_from_soil = 0;
cfe_ptr->vol_struct.vol_et_from_rain = 0;
cfe_ptr->vol_struct.vol_et_from_retention_depth = 0;
cfe_ptr->vol_struct.vol_et_to_atm = 0;
}

Expand Down Expand Up @@ -3325,11 +3328,12 @@ extern void mass_balance_check(cfe_state_struct* cfe_ptr){
giuh_residual = cfe_ptr->vol_struct.vol_runoff - cfe_ptr->vol_struct.vol_out_surface - cfe_ptr->vol_struct.vol_end_surface -
cfe_ptr->vol_struct.vol_runon_infilt;
printf("********************* SURFACE MASS BALANCE *********************\n");
printf(" Volume into surface = %8.4lf m\n",cfe_ptr->vol_struct.vol_runoff);
printf(" Volume out surface = %8.4lf m\n",cfe_ptr->vol_struct.vol_out_surface);
printf(" Volume end surface = %8.4lf m\n",cfe_ptr->vol_struct.vol_end_surface);
printf(" Runon infiltration = %8.4lf m\n",cfe_ptr->vol_struct.vol_runon_infilt);
printf(" Surface residual = %6.4e m\n",giuh_residual); // should equal zero
printf(" Volume into surface = %8.4lf m\n",cfe_ptr->vol_struct.vol_runoff);
printf(" Volume out surface = %8.4lf m\n",cfe_ptr->vol_struct.vol_out_surface);
printf(" Volume end surface = %8.4lf m\n",cfe_ptr->vol_struct.vol_end_surface);
printf(" Runon infiltration = %8.4lf m\n",cfe_ptr->vol_struct.vol_runon_infilt);
printf(" Vol_et_from_retention_depth = %8.4lf m\n",cfe_ptr->vol_struct.vol_et_from_retention_depth);
printf(" Surface 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");

Expand Down
63 changes: 54 additions & 9 deletions src/cfe.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -73,21 +73,44 @@ extern void cfe(
evap_struct->reduced_potential_et_m_per_timestep = evap_struct->potential_et_m_per_s * time_step_size;

evap_struct->actual_et_from_rain_m_per_timestep = 0;
if(timestep_rainfall_input_m > 0) {et_from_rainfall(&timestep_rainfall_input_m,evap_struct);}
if (timestep_rainfall_input_m > 0) {
et_from_rainfall(&timestep_rainfall_input_m,evap_struct);
}

massbal_struct->vol_et_from_rain = massbal_struct->vol_et_from_rain + evap_struct->actual_et_from_rain_m_per_timestep;
massbal_struct->vol_et_to_atm = massbal_struct->vol_et_to_atm + evap_struct->actual_et_from_rain_m_per_timestep;
massbal_struct->volout=massbal_struct->volout+evap_struct->actual_et_from_rain_m_per_timestep;
massbal_struct->vol_et_to_atm = massbal_struct->vol_et_to_atm + evap_struct->actual_et_from_rain_m_per_timestep;
massbal_struct->volout = massbal_struct->volout + evap_struct->actual_et_from_rain_m_per_timestep;

// AET from surface retention depth
evap_struct->actual_et_from_retention_depth_m_per_timestep = 0;
// NJF Need some way to know if nash surface is used nash_storage
// is properly allocated, otherwise this can segfault!
if (nash_surface_params->nash_storage != NULL &&
nash_surface_params->nash_storage[0] > 0.0 &&
evap_struct->reduced_potential_et_m_per_timestep > 0.0) {
et_from_retention_depth(nash_surface_params, evap_struct);
}

// LKC: Change this. Now evaporation happens before runoff calculation. This was creating issues since it modifies storage_m and not storage_deficit
massbal_struct->vol_et_from_retention_depth += evap_struct->actual_et_from_retention_depth_m_per_timestep;
massbal_struct->vol_et_to_atm += evap_struct->actual_et_from_retention_depth_m_per_timestep;
massbal_struct->volout += evap_struct->actual_et_from_retention_depth_m_per_timestep;
massbal_struct->vol_out_surface += evap_struct->actual_et_from_retention_depth_m_per_timestep;

// LKC: Change this. Now evaporation happens before runoff calculation. This was creating issues since it modifies
// storage_m and not storage_deficit
evap_struct->actual_et_from_soil_m_per_timestep = 0;
if(soil_reservoir_struct->storage_m > NWM_soil_params_struct.wilting_point_m)
{et_from_soil(soil_reservoir_struct, evap_struct, &NWM_soil_params_struct);}
if(soil_reservoir_struct->storage_m > NWM_soil_params_struct.wilting_point_m) {
et_from_soil(soil_reservoir_struct, evap_struct, &NWM_soil_params_struct);
}

massbal_struct->vol_et_from_soil = massbal_struct->vol_et_from_soil + evap_struct->actual_et_from_soil_m_per_timestep;
massbal_struct->vol_et_to_atm = massbal_struct->vol_et_to_atm + evap_struct->actual_et_from_soil_m_per_timestep;
massbal_struct->volout=massbal_struct->volout+evap_struct->actual_et_from_soil_m_per_timestep;
massbal_struct->vol_et_to_atm = massbal_struct->vol_et_to_atm + evap_struct->actual_et_from_soil_m_per_timestep;
massbal_struct->volout = massbal_struct->volout + evap_struct->actual_et_from_soil_m_per_timestep;

evap_struct->actual_et_m_per_timestep = evap_struct->actual_et_from_rain_m_per_timestep +
evap_struct->actual_et_from_retention_depth_m_per_timestep +
evap_struct->actual_et_from_soil_m_per_timestep;

evap_struct->actual_et_m_per_timestep=evap_struct->actual_et_from_rain_m_per_timestep+evap_struct->actual_et_from_soil_m_per_timestep;
// LKC: This needs to be calcualted here after et_from_soil since soil_reservoir_struct->storage_m changes
soil_reservoir_storage_deficit_m=(NWM_soil_params_struct.smcmax*NWM_soil_params_struct.D-soil_reservoir_struct->storage_m);

Expand Down Expand Up @@ -577,6 +600,28 @@ void et_from_rainfall(double *timestep_rainfall_input_m, struct evapotranspirati
}
}

//##############################################################
//########### ET FROM SURFACE RETENTION DEPTH ##############
//##############################################################
void et_from_retention_depth(struct nash_cascade_parameters *nash_surface_params, struct evapotranspiration_structure *et_struct)
{
// NJF try not to use nash_storage if it hasn't been initialized...
if(nash_surface_params->nash_storage == NULL) return;

if (et_struct->reduced_potential_et_m_per_timestep >= nash_surface_params->nash_storage[0]) {
et_struct->actual_et_from_retention_depth_m_per_timestep = nash_surface_params->nash_storage[0];
nash_surface_params->nash_storage[0] = 0.0;
}
else {
et_struct->actual_et_from_retention_depth_m_per_timestep = et_struct->reduced_potential_et_m_per_timestep;
//et_struct->reduced_potential_et_m_per_timestep = 0.0;
nash_surface_params->nash_storage[0] -= et_struct->actual_et_from_retention_depth_m_per_timestep;
}

et_struct->reduced_potential_et_m_per_timestep -= et_struct->actual_et_from_retention_depth_m_per_timestep;

}

//##############################################################
//#################### ET FROM SOIL ########################
//##############################################################
Expand Down