-
Notifications
You must be signed in to change notification settings - Fork 19
/
soil_water_stress.c
117 lines (94 loc) · 3.45 KB
/
soil_water_stress.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
/// @file soil_water_stress.c
/// @brief Compute soil water stress factor.
/// @note Please refer to file soil_water_factor.cpp for the original version. LHE.
/// @version 2.0
/// @authors Rewritten by: Liming He. Jan 29, 2013.
/// @authors Modified by: Mustapha El Maayar. March 2008
/// @authors Written by: Weimin Ju
/// @authors Last revision by: Liming He.
/// @date May 22, 2015.
#include <stdio.h>
#include <math.h>
#include "debug.h"
#include "soil.h"
/// @brief Function to compute soil water stress factor
/// @details [output] dt, fw-soil water stress
/// @param p soil conditions struct
/// @return void
void soil_water_factor_v2(struct Soil p[])
{
double ft[MAX_LAYERS], fpsisr[MAX_LAYERS];
double dtt[MAX_LAYERS];
double t1, t2;
// double fw;
int i;
t1 = -0.02;
t2 = 2.0;
if (p->psim[0] <= 0.000001) // just in case that this function is called before "updatesoilmoisture". LHE
for (i = 0; i < p->n_layer; i++)
{
p->psim[i] = p->psi_sat[i] * pow(p->thetam[i] / p->fei[i], -p->b[i]);
p->psim[i] = max(p->psi_sat[i], p->psim[i]); // I see no necessity to use this line unless thetam > fei. LHE May 20, 2015
}
for (i = 0; i < p->n_layer; i++)
{
if (p->psim[i] > p->psi_min) // changed 10.0 to psi_min. LHE. Feb. 13, 2013.
fpsisr[i] = 1.0 / (1 + pow((p->psim[i] - p->psi_min) / p->psi_min, p->alpha)); //psi_sr in m H2O! This is the old version. LHE.
else
fpsisr[i] = 1.0;
if (p->temp_soil_p[i] > 0.0) // using a previous temperature | soil water factor is calculated before Tm is updated. LHE.
ft[i] = (1.0 - exp(t1 * pow(p->temp_soil_p[i], t2))); // 1/x
else
ft[i] = 0.0;
// ft[i] = min(ft[i], 1.0); // LHE.
//ft[i] = max(0.001, ft[i]); // changed 0.1 to 0.001. LHE.
fpsisr[i] = fpsisr[i] * ft[i];
}
if (FW_VERSION == 1)
{
for (i = 0; i < p->n_layer; i++)
dtt[i] = p->f_root[i] * fpsisr[i]; /* eq. 14 in Ju 2006 */
}
else
{
for (i = 0; i < p->n_layer; i++)
dtt[i] = p->f_root[i];
}
//dtt[0] = 0.0;
double dtt_sum = 0.0;
for (i = 0; i < p->n_layer; i++)
dtt_sum = dtt_sum + dtt[i];/* eq. 14 in JU2006 */
if (dtt_sum < 0.000001)
p->f_soilwater = 0.1; // when soil temperatures in all layers are <=0. LHe
else
{
for (i = 0; i < p->n_layer; i++)
{
p->dt[i] = dtt[i] / dtt_sum;
//if (p->dt[0] < 0.0000001)
//printf("%f\n", p->dt[0]);
if (isnan(p->dt[i]))
printf("%f\n", p->dt[0]);
}
//fpsisr[0] = 0;
double fpsisr_sum = 0;
for (i = 0; i < p->n_layer; i++)
{
fpsisr_sum = fpsisr_sum + fpsisr[i] * p->dt[i]; /* eq. 12, in Chen 2012 GBC; eq 15 in JU */
}
p->f_soilwater = max(0.1, fpsisr_sum);
}
// in-process value check
//printf("fpsisr[0]=%f\n",fpsisr[0]);
}
// a very simple fw by LHE July 29, 2014 @home
//double tan_theta = (1.5 - 0.3) / (soil->Layer[2].theta_vfc - soil->Layer[2].theta_vwp);
//fw = 0.3 + (soil->Layer[2].thetam - soil->Layer[2].theta_vwp) * tan_theta;
//fw = max(0.3, fw);
//double tan_theta = (1.0 - 0.3) / (soil->Layer[2].theta_vfc - soil->Layer[2].theta_vwp);
//fw = 0.3 + (soil->Layer[2].thetam - soil->Layer[2].theta_vwp) * tan_theta;
//fw = max(0.3, fw);
//
//fw = min(3.0, fw);
//
//soil->f_soilwater = fw;