Skip to content
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

Ptl mbal testing #21

Merged
merged 8 commits into from
Apr 4, 2024
Merged

Ptl mbal testing #21

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
82d63c3
tested LASAM over approximately 40k parameter sets. Resolved mass bal…
Mar 25, 2024
2b0bb13
tested LASAM over approximately 40k parameter sets. Resolved mass bal…
Mar 25, 2024
48b0f56
updated unit test
Mar 25, 2024
032072d
updating with suggestions from PR review
Apr 3, 2024
afc7673
Updated lgar.cxx in response to a few comments
Apr 4, 2024
ed41fc5
Removed a few comments that were not necessary
Apr 4, 2024
5cc3d70
when theta updates, it is now enforced to be greater than theta_r and…
Apr 4, 2024
c2f0c79
increased clarity of comment detailing how lgar_theta_mass_balance wi…
Apr 4, 2024
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
6 changes: 3 additions & 3 deletions include/all.hxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -272,12 +272,12 @@ extern int lgar_read_vG_param_file(char const* vG_param_file_name, int num_soil_
struct soil_properties_ *soil_properties);

// creates a surficial front (new top most wetting front)
extern void lgar_create_surficial_front(double *ponded_depth_cm, double *volin, double dry_depth,
extern void lgar_create_surficial_front(int num_layers, double *ponded_depth_cm, double *volin, double dry_depth,
double theta1, int *soil_type, double *cum_layer_thickness_cm,
double *frozen_factor, struct wetting_front **head, struct soil_properties_ *soil_properties);

// computes the infiltration capacity, fp, of the soil
extern double lgar_insert_water(bool use_closed_form_G, int nint, double timestep_h, double *ponded_depth,
extern double lgar_insert_water(bool use_closed_form_G, int nint, double timestep_h, double AET_demand_cm, double *ponded_depth,
double *volin_this_timestep, double precip_timestep_cm, int wf_free_drainge_demand,
int num_layers, double ponded_depth_max_cm, int *soil_type, double *cum_layer_thickness_cm,
double *frozen_factor, struct wetting_front* head, struct soil_properties_ *soil_properties);
Expand Down Expand Up @@ -315,7 +315,7 @@ extern int wetting_front_free_drainage(struct wetting_front* head);

// computes updated theta (soil moisture content) after moving down a wetting front; called for each wetting front to ensure mass is conserved
extern double lgar_theta_mass_balance(int layer_num, int soil_num, double psi_cm, double new_mass,
double prior_mass, double *delta_theta, double *layer_thickness_cm,
double prior_mass, double *AET_demand_cm, double *delta_theta, double *layer_thickness_cm,
int *soil_type, struct soil_properties_ *soil_properties);

/********************************************************************/
Expand Down
7 changes: 4 additions & 3 deletions src/aet.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,8 +28,6 @@ extern double calc_aet(double PET_timestep_cm, double time_step_h, double wiltin
struct wetting_front *current;

double theta_wp;

double relative_moisture_at_which_PET_equals_AET = 0.75;

double Se,theta_e,theta_r;
double vg_a, vg_m, vg_n;
Expand All @@ -46,10 +44,13 @@ extern double calc_aet(double PET_timestep_cm, double time_step_h, double wiltin
vg_n = soil_properties[soil_num].vg_n;

// compute theta field capacity
double theta_fc = (theta_e - theta_r) * relative_moisture_at_which_PET_equals_AET + theta_r;
double head_at_which_PET_equals_AET_cm = 340.9;//*10/33; //340.9 is 0.33 atm, expressed in water depth, which is a good field capacity for most soils.
//Coarser soils like sand will have a field capacity of 0.1 atm or so.
double theta_fc = calc_theta_from_h(head_at_which_PET_equals_AET_cm, vg_a,vg_m, vg_n, theta_e, theta_r);

double wp_head_theta = calc_theta_from_h(wilting_point_psi_cm, vg_a,vg_m, vg_n, theta_e, theta_r);


theta_wp = (theta_fc - wp_head_theta)*1/2 + wp_head_theta; // theta_50 in python

Se = calc_Se_from_theta(theta_wp,theta_e,theta_r);
Expand Down
24 changes: 20 additions & 4 deletions src/bmi_lgar.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -238,20 +238,37 @@ Update()
double temp_pd = 0.0; // necessary to assign zero precip due to the creation of new wetting front; AET will still be taken out of the layers

// move the wetting fronts without adding any water; this is done to close the mass balance
// and also to merge / cross if necessary
lgar_move_wetting_fronts(subtimestep_h, &temp_pd, wf_free_drainage_demand, volend_subtimestep_cm,
num_layers, &AET_subtimestep_cm, state->lgar_bmi_params.cum_layer_thickness_cm,
state->lgar_bmi_params.layer_soil_type, state->lgar_bmi_params.frozen_factor,
&state->head, state->state_previous, state->soil_properties);

if (temp_pd != 0.0){ //if temp_pd != 0.0, that means that some water left the model through the lower model bdy
volrech_subtimestep_cm = temp_pd;
volrech_timestep_cm += volrech_subtimestep_cm;
temp_pd = 0.0;
}

// depth of the surficial front to be created
dry_depth = lgar_calc_dry_depth(use_closed_form_G, nint, subtimestep_h, &delta_theta, state->lgar_bmi_params.layer_soil_type,
state->lgar_bmi_params.cum_layer_thickness_cm, state->lgar_bmi_params.frozen_factor,
state->head, state->soil_properties);

if (verbosity.compare("high") == 0) {
printf("State before moving creating new WF...\n");
listPrint(state->head);
}

lgar_create_surficial_front(&ponded_depth_subtimestep_cm, &volin_subtimestep_cm, dry_depth, state->head->theta,
lgar_create_surficial_front(num_layers, &ponded_depth_subtimestep_cm, &volin_subtimestep_cm, dry_depth, state->head->theta,
state->lgar_bmi_params.layer_soil_type, state->lgar_bmi_params.cum_layer_thickness_cm,
state->lgar_bmi_params.frozen_factor, &state->head, state->soil_properties);

if (verbosity.compare("high") == 0) {
printf("State after moving creating new WF...\n");
listPrint(state->head);
}

state->state_previous = NULL;
state->state_previous = listCopy(state->head);

Expand All @@ -269,7 +286,7 @@ Update()

if (ponded_depth_subtimestep_cm > 0 && !create_surficial_front) {

volrunoff_subtimestep_cm = lgar_insert_water(use_closed_form_G, nint, subtimestep_h, &ponded_depth_subtimestep_cm,
volrunoff_subtimestep_cm = lgar_insert_water(use_closed_form_G, nint, subtimestep_h, AET_subtimestep_cm, &ponded_depth_subtimestep_cm,
&volin_subtimestep_cm, precip_subtimestep_cm_per_h,
wf_free_drainage_demand, num_layers,
ponded_depth_max_cm, state->lgar_bmi_params.layer_soil_type,
Expand Down Expand Up @@ -319,7 +336,6 @@ Update()

volin_subtimestep_cm = volin_subtimestep_cm_temp;
}

/*----------------------------------------------------------------------*/
// calculate derivative (dz/dt) for all wetting fronts
lgar_dzdt_calc(use_closed_form_G, nint, ponded_depth_subtimestep_cm, state->lgar_bmi_params.layer_soil_type,
Expand Down Expand Up @@ -360,7 +376,7 @@ Update()
listPrint(state->head);
}

bool unexpected_local_error = fabs(local_mb) > 1.0E-6 ? true : false;
bool unexpected_local_error = fabs(local_mb) > 1.0E-4 ? true : false;

if (verbosity.compare("high") == 0 || verbosity.compare("low") == 0 || unexpected_local_error) {
printf("\nLocal mass balance at this timestep... \n\
Expand Down
Loading
Loading