Skip to content

Commit

Permalink
Removed useVegType array and replaced all usage with the proper fract…
Browse files Browse the repository at this point in the history 
…ion* value. Resolution to issue #110
  • Loading branch information
@BrendenBe1
BrendenBe1 committed Jan 9, 2018
1 parent 66d7783 commit 329aa69
Show file tree
Hide file tree
Showing 3 changed files with 124 additions and 80 deletions.
192 changes: 123 additions & 69 deletions SW_Output.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2316,6 +2316,10 @@ static void get_swa(void)
#ifdef STEPWAT
TimeInt p;
memset(SXW.sum_dSWA_repartitioned, 0, sizeof(SXW.sum_dSWA_repartitioned)); // need to reset sum_dSWA_repartitioned each year
RealD val_forb = SW_MISSING;
RealD val_tree = SW_MISSING;
RealD val_shrub = SW_MISSING;
RealD val_grass = SW_MISSING;
#endif

LyrIndex i;
Expand Down Expand Up @@ -2391,27 +2395,19 @@ static void get_swa(void)
curr_crit_val = SW_VegProd.critSoilWater[j]; // get critical value for current veg type
// go through each critical value to see which ones need to be set for each veg_type
for(k=0; k<4; k++){
if(SW_VegProd.useVegType[j][0] == 0){ // set anything that is at a value not set to 0
swa_master[j][k][i] = 0.;
}
else if(SW_VegProd.useVegType[k][0] == 0){
swa_master[j][k][i] = 0.;
if(k == j){
// dont need to check for its own critical value
}
else{
if(k == j){
// dont need to check for its own critical value
new_crit_val = SW_VegProd.critSoilWater[k];
if(curr_crit_val < new_crit_val){ // need to store this value since it has access to it
swa_master[j][k][i] = swa_master[k][k][i]; // itclp(veg_type, new_critical_value, layer, timeperiod)
}
else{
new_crit_val = SW_VegProd.critSoilWater[k];
if(curr_crit_val < new_crit_val){ // need to store this value since it has access to it
swa_master[j][k][i] = swa_master[k][k][i]; // itclp(veg_type, new_critical_value, layer, timeperiod)
}
if(curr_crit_val > new_crit_val){ // need to set this value to 0 since it does not have access to it
swa_master[j][k][i] = 0.; // itclp(veg_type, new_critical_value, layer, timeperiod)
}
if(curr_crit_val == new_crit_val){
// do nothing, dont want to swap values
}
if(curr_crit_val > new_crit_val){ // need to set this value to 0 since it does not have access to it
swa_master[j][k][i] = 0.; // itclp(veg_type, new_critical_value, layer, timeperiod)
}
if(curr_crit_val == new_crit_val){
// do nothing, dont want to swap values
}
}
}
Expand Down Expand Up @@ -2449,6 +2445,26 @@ static void get_swa(void)
break;
}

// first set each veg type to its crit value defined by inputs
// calculate the SWA values
// for all veg types checking to see if set to be used or not in the .in files
if(SW_VegProd.fractionTree != 0)
val_tree = fmax(0., val - SW_Site.lyr[i]->swcBulk_atSWPcrit_tree);
else
val_tree = 0.;
if(SW_VegProd.fractionShrub != 0)
val_shrub = fmax(0., val - SW_Site.lyr[i]->swcBulk_atSWPcrit_shrub);
else
val_shrub = 0.;
if(SW_VegProd.fractionGrass != 0)
val_grass = fmax(0., val - SW_Site.lyr[i]->swcBulk_atSWPcrit_grass);
else
val_grass = 0.;
if(SW_VegProd.fractionForb != 0)
val_forb = fmax(0., val - SW_Site.lyr[i]->swcBulk_atSWPcrit_forb);
else
val_forb = 0.;

// ########################################################################
// # Getting values for SWA_master at each critical value it has acess to #
// ########################################################################
Expand All @@ -2458,10 +2474,10 @@ static void get_swa(void)
*/
//yesterday_swa = SXW.SWA_master[Itclp(1,1,i,p)];
// first set each veg type to its crit value defined by inputs
SXW.SWA_master[Itclp(0,0,i,p)] = fmax(0., val - SW_Site.lyr[i]->swcBulk_atSWPcrit_tree);
SXW.SWA_master[Itclp(1,1,i,p)] = fmax(0., val - SW_Site.lyr[i]->swcBulk_atSWPcrit_shrub);
SXW.SWA_master[Itclp(2,2,i,p)] = fmax(0., val - SW_Site.lyr[i]->swcBulk_atSWPcrit_grass);
SXW.SWA_master[Itclp(3,3,i,p)] = fmax(0., val - SW_Site.lyr[i]->swcBulk_atSWPcrit_forb);
SXW.SWA_master[Itclp(0,0,i,p)] = val_tree;
SXW.SWA_master[Itclp(1,1,i,p)] = val_shrub;
SXW.SWA_master[Itclp(2,2,i,p)] = val_grass;
SXW.SWA_master[Itclp(3,3,i,p)] = val_forb;

// need to check which other critical value each veg_type has access to aside from its own
//(i.e. if shrub=-3.9 then it also has access to -3.5 and -2.0)
Expand All @@ -2473,27 +2489,19 @@ static void get_swa(void)
curr_crit_val = SW_VegProd.critSoilWater[j]; // get critical value for current veg type
// go through each critical value to see which ones need to be set for each veg_type
for(k=0; k<4; k++){
if(SW_VegProd.useVegType[j][0] == 0){ // set anything that is at a value not set to 0
SXW.SWA_master[Itclp(j,k,i,p)] = 0.;
}
else if(SW_VegProd.useVegType[k][0] == 0){
SXW.SWA_master[Itclp(j,k,i,p)] = 0.;
if(k == j){
// dont need to check for its own critical value
}
else{
if(k == j){
// dont need to check for its own critical value
new_crit_val = SW_VegProd.critSoilWater[k];
if(curr_crit_val < new_crit_val){ // need to store this value since it has access to it
SXW.SWA_master[Itclp(j,k,i,p)] = SXW.SWA_master[Itclp(k,k,i,p)]; // itclp(veg_type, new_critical_value, layer, timeperiod)
}
else{
new_crit_val = SW_VegProd.critSoilWater[k];
if(curr_crit_val < new_crit_val){ // need to store this value since it has access to it
SXW.SWA_master[Itclp(j,k,i,p)] = SXW.SWA_master[Itclp(k,k,i,p)]; // itclp(veg_type, new_critical_value, layer, timeperiod)
}
if(curr_crit_val > new_crit_val){ // need to set this value to 0 since it does not have access to it
SXW.SWA_master[Itclp(j,k,i,p)] = 0.; // itclp(veg_type, new_critical_value, layer, timeperiod)
}
if(curr_crit_val == new_crit_val){
// do nothing, dont want to swap values
}
if(curr_crit_val > new_crit_val){ // need to set this value to 0 since it does not have access to it
SXW.SWA_master[Itclp(j,k,i,p)] = 0.; // itclp(veg_type, new_critical_value, layer, timeperiod)
}
if(curr_crit_val == new_crit_val){
// do nothing, dont want to swap values
}
}
}
Expand Down Expand Up @@ -2573,14 +2581,27 @@ static void get_dSWAbulk(RealF swa_master[16][16][20], int i, int p){
RealD val = SW_MISSING;

#if !defined(STEPWAT) && !defined(RSOILWAT)
float veg_type_in_use; // set to current veg type fraction value to avoid multiple if loops. should just need 1 instead of 3 now.
float inner_loop_veg_type; // set to inner loop veg type
smallestCritVal = SW_VegProd.critSoilWater[SW_VegProd.rank_SWPcrits[0]];
largestCritVal = SW_VegProd.critSoilWater[SW_VegProd.rank_SWPcrits[3]];
RealF dSWA_bulk[16][16][20];
RealF dSWA_bulk_repartioned[16][16][20];
RealF dSWA_repartitioned_sum[16][16];

// loop through each veg type to get dSWAbulk
for(curr_vegType = 3; curr_vegType >= 0; curr_vegType--){ // go through each veg type and recalculate if necessary. starts at smallest
curr_crit_rank_index = SW_VegProd.rank_SWPcrits[curr_vegType]; // get rank index for start of next loop
// set veg type fraction here
if(curr_crit_rank_index == 0)
veg_type_in_use = SW_VegProd.fractionTree;
else if(curr_crit_rank_index == 1)
veg_type_in_use = SW_VegProd.fractionShrub;
else if(curr_crit_rank_index == 2)
veg_type_in_use = SW_VegProd.fractionGrass;
else
veg_type_in_use = SW_VegProd.fractionForb;

for(kv=curr_vegType; kv>=0; kv--){
crit_val = SW_VegProd.critSoilWater[SW_VegProd.rank_SWPcrits[kv]]; // get crit value at current index
kv_veg_type = SW_VegProd.rank_SWPcrits[kv]; // get index for veg_type. dont want to access swa_master at rank_SWPcrits index
Expand All @@ -2592,7 +2613,7 @@ static void get_dSWAbulk(RealF swa_master[16][16][20], int i, int p){
prev_crit_val = SW_VegProd.critSoilWater[SW_VegProd.rank_SWPcrits[kv]]; // get crit value for index lower
prev_crit_veg_type = SW_VegProd.rank_SWPcrits[kv]; // get veg type that belongs to the corresponding critical value
}
if(SW_VegProd.useVegType[curr_crit_rank_index][0] == 0){ // [0=tree(-2.0,off), 1=shrub(-3.9,on), 2=grass(-3.5,on), 3=forb(-2.0,on)]
if(veg_type_in_use == 0){ // [0=tree(-2.0,off), 1=shrub(-3.9,on), 2=grass(-3.5,on), 3=forb(-2.0,on)]
dSWA_bulk[curr_crit_rank_index][kv_veg_type][i] = 0.; // set to 0 to ensure no absent values
swa_master[curr_crit_rank_index][kv_veg_type][i] = 0.;
dSWA_bulk_repartioned[curr_crit_rank_index][kv_veg_type][i] = 0.;
Expand Down Expand Up @@ -2622,18 +2643,28 @@ static void get_dSWAbulk(RealF swa_master[16][16][20], int i, int p){
else{ // all values other than largest well need repartitioning
if(crit_val == smallestCritVal){ // if smallest value then all veg_types have access to it so just need to multiply by its fraction
dSWA_bulk_repartioned[curr_crit_rank_index][kv_veg_type][i] =
dSWA_bulk[curr_crit_rank_index][kv_veg_type][i] * SW_VegProd.useVegType[curr_crit_rank_index][1];
dSWA_bulk[curr_crit_rank_index][kv_veg_type][i] * veg_type_in_use;
// multiply by fraction for index of curr_vegType not kv
}
else{ // critical values that more than one veg type have access to but less than all veg types
vegFractionSum = 0; // set to 0 each time to make sure it gets correct value
// need to calculute new fraction for these since sum of them no longer adds up to 1. do this by adding fractions values
// of veg types who have access and then dividing these fraction values by the total sum. keeps the ratio and adds up to 1
for(j=3;j>=0;j--){ // go through all critical values and sum the fractions of the veg types who have access
// set veg type fraction here
if(j == 0)
inner_loop_veg_type = SW_VegProd.fractionTree;
else if(j == 1)
inner_loop_veg_type = SW_VegProd.fractionShrub;
else if(j == 2)
inner_loop_veg_type = SW_VegProd.fractionGrass;
else
inner_loop_veg_type = SW_VegProd.fractionForb;

if(SW_VegProd.critSoilWater[j] <= crit_val)
vegFractionSum += SW_VegProd.useVegType[j][1];
vegFractionSum += inner_loop_veg_type;
}
newFraction = SW_VegProd.useVegType[curr_crit_rank_index][1] / vegFractionSum; // divide veg fraction by sum to get new fraction value
newFraction = veg_type_in_use / vegFractionSum; // divide veg fraction by sum to get new fraction value
dSWA_bulk_repartioned[curr_crit_rank_index][kv_veg_type][i] =
dSWA_bulk[curr_crit_rank_index][kv_veg_type][i] * newFraction;
}
Expand Down Expand Up @@ -2667,6 +2698,8 @@ static void get_dSWAbulk(RealF swa_master[16][16][20], int i, int p){
printf("---------------------\n\n");*/

#elif defined(STEPWAT)
float veg_type_in_use; // set to current veg type fraction value to avoid multiple if loops. should just need 1 instead of 3 now.
float inner_loop_veg_type; // set to inner loop veg type
smallestCritVal = SW_VegProd.critSoilWater[SXW.rank_SWPcrits[0]];
largestCritVal = SW_VegProd.critSoilWater[SXW.rank_SWPcrits[3]];

Expand Down Expand Up @@ -2708,6 +2741,17 @@ static void get_dSWAbulk(RealF swa_master[16][16][20], int i, int p){
// loop through each veg type to get dSWAbulk
for(curr_vegType = 3; curr_vegType >= 0; curr_vegType--){ // go through each veg type and recalculate if necessary. starts at smallest
curr_crit_rank_index = SXW.rank_SWPcrits[curr_vegType]; // get rank index for start of next loop

// set veg type fraction here
if(curr_crit_rank_index == 0)
veg_type_in_use = SW_VegProd.fractionTree;
else if(curr_crit_rank_index == 1)
veg_type_in_use = SW_VegProd.fractionShrub;
else if(curr_crit_rank_index == 2)
veg_type_in_use = SW_VegProd.fractionGrass;
else
veg_type_in_use = SW_VegProd.fractionForb;

for(kv=curr_vegType; kv>=0; kv--){
crit_val = SW_VegProd.critSoilWater[SXW.rank_SWPcrits[kv]]; // get crit value at current index
kv_veg_type = SXW.rank_SWPcrits[kv]; // get index for veg_type. dont want to access swa_master at rank_SWPcrits index
Expand All @@ -2722,30 +2766,30 @@ static void get_dSWAbulk(RealF swa_master[16][16][20], int i, int p){

//printf("%f,%f\n\n", crit_val, prev_crit_val);

if(SW_VegProd.useVegType[curr_crit_rank_index][0] == 0){ // [0=tree(-2.0,off), 1=shrub(-3.9,on), 2=grass(-3.5,on), 3=forb(-2.0,on)]
// do nothing since veg type is turned off
//printf("%d,%d\n", curr_crit_rank_index, kv_veg_type);
SXW.SWA_master[Itclp(curr_crit_rank_index,kv_veg_type,i,p)] = 0.;
SXW.dSWAbulk[Itclp(curr_crit_rank_index,kv_veg_type,i,p)] = 0.;
SXW.dSWA_repartitioned[Itclp(curr_crit_rank_index,kv_veg_type,i,p)] = 0.;
}
else{ // check if need to recalculate for veg types in use
if(crit_val < prev_crit_val){ // if true then we need to recalculate
if(SXW.SWA_master[Itclp(curr_crit_rank_index,kv_veg_type,i,p)] == 0){
SXW.dSWAbulk[Itclp(curr_crit_rank_index,kv_veg_type,i,p)] = 0.;
if(veg_type_in_use == 0){ // [0=tree(-2.0,off), 1=shrub(-3.9,on), 2=grass(-3.5,on), 3=forb(-2.0,on)]
// do nothing since veg type is turned off
//printf("%d,%d\n", curr_crit_rank_index, kv_veg_type);
SXW.SWA_master[Itclp(curr_crit_rank_index,kv_veg_type,i,p)] = 0.;
SXW.dSWAbulk[Itclp(curr_crit_rank_index,kv_veg_type,i,p)] = 0.;
SXW.dSWA_repartitioned[Itclp(curr_crit_rank_index,kv_veg_type,i,p)] = 0.;
}
else{ // check if need to recalculate for veg types in use
if(crit_val < prev_crit_val){ // if true then we need to recalculate
if(SXW.SWA_master[Itclp(curr_crit_rank_index,kv_veg_type,i,p)] == 0){
SXW.dSWAbulk[Itclp(curr_crit_rank_index,kv_veg_type,i,p)] = 0.;
}
else{
SXW.dSWAbulk[Itclp(curr_crit_rank_index,kv_veg_type,i,p)] =
SXW.SWA_master[Itclp(curr_crit_rank_index,kv_veg_type,i,p)]-SXW.SWA_master[Itclp(curr_crit_rank_index,prev_crit_veg_type,i,p)];
}
}
else if(crit_val == prev_crit_val){ // critical values equal just set to itself
SXW.dSWAbulk[Itclp(curr_crit_rank_index,kv_veg_type,i,p)] = SXW.SWA_master[Itclp(curr_crit_rank_index,kv_veg_type,i,p)];
//printf("%f,%f\n\n", crit_val, prev_crit_val);
}
else{
SXW.dSWAbulk[Itclp(curr_crit_rank_index,kv_veg_type,i,p)] =
SXW.SWA_master[Itclp(curr_crit_rank_index,kv_veg_type,i,p)]-SXW.SWA_master[Itclp(curr_crit_rank_index,prev_crit_veg_type,i,p)];
}
}
else if(crit_val == prev_crit_val){ // critical values equal just set to itself
SXW.dSWAbulk[Itclp(curr_crit_rank_index,kv_veg_type,i,p)] = SXW.SWA_master[Itclp(curr_crit_rank_index,kv_veg_type,i,p)];
//printf("%f,%f\n\n", crit_val, prev_crit_val);
}
else{
// do nothing if crit val >. this will be handled later
}
// do nothing if crit val >. this will be handled later
}
/* ##############################################################
below if else blocks are for redistributing dSWAbulk values
############################################################## */
Expand All @@ -2754,18 +2798,28 @@ static void get_dSWAbulk(RealF swa_master[16][16][20], int i, int p){
else{ // all values other than largest well need repartitioning
if(crit_val == smallestCritVal){ // if smallest value then all veg_types have access to it so just need to multiply by its fraction
SXW.dSWA_repartitioned[Itclp(curr_crit_rank_index,kv_veg_type,i,p)] =
SXW.dSWAbulk[Itclp(curr_crit_rank_index,kv_veg_type,i,p)] * SW_VegProd.useVegType[curr_crit_rank_index][1];
SXW.dSWAbulk[Itclp(curr_crit_rank_index,kv_veg_type,i,p)] * veg_type_in_use;
// multiply by fraction for index of curr_vegType not kv
}
else{ // critical values that more than one veg type have access to but less than all veg types
vegFractionSum = 0; // set to 0 each time to make sure it gets correct value
// need to calculute new fraction for these since sum of them no longer adds up to 1. do this by adding fractions values
// of veg types who have access and then dividing these fraction values by the total sum. keeps the ratio and adds up to 1
for(j=3;j>=0;j--){ // go through all critical values and sum the fractions of the veg types who have access
// set veg type fraction here
if(j == 0)
inner_loop_veg_type = SW_VegProd.fractionTree;
else if(j == 1)
inner_loop_veg_type = SW_VegProd.fractionShrub;
else if(j == 2)
inner_loop_veg_type = SW_VegProd.fractionGrass;
else
inner_loop_veg_type = SW_VegProd.fractionForb;

if(SW_VegProd.critSoilWater[j] <= crit_val)
vegFractionSum += SW_VegProd.useVegType[j][1];
vegFractionSum += inner_loop_veg_type;
}
newFraction = SW_VegProd.useVegType[curr_crit_rank_index][1] / vegFractionSum; // divide veg fraction by sum to get new fraction value
newFraction = veg_type_in_use / vegFractionSum; // divide veg fraction by sum to get new fraction value
SXW.dSWA_repartitioned[Itclp(curr_crit_rank_index,kv_veg_type,i,p)] =
SXW.dSWAbulk[Itclp(curr_crit_rank_index,kv_veg_type,i,p)] * newFraction; // get repartioned soilwater using new fraction value
}
Expand Down
5 changes: 0 additions & 5 deletions SW_VegProd.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -117,11 +117,6 @@ void SW_VPD_read(void) {
v->fractionForb = help_forb;
v->fractionBareGround = help_bareGround;

SW_VegProd.useVegType[0][0] = help_tree;
SW_VegProd.useVegType[1][0] = help_shrub;
SW_VegProd.useVegType[2][0] = help_grass;
SW_VegProd.useVegType[3][0] = help_forb;

break;

/* albedo */
Expand Down
7 changes: 1 addition & 6 deletions SW_VegProd.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -96,12 +96,7 @@ typedef struct {
fractionBareGround; /* bare ground component fraction of total vegetation */

RealD critSoilWater[4]; // storing values in same order as defined in rgroup.in (0=tree, 1=shrub, 2=grass, 3=forb)
RealD useVegType[4][2]; // storing which veg types are set to be used and the fraction applied to them. [*][0] = .in file val. [*][1] = updated fractions
// the first element for each type is the fixed input and the second element is the updated fraction value.
// So for example SW_VegProd.useVegType[1][0] would always be -3.9 while SW_VegProd.useVegType[1][1] would be the updated value calculated in update_productivity.
// The reason for doing it this way was so all the values would be stored in the same array so we could just use the proper index in the for loop rather
// than doing an if statement to see if we need to use v->fractionShrub or v->fractionGrass.


int rank_SWPcrits[5]; // array to store the SWP crits in order of lest negative to most negative (used in sxw_resource)

RealD bareGround_albedo; /* create this here instead of creating a bareGround VegType, because it only needs albedo and no other data member */
Expand Down

0 comments on commit 329aa69

Please sign in to comment.