Skip to content
GitLab

Source

Target

Hide whitespace changes
Inline Side-by-side
Showing
with 1991 additions and 1735 deletions
+1991 -1735
R/DataAltiExtrapolation_Valery.R
View file @ 4bccfee2
...@@ -3,543 +3,171 @@ DataAltiExtrapolation_Valery <- function(DatesR, ...@@ -3,543 +3,171 @@ DataAltiExtrapolation_Valery <- function(DatesR,
TempMean, TempMin = NULL, TempMax = NULL, TempMean, TempMin = NULL, TempMax = NULL,
ZInputs, HypsoData, NLayers, ZInputs, HypsoData, NLayers,
verbose = TRUE) { verbose = TRUE) {
##Altitudinal_gradient_functions_______________________________________________________________ ##Altitudinal_gradient_functions_______________________________________________________________
##unique_gradient_for_precipitation ##unique_gradient_for_precipitation
GradP_Valery2010 <- function() { GradP_Valery2010 <- 0.00041 ### value from Valery PhD thesis page 126
return(0.00041) ### value from Valerie PhD thesis page 126
}
##daily_gradients_for_mean_min_and_max_air_temperature
GradT_Valery2010 <- function() {
RESULT <- matrix(c(
01, 01, 0.434, 0.366, 0.498,
02, 01, 0.434, 0.366, 0.500,
03, 01, 0.435, 0.367, 0.501,
04, 01, 0.436, 0.367, 0.503,
05, 01, 0.437, 0.367, 0.504,
06, 01, 0.439, 0.367, 0.506,
07, 01, 0.440, 0.367, 0.508,
08, 01, 0.441, 0.368, 0.510,
09, 01, 0.442, 0.368, 0.512,
10, 01, 0.444, 0.368, 0.514,
11, 01, 0.445, 0.368, 0.517,
12, 01, 0.446, 0.368, 0.519,
13, 01, 0.448, 0.369, 0.522,
14, 01, 0.450, 0.369, 0.525,
15, 01, 0.451, 0.369, 0.527,
16, 01, 0.453, 0.370, 0.530,
17, 01, 0.455, 0.370, 0.533,
18, 01, 0.456, 0.370, 0.537,
19, 01, 0.458, 0.371, 0.540,
20, 01, 0.460, 0.371, 0.543,
21, 01, 0.462, 0.371, 0.547,
22, 01, 0.464, 0.372, 0.550,
23, 01, 0.466, 0.372, 0.554,
24, 01, 0.468, 0.373, 0.558,
25, 01, 0.470, 0.373, 0.561,
26, 01, 0.472, 0.374, 0.565,
27, 01, 0.474, 0.374, 0.569,
28, 01, 0.476, 0.375, 0.573,
29, 01, 0.478, 0.375, 0.577,
30, 01, 0.480, 0.376, 0.582,
31, 01, 0.483, 0.376, 0.586,
01, 02, 0.485, 0.377, 0.590,
02, 02, 0.487, 0.377, 0.594,
03, 02, 0.489, 0.378, 0.599,
04, 02, 0.492, 0.379, 0.603,
05, 02, 0.494, 0.379, 0.607,
06, 02, 0.496, 0.380, 0.612,
07, 02, 0.498, 0.381, 0.616,
08, 02, 0.501, 0.381, 0.621,
09, 02, 0.503, 0.382, 0.625,
10, 02, 0.505, 0.383, 0.630,
11, 02, 0.508, 0.384, 0.634,
12, 02, 0.510, 0.384, 0.639,
13, 02, 0.512, 0.385, 0.643,
14, 02, 0.515, 0.386, 0.648,
15, 02, 0.517, 0.387, 0.652,
16, 02, 0.519, 0.387, 0.657,
17, 02, 0.522, 0.388, 0.661,
18, 02, 0.524, 0.389, 0.666,
19, 02, 0.526, 0.390, 0.670,
20, 02, 0.528, 0.391, 0.674,
21, 02, 0.530, 0.392, 0.679,
22, 02, 0.533, 0.393, 0.683,
23, 02, 0.535, 0.393, 0.687,
24, 02, 0.537, 0.394, 0.691,
25, 02, 0.539, 0.395, 0.695,
26, 02, 0.541, 0.396, 0.699,
27, 02, 0.543, 0.397, 0.703,
28, 02, 0.545, 0.398, 0.707,
29, 02, 0.546, 0.399, 0.709,
01, 03, 0.547, 0.399, 0.711,
02, 03, 0.549, 0.400, 0.715,
03, 03, 0.551, 0.401, 0.718,
04, 03, 0.553, 0.402, 0.722,
05, 03, 0.555, 0.403, 0.726,
06, 03, 0.557, 0.404, 0.729,
07, 03, 0.559, 0.405, 0.732,
08, 03, 0.560, 0.406, 0.736,
09, 03, 0.562, 0.406, 0.739,
10, 03, 0.564, 0.407, 0.742,
11, 03, 0.566, 0.408, 0.745,
12, 03, 0.567, 0.409, 0.748,
13, 03, 0.569, 0.410, 0.750,
14, 03, 0.570, 0.411, 0.753,
15, 03, 0.572, 0.412, 0.756,
16, 03, 0.573, 0.413, 0.758,
17, 03, 0.575, 0.414, 0.761,
18, 03, 0.576, 0.415, 0.763,
19, 03, 0.577, 0.416, 0.765,
20, 03, 0.579, 0.417, 0.767,
21, 03, 0.580, 0.417, 0.769,
22, 03, 0.581, 0.418, 0.771,
23, 03, 0.582, 0.419, 0.773,
24, 03, 0.583, 0.420, 0.774,
25, 03, 0.584, 0.421, 0.776,
26, 03, 0.585, 0.422, 0.777,
27, 03, 0.586, 0.422, 0.779,
28, 03, 0.587, 0.423, 0.780,
29, 03, 0.588, 0.424, 0.781,
30, 03, 0.589, 0.425, 0.782,
31, 03, 0.590, 0.425, 0.783,
01, 04, 0.591, 0.426, 0.784,
02, 04, 0.591, 0.427, 0.785,
03, 04, 0.592, 0.427, 0.785,
04, 04, 0.593, 0.428, 0.786,
05, 04, 0.593, 0.429, 0.787,
06, 04, 0.594, 0.429, 0.787,
07, 04, 0.595, 0.430, 0.787,
08, 04, 0.595, 0.431, 0.788,
09, 04, 0.596, 0.431, 0.788,
10, 04, 0.596, 0.432, 0.788,
11, 04, 0.597, 0.432, 0.788,
12, 04, 0.597, 0.433, 0.788,
13, 04, 0.597, 0.433, 0.788,
14, 04, 0.598, 0.434, 0.788,
15, 04, 0.598, 0.434, 0.788,
16, 04, 0.598, 0.435, 0.787,
17, 04, 0.599, 0.435, 0.787,
18, 04, 0.599, 0.436, 0.787,
19, 04, 0.599, 0.436, 0.786,
20, 04, 0.599, 0.436, 0.786,
21, 04, 0.600, 0.437, 0.785,
22, 04, 0.600, 0.437, 0.785,
23, 04, 0.600, 0.437, 0.784,
24, 04, 0.600, 0.438, 0.784,
25, 04, 0.600, 0.438, 0.783,
26, 04, 0.601, 0.438, 0.783,
27, 04, 0.601, 0.438, 0.782,
28, 04, 0.601, 0.439, 0.781,
29, 04, 0.601, 0.439, 0.781,
30, 04, 0.601, 0.439, 0.780,
01, 05, 0.601, 0.439, 0.779,
02, 05, 0.601, 0.439, 0.778,
03, 05, 0.601, 0.439, 0.778,
04, 05, 0.601, 0.440, 0.777,
05, 05, 0.601, 0.440, 0.776,
06, 05, 0.601, 0.440, 0.775,
07, 05, 0.601, 0.440, 0.775,
08, 05, 0.601, 0.440, 0.774,
09, 05, 0.601, 0.440, 0.773,
10, 05, 0.602, 0.440, 0.772,
11, 05, 0.602, 0.440, 0.772,
12, 05, 0.602, 0.440, 0.771,
13, 05, 0.602, 0.440, 0.770,
14, 05, 0.602, 0.440, 0.770,
15, 05, 0.602, 0.440, 0.769,
16, 05, 0.602, 0.440, 0.768,
17, 05, 0.602, 0.440, 0.768,
18, 05, 0.602, 0.440, 0.767,
19, 05, 0.602, 0.440, 0.767,
20, 05, 0.602, 0.440, 0.766,
21, 05, 0.602, 0.440, 0.766,
22, 05, 0.602, 0.440, 0.765,
23, 05, 0.602, 0.440, 0.765,
24, 05, 0.602, 0.440, 0.764,
25, 05, 0.602, 0.440, 0.764,
26, 05, 0.602, 0.440, 0.764,
27, 05, 0.602, 0.439, 0.763,
28, 05, 0.602, 0.439, 0.763,
29, 05, 0.602, 0.439, 0.763,
30, 05, 0.602, 0.439, 0.762,
31, 05, 0.602, 0.439, 0.762,
01, 06, 0.602, 0.439, 0.762,
02, 06, 0.602, 0.439, 0.762,
03, 06, 0.602, 0.439, 0.762,
04, 06, 0.602, 0.439, 0.762,
05, 06, 0.602, 0.439, 0.762,
06, 06, 0.602, 0.438, 0.761,
07, 06, 0.602, 0.438, 0.761,
08, 06, 0.602, 0.438, 0.761,
09, 06, 0.602, 0.438, 0.761,
10, 06, 0.602, 0.438, 0.761,
11, 06, 0.602, 0.438, 0.762,
12, 06, 0.602, 0.438, 0.762,
13, 06, 0.602, 0.438, 0.762,
14, 06, 0.602, 0.438, 0.762,
15, 06, 0.602, 0.437, 0.762,
16, 06, 0.602, 0.437, 0.762,
17, 06, 0.602, 0.437, 0.762,
18, 06, 0.602, 0.437, 0.762,
19, 06, 0.602, 0.437, 0.763,
20, 06, 0.602, 0.437, 0.763,
21, 06, 0.602, 0.437, 0.763,
22, 06, 0.602, 0.436, 0.763,
23, 06, 0.602, 0.436, 0.763,
24, 06, 0.602, 0.436, 0.764,
25, 06, 0.602, 0.436, 0.764,
26, 06, 0.601, 0.436, 0.764,
27, 06, 0.601, 0.436, 0.764,
28, 06, 0.601, 0.436, 0.764,
29, 06, 0.601, 0.435, 0.765,
30, 06, 0.601, 0.435, 0.765,
01, 07, 0.601, 0.435, 0.765,
02, 07, 0.600, 0.435, 0.765,
03, 07, 0.600, 0.435, 0.765,
04, 07, 0.600, 0.434, 0.766,
05, 07, 0.600, 0.434, 0.766,
06, 07, 0.599, 0.434, 0.766,
07, 07, 0.599, 0.434, 0.766,
08, 07, 0.599, 0.434, 0.766,
09, 07, 0.598, 0.433, 0.766,
10, 07, 0.598, 0.433, 0.766,
11, 07, 0.598, 0.433, 0.766,
12, 07, 0.597, 0.433, 0.766,
13, 07, 0.597, 0.432, 0.767,
14, 07, 0.597, 0.432, 0.767,
15, 07, 0.596, 0.432, 0.767,
16, 07, 0.596, 0.432, 0.766,
17, 07, 0.595, 0.431, 0.766,
18, 07, 0.595, 0.431, 0.766,
19, 07, 0.594, 0.431, 0.766,
20, 07, 0.594, 0.430, 0.766,
21, 07, 0.593, 0.430, 0.766,
22, 07, 0.593, 0.430, 0.766,
23, 07, 0.592, 0.429, 0.765,
24, 07, 0.592, 0.429, 0.765,
25, 07, 0.591, 0.428, 0.765,
26, 07, 0.590, 0.428, 0.765,
27, 07, 0.590, 0.428, 0.764,
28, 07, 0.589, 0.427, 0.764,
29, 07, 0.588, 0.427, 0.764,
30, 07, 0.588, 0.426, 0.763,
31, 07, 0.587, 0.426, 0.763,
01, 08, 0.586, 0.425, 0.762,
02, 08, 0.586, 0.425, 0.762,
03, 08, 0.585, 0.424, 0.761,
04, 08, 0.584, 0.424, 0.761,
05, 08, 0.583, 0.423, 0.760,
06, 08, 0.583, 0.423, 0.760,
07, 08, 0.582, 0.422, 0.759,
08, 08, 0.581, 0.421, 0.758,
09, 08, 0.580, 0.421, 0.758,
10, 08, 0.579, 0.420, 0.757,
11, 08, 0.578, 0.420, 0.756,
12, 08, 0.578, 0.419, 0.755,
13, 08, 0.577, 0.418, 0.754,
14, 08, 0.576, 0.418, 0.754,
15, 08, 0.575, 0.417, 0.753,
16, 08, 0.574, 0.416, 0.752,
17, 08, 0.573, 0.415, 0.751,
18, 08, 0.572, 0.415, 0.750,
19, 08, 0.571, 0.414, 0.749,
20, 08, 0.570, 0.413, 0.748,
21, 08, 0.569, 0.413, 0.747,
22, 08, 0.569, 0.412, 0.746,
23, 08, 0.568, 0.411, 0.745,
24, 08, 0.567, 0.410, 0.744,
25, 08, 0.566, 0.409, 0.743,
26, 08, 0.565, 0.409, 0.742,
27, 08, 0.564, 0.408, 0.741,
28, 08, 0.563, 0.407, 0.740,
29, 08, 0.562, 0.406, 0.738,
30, 08, 0.561, 0.405, 0.737,
31, 08, 0.560, 0.405, 0.736,
01, 09, 0.558, 0.404, 0.735,
02, 09, 0.557, 0.403, 0.734,
03, 09, 0.556, 0.402, 0.732,
04, 09, 0.555, 0.401, 0.731,
05, 09, 0.554, 0.401, 0.730,
06, 09, 0.553, 0.400, 0.728,
07, 09, 0.552, 0.399, 0.727,
08, 09, 0.551, 0.398, 0.725,
09, 09, 0.550, 0.397, 0.724,
10, 09, 0.549, 0.396, 0.723,
11, 09, 0.548, 0.396, 0.721,
12, 09, 0.546, 0.395, 0.720,
13, 09, 0.545, 0.394, 0.718,
14, 09, 0.544, 0.393, 0.717,
15, 09, 0.543, 0.392, 0.715,
16, 09, 0.542, 0.391, 0.713,
17, 09, 0.541, 0.391, 0.712,
18, 09, 0.540, 0.390, 0.710,
19, 09, 0.538, 0.389, 0.709,
20, 09, 0.537, 0.388, 0.707,
21, 09, 0.536, 0.388, 0.705,
22, 09, 0.535, 0.387, 0.703,
23, 09, 0.533, 0.386, 0.702,
24, 09, 0.532, 0.385, 0.700,
25, 09, 0.531, 0.385, 0.698,
26, 09, 0.530, 0.384, 0.696,
27, 09, 0.528, 0.383, 0.694,
28, 09, 0.527, 0.383, 0.692,
29, 09, 0.526, 0.382, 0.690,
30, 09, 0.525, 0.381, 0.688,
01, 10, 0.523, 0.381, 0.686,
02, 10, 0.522, 0.380, 0.684,
03, 10, 0.521, 0.379, 0.682,
04, 10, 0.519, 0.379, 0.680,
05, 10, 0.518, 0.378, 0.678,
06, 10, 0.517, 0.377, 0.676,
07, 10, 0.515, 0.377, 0.674,
08, 10, 0.514, 0.376, 0.671,
09, 10, 0.512, 0.376, 0.669,
10, 10, 0.511, 0.375, 0.667,
11, 10, 0.510, 0.375, 0.664,
12, 10, 0.508, 0.374, 0.662,
13, 10, 0.507, 0.374, 0.659,
14, 10, 0.505, 0.373, 0.657,
15, 10, 0.504, 0.373, 0.654,
16, 10, 0.502, 0.372, 0.652,
17, 10, 0.501, 0.372, 0.649,
18, 10, 0.499, 0.372, 0.647,
19, 10, 0.498, 0.371, 0.644,
20, 10, 0.496, 0.371, 0.641,
21, 10, 0.495, 0.371, 0.639,
22, 10, 0.493, 0.370, 0.636,
23, 10, 0.492, 0.370, 0.633,
24, 10, 0.490, 0.370, 0.630,
25, 10, 0.489, 0.369, 0.628,
26, 10, 0.487, 0.369, 0.625,
27, 10, 0.485, 0.369, 0.622,
28, 10, 0.484, 0.368, 0.619,
29, 10, 0.482, 0.368, 0.616,
30, 10, 0.481, 0.368, 0.613,
31, 10, 0.479, 0.368, 0.610,
01, 11, 0.478, 0.368, 0.607,
02, 11, 0.476, 0.367, 0.604,
03, 11, 0.475, 0.367, 0.601,
04, 11, 0.473, 0.367, 0.598,
05, 11, 0.471, 0.367, 0.595,
06, 11, 0.470, 0.367, 0.592,
07, 11, 0.468, 0.367, 0.589,
08, 11, 0.467, 0.366, 0.586,
09, 11, 0.465, 0.366, 0.583,
10, 11, 0.464, 0.366, 0.580,
11, 11, 0.462, 0.366, 0.577,
12, 11, 0.461, 0.366, 0.574,
13, 11, 0.459, 0.366, 0.571,
14, 11, 0.458, 0.366, 0.568,
15, 11, 0.456, 0.366, 0.565,
16, 11, 0.455, 0.366, 0.562,
17, 11, 0.454, 0.366, 0.559,
18, 11, 0.452, 0.365, 0.556,
19, 11, 0.451, 0.365, 0.553,
20, 11, 0.450, 0.365, 0.550,
21, 11, 0.448, 0.365, 0.547,
22, 11, 0.447, 0.365, 0.544,
23, 11, 0.446, 0.365, 0.542,
24, 11, 0.445, 0.365, 0.539,
25, 11, 0.443, 0.365, 0.536,
26, 11, 0.442, 0.365, 0.533,
27, 11, 0.441, 0.365, 0.531,
28, 11, 0.440, 0.365, 0.528,
29, 11, 0.439, 0.365, 0.526,
30, 11, 0.438, 0.365, 0.523,
01, 12, 0.437, 0.365, 0.521,
02, 12, 0.436, 0.365, 0.519,
03, 12, 0.435, 0.365, 0.517,
04, 12, 0.434, 0.365, 0.515,
05, 12, 0.434, 0.365, 0.513,
06, 12, 0.433, 0.365, 0.511,
07, 12, 0.432, 0.365, 0.509,
08, 12, 0.431, 0.365, 0.507,
09, 12, 0.431, 0.365, 0.505,
10, 12, 0.430, 0.365, 0.504,
11, 12, 0.430, 0.365, 0.502,
12, 12, 0.429, 0.365, 0.501,
13, 12, 0.429, 0.365, 0.500,
14, 12, 0.429, 0.365, 0.498,
15, 12, 0.428, 0.365, 0.497,
16, 12, 0.428, 0.365, 0.496,
17, 12, 0.428, 0.365, 0.496,
18, 12, 0.428, 0.365, 0.495,
19, 12, 0.428, 0.365, 0.494,
20, 12, 0.428, 0.365, 0.494,
21, 12, 0.428, 0.365, 0.494,
22, 12, 0.428, 0.365, 0.493,
23, 12, 0.429, 0.365, 0.493,
24, 12, 0.429, 0.366, 0.493,
25, 12, 0.429, 0.366, 0.493,
26, 12, 0.430, 0.366, 0.494,
27, 12, 0.430, 0.366, 0.494,
28, 12, 0.431, 0.366, 0.495,
29, 12, 0.431, 0.366, 0.495,
30, 12, 0.432, 0.366, 0.496,
31, 12, 0.433, 0.366, 0.497), ncol = 5, byrow = TRUE)
dimnames(RESULT) <- list(1:366, c("day", "month", "grad_Tmean", "grad_Tmin", "grad_Tmax"))
return(RESULT)
}
##Format_______________________________________________________________________________________ ##Format_______________________________________________________________________________________
HypsoData <- as.double(HypsoData) HypsoData <- as.double(HypsoData)
ZInputs <- as.double(ZInputs) ZInputs <- as.double(ZInputs)
##ElevationLayers_Creation_____________________________________________________________________ ##ElevationLayers_Creation_____________________________________________________________________
ZLayers <- as.double(rep(NA, NLayers)) ZLayers <- as.double(rep(NA, NLayers))
if (!identical(HypsoData, as.double(rep(NA, 101)))) { if (!identical(HypsoData, as.double(rep(NA, 101)))) {
nmoy <- 100 %/% NLayers nmoy <- 100 %/% NLayers
nreste <- 100 %% NLayers nreste <- 100 %% NLayers
ncont <- 0 ncont <- 0
for (iLayer in 1:NLayers) { for (iLayer in 1:NLayers) {
if (nreste > 0) { if (nreste > 0) {
nn <- nmoy + 1 nn <- nmoy + 1
nreste <- nreste - 1 nreste <- nreste - 1
} else { } else {
nn <- nmoy nn <- nmoy
} }
if (nn == 1) { if (nn == 1) {
ZLayers[iLayer] <- HypsoData[ncont + 1] ZLayers[iLayer] <- HypsoData[ncont + 1]
}
if (nn == 2) {
ZLayers[iLayer] <- 0.5 * (HypsoData[ncont + 1] + HypsoData[ncont + 2])
}
if (nn > 2) {
ZLayers[iLayer] <- HypsoData[ncont + nn / 2]
}
ncont <- ncont + nn
} }
if (nn == 2) {
ZLayers[iLayer] <- 0.5 * (HypsoData[ncont + 1] + HypsoData[ncont + 2])
}
if (nn > 2) {
ZLayers[iLayer] <- HypsoData[ncont + nn / 2 + 1]
}
ncont <- ncont + nn
} }
}
##Precipitation_extrapolation__________________________________________________________________ ##Precipitation_extrapolation__________________________________________________________________
##Initialisation ##Initialisation
if (identical(ZInputs, HypsoData[51]) & NLayers == 1) { if (identical(ZInputs, HypsoData[51]) & NLayers == 1) {
LayerPrecip <- list(as.double(Precip)) LayerPrecip <- list(as.double(Precip))
} else { } else {
##Elevation_gradients_for_daily_mean_precipitation ##Elevation_gradients_for_daily_mean_precipitation
GradP <- GradP_Valery2010() ### single value GradP <- GradP_Valery2010 ### single value
TabGradP <- rep(GradP, length(Precip)) TabGradP <- rep(GradP, length(Precip))
##Extrapolation ##Extrapolation
##Thresold_of_inputs_median_elevation ##Thresold_of_inputs_median_elevation
Zthreshold <- 4000 Zthreshold <- 4000
LayerPrecip_df <- sapply(1:NLayers, function(iLayer) { LayerPrecip_mat <- sapply(1:NLayers, function(iLayer) {
##If_layer_elevation_smaller_than_Zthreshold ##If_layer_elevation_smaller_than_Zthreshold
if (ZLayers[iLayer] <= Zthreshold) { if (ZLayers[iLayer] <= Zthreshold) {
prcp <- as.double(Precip * exp(TabGradP * (ZLayers[iLayer] - ZInputs))) prcp <- as.double(Precip * exp(TabGradP * (ZLayers[iLayer] - ZInputs)))
##If_layer_elevation_greater_than_Zthreshold ##If_layer_elevation_greater_than_Zthreshold
} else {
##If_inputs_median_elevation_smaller_than_Zthreshold
if (ZInputs <= Zthreshold) {
prcp <- as.double(Precip * exp(TabGradP * (Zthreshold - ZInputs)))
##If_inputs_median_elevation_greater_then_Zthreshold
} else { } else {
##If_inputs_median_elevation_smaller_than_Zthreshold prcp <- as.double(Precip)
if (ZInputs <= Zthreshold) {
prcp <- as.double(Precip * exp(TabGradP * (Zthreshold - ZInputs)))
##If_inputs_median_elevation_greater_then_Zthreshold
} else {
prcp <- as.double(Precip)
}
} }
return(prcp)
})
if (PrecipScale) {
LayerPrecip_mat <- LayerPrecip_df / rowMeans(LayerPrecip_df) * Precip
LayerPrecip_mat[is.nan(LayerPrecip_mat)] <- 0
} }
LayerPrecip <- as.list(as.data.frame(LayerPrecip_mat)) return(prcp)
})
if (PrecipScale) {
LayerPrecip_mat <- LayerPrecip_mat / rowMeans(LayerPrecip_mat) * Precip
LayerPrecip_mat[is.nan(LayerPrecip_mat)] <- 0
} }
LayerPrecip <- as.list(as.data.frame(LayerPrecip_mat))
}
##Temperature_extrapolation____________________________________________________________________ ##Temperature_extrapolation____________________________________________________________________
##Initialisation ##Initialisation
LayerTempMean <- LayerTempMean <- list()
list() LayerTempMin <- list()
LayerTempMin <- list() LayerTempMax <- list()
LayerTempMax <- list()
if (identical(ZInputs, HypsoData[51]) & NLayers == 1) {
if (identical(ZInputs, HypsoData[51]) & NLayers == 1) { LayerTempMean[[1]] <- as.double(TempMean)
LayerTempMean[[1]] <- as.double(TempMean)
if (!is.null(TempMin) & !is.null(TempMax)) {
if (!is.null(TempMin) & !is.null(TempMax)) { LayerTempMin[[1]] <- as.double(TempMin)
LayerTempMin[[1]] <- as.double(TempMin) LayerTempMax[[1]] <- as.double(TempMax)
LayerTempMax[[1]] <- as.double(TempMax)
}
} else {
##Elevation_gradients_for_daily_mean_min_and_max_temperature
GradT <- as.data.frame(GradT_Valery2010())
iday <- match(format(DatesR, format = "%d%m"),
sprintf("%02i%02i", GradT[, "day"], GradT[, "month"]))
TabGradT <-
GradT[iday, c("grad_Tmean", "grad_Tmin", "grad_Tmax")]
##Extrapolation
##On_each_elevation_layer...
for (iLayer in 1:NLayers) {
LayerTempMean[[iLayer]] <- as.double(TempMean + (ZInputs - ZLayers[iLayer]) * abs(TabGradT[, "grad_Tmean"]) / 100)
if (!is.null(TempMin) & !is.null(TempMax)) {
LayerTempMin[[iLayer]] <- as.double(TempMin + (ZInputs - ZLayers[iLayer]) * abs(TabGradT[, "grad_Tmin"]) / 100)
LayerTempMax[[iLayer]] <- as.double(TempMax + (ZInputs - ZLayers[iLayer]) * abs(TabGradT[, "grad_Tmax"]) / 100)
}
}
} }
} else {
##Elevation_gradients_for_daily_mean_min_and_max_temperature
GradT <- .GradT_Valery2010
##Solid_Fraction_for_each_elevation_layer______________________________________________________ iday <- match(format(DatesR, format = "%d%m"),
LayerFracSolidPrecip <- list() sprintf("%02i%02i", GradT[, "day"], GradT[, "month"]))
TabGradT <- GradT[iday, c("grad_Tmean", "grad_Tmin", "grad_Tmax")]
##Thresold_of_inputs_median_elevation ##Extrapolation
Zthreshold <- 1500
##On_each_elevation_layer... ##On_each_elevation_layer...
for (iLayer in 1:NLayers) { for (iLayer in 1:NLayers) {
Option <- "USACE" LayerTempMean[[iLayer]] <- as.double(TempMean + (ZInputs - ZLayers[iLayer]) * abs(TabGradT[, "grad_Tmean"]) / 100)
if (!is.null(TempMin) & !is.null(TempMax)) {
if (!is.na(ZInputs)) { LayerTempMin[[iLayer]] <- as.double(TempMin + (ZInputs - ZLayers[iLayer]) * abs(TabGradT[, "grad_Tmin"]) / 100)
if (ZInputs < Zthreshold & !is.null(TempMin) & !is.null(TempMax)) { LayerTempMax[[iLayer]] <- as.double(TempMax + (ZInputs - ZLayers[iLayer]) * abs(TabGradT[, "grad_Tmax"]) / 100)
Option <- "Hydrotel"
}
}
##Turcotte_formula_from_Hydrotel
if (Option == "Hydrotel") {
TempMin <- LayerTempMin[[iLayer]]
TempMax <- LayerTempMax[[iLayer]]
SolidFraction <- 1 - TempMax / (TempMax - TempMin)
SolidFraction[TempMin >= 0] <- 0
SolidFraction[TempMax <= 0] <- 1
}
##USACE_formula
if (Option == "USACE") {
USACE_Tmin <- -1.0
USACE_Tmax <- 3.0
TempMean <- LayerTempMean[[iLayer]]
SolidFraction <- 1 - (TempMean - USACE_Tmin) / (USACE_Tmax - USACE_Tmin)
SolidFraction[TempMean > USACE_Tmax] <- 0
SolidFraction[TempMean < USACE_Tmin] <- 1
} }
LayerFracSolidPrecip[[iLayer]] <- as.double(SolidFraction)
} }
}
##Solid_Fraction_for_each_elevation_layer______________________________________________________
LayerFracSolidPrecip <- list()
##Thresold_of_inputs_median_elevation
Zthreshold <- 1500
##Option
Option <- "USACE"
if (!is.na(ZInputs)) {
if (ZInputs < Zthreshold & !is.null(TempMin) & !is.null(TempMax)) {
Option <- "Hydrotel"
}
}
##On_each_elevation_layer...
for (iLayer in 1:NLayers) {
##Turcotte_formula_from_Hydrotel
if (Option == "Hydrotel") {
TempMin <- LayerTempMin[[iLayer]]
TempMax <- LayerTempMax[[iLayer]]
SolidFraction <- 1 - TempMax / (TempMax - TempMin)
SolidFraction[TempMin >= 0] <- 0
SolidFraction[TempMax <= 0] <- 1
}
##USACE_formula
if (Option == "USACE") {
USACE_Tmin <- -1.0
USACE_Tmax <- 3.0
TempMean <- LayerTempMean[[iLayer]]
SolidFraction <- 1 - (TempMean - USACE_Tmin) / (USACE_Tmax - USACE_Tmin)
SolidFraction[TempMean > USACE_Tmax] <- 0
SolidFraction[TempMean < USACE_Tmin] <- 1
}
LayerFracSolidPrecip[[iLayer]] <- as.double(SolidFraction)
}
namesLayer <- sprintf("L%i", seq_along(LayerPrecip))
names(LayerPrecip) <- namesLayer
names(LayerTempMean) <- namesLayer
if (!is.null(TempMin) & !is.null(TempMax)) {
names(LayerTempMin) <- namesLayer
names(LayerTempMax) <- namesLayer
}
names(LayerFracSolidPrecip) <- namesLayer
##END__________________________________________________________________________________________ ##END__________________________________________________________________________________________
return(list(LayerPrecip = LayerPrecip, return(list(LayerPrecip = LayerPrecip,
LayerTempMean = LayerTempMean, LayerTempMean = LayerTempMean,
LayerTempMin = LayerTempMin, LayerTempMin = LayerTempMin,
LayerTempMax = LayerTempMax, LayerTempMax = LayerTempMax,
LayerFracSolidPrecip = LayerFracSolidPrecip, LayerFracSolidPrecip = LayerFracSolidPrecip,
ZLayers = ZLayers)) ZLayers = ZLayers))
}
}
R/ErrorCrit.R
View file @ 4bccfee2
ErrorCrit <- function(InputsCrit,OutputsModel,FUN_CRIT, warnings = TRUE, verbose = TRUE){ ErrorCrit <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE) {
return( FUN_CRIT(InputsCrit,OutputsModel, warnings = warnings, verbose = verbose) )
## ---------- Arguments check
if (!inherits(InputsCrit, "InputsCrit")) {
stop("InputsCrit must be of class 'InputsCrit'")
}
if (!inherits(OutputsModel, "OutputsModel")) {
stop("OutputsModel must be of class 'OutputsModel'")
}
## ---------- Criterion computation
## ----- Single criterion
if (inherits(InputsCrit, "Single")) {
FUN_CRIT <- match.fun(InputsCrit$FUN_CRIT)
OutputsCrit <- FUN_CRIT(InputsCrit = InputsCrit,
OutputsModel = OutputsModel,
warnings = warnings,
verbose = verbose)
}
## ----- Multiple criteria or Composite criterion
if (inherits(InputsCrit, "Multi") | inherits(InputsCrit, "Compo")) {
listOutputsCrit <- lapply(InputsCrit, FUN = function(iInputsCrit) {
FUN_CRIT <- match.fun(iInputsCrit$FUN_CRIT)
FUN_CRIT(InputsCrit = iInputsCrit,
OutputsModel = OutputsModel,
warnings = warnings,
verbose = verbose)
})
listValCrit <- sapply(listOutputsCrit, function(x) x[["CritValue"]])
listNameCrit <- sapply(listOutputsCrit, function(x) x[["CritName"]])
listweights <- unlist(lapply(InputsCrit, function(x) x[["Weights"]]))
listweights <- listweights / sum(listweights)
if (inherits(InputsCrit, "Compo")) {
CritValue <- sum(listValCrit * listweights)
OutputsCritCompo <- list(MultiCritValues = listValCrit,
MultiCritNames = listNameCrit,
MultiCritWeights = listweights)
OutputsCrit <- list(CritValue = CritValue,
CritName = "Composite",
CritBestValue = +1,
Multiplier = -1,
Ind_notcomputed = NULL,
CritCompo = OutputsCritCompo,
MultiCrit = listOutputsCrit)
class(OutputsCrit) <- c("Compo", "ErrorCrit")
if (verbose) {
message("------------------------------------\n")
message("Crit. Composite = ", sprintf("%.4f", CritValue))
msgForm <- paste(sprintf("%.2f", listweights), listNameCrit, sep = " * ", collapse = ", ")
msgForm <- unlist(strsplit(msgForm, split = ","))
msgFormSep <- rep(c(",", ",", ",\n\t\t "), times = ceiling(length(msgForm)/3))[1:length(msgForm)]
msgForm <- paste(msgForm, msgFormSep, sep = "", collapse = "")
msgForm <- gsub("\\,\\\n\\\t\\\t $|\\,$", "", msgForm)
message("\tFormula: sum(", msgForm, ")\n")
}
} else {
OutputsCrit <- listOutputsCrit
class(OutputsCrit) <- c("Multi", "ErrorCrit")
}
}
return(OutputsCrit)
} }
R/ErrorCrit_KGE.R
View file @ 4bccfee2
ErrorCrit_KGE <- function(InputsCrit,OutputsModel, warnings = TRUE, verbose = TRUE){ ErrorCrit_KGE <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE) {
## Arguments check
##Arguments_check________________________________ if (!inherits(OutputsModel, "OutputsModel")) {
if(inherits(InputsCrit,"InputsCrit")==FALSE){ stop("InputsCrit must be of class 'InputsCrit' \n"); return(NULL); } stop("'OutputsModel' must be of class 'OutputsModel'")
if(inherits(OutputsModel,"OutputsModel")==FALSE){ stop("OutputsModel must be of class 'OutputsModel' \n"); return(NULL); }
##Initialisation_________________________________
CritName <- NA;
if(InputsCrit$transfo=="" ){ CritName <- "KGE[Q]" ; }
if(InputsCrit$transfo=="sqrt"){ CritName <- "KGE[sqrt(Q)]"; }
if(InputsCrit$transfo=="log" ){ CritName <- "KGE[log(Q)]" ; }
if(InputsCrit$transfo=="inv" ){ CritName <- "KGE[1/Q]" ; }
if(InputsCrit$transfo=="sort"){ CritName <- "KGE[sort(Q)]"; }
CritValue <- NA;
CritBestValue <- +1;
Multiplier <- -1; ### must be equal to -1 or +1 only
##Data_preparation_______________________________
VarObs <- InputsCrit$Qobs ; VarObs[!InputsCrit$BoolCrit] <- NA;
VarSim <- OutputsModel$Qsim; VarSim[!InputsCrit$BoolCrit] <- NA;
##Data_transformation
if("Ind_zeroes" %in% names(InputsCrit) & "epsilon" %in% names(InputsCrit)){ if(length(InputsCrit$Ind_zeroes)>0){
VarObs <- VarObs + InputsCrit$epsilon;
VarSim <- VarSim + InputsCrit$epsilon;
} }
if(InputsCrit$transfo=="sqrt"){ VarObs <- sqrt(VarObs); VarSim <- sqrt(VarSim); }
if(InputsCrit$transfo=="log" ){ VarObs <- log(VarObs) ; VarSim <- log(VarSim) ; VarSim[VarSim < -1E100] <- NA; }
if(InputsCrit$transfo=="inv" ){ VarObs <- 1/VarObs ; VarSim <- 1/VarSim ; VarSim[abs(VarSim) > 1E+100] <- NA; }
if(InputsCrit$transfo=="sort"){
VarSim[is.na(VarObs)] <- NA
VarSim <- sort(VarSim, na.last = TRUE)
VarObs <- sort(VarObs, na.last = TRUE)
InputsCrit$BoolCrit <- sort(InputsCrit$BoolCrit, decreasing = TRUE)
} }
##TS_ignore
TS_ignore <- !is.finite(VarObs) | !is.finite(VarSim) | !InputsCrit$BoolCrit ; EC <- .ErrorCrit(InputsCrit = InputsCrit, Crit = "KGE", OutputsModel = OutputsModel, warnings = warnings)
Ind_TS_ignore <- which(TS_ignore); if(length(Ind_TS_ignore)==0){ Ind_TS_ignore <- NULL; }
if(sum(!TS_ignore)==0){ OutputsCrit <- list(NA); names(OutputsCrit) <- c("CritValue"); return(OutputsCrit); } CritValue <- NA
if(sum(!TS_ignore)==1){ OutputsCrit <- list(NA); names(OutputsCrit) <- c("CritValue"); return(OutputsCrit); } ### to avoid a problem in standard deviation computation SubCritValues <- rep(NA, 3)
if(inherits(OutputsModel,"hourly" )){ WarningTS <- 365; } SubCritNames <- c("r", "alpha", "beta")
if(inherits(OutputsModel,"daily" )){ WarningTS <- 365; } SubCritPrint <- rep(NA, 3)
if(inherits(OutputsModel,"monthly")){ WarningTS <- 12; }
if(inherits(OutputsModel,"yearly" )){ WarningTS <- 3; } if (EC$CritCompute) {
if(sum(!TS_ignore)<WarningTS & warnings){ warning("\t criterion computed on less than ", WarningTS, " time-steps ") } ## Other variables preparation
##Other_variables_preparation meanVarObs <- mean(EC$VarObs[!EC$TS_ignore])
meanVarObs <- mean(VarObs[!TS_ignore]); meanVarSim <- mean(EC$VarSim[!EC$TS_ignore])
meanVarSim <- mean(VarSim[!TS_ignore]);
iCrit <- 0 ## SubErrorCrit KGE rPearson
SubCritPrint <- NULL SubCritPrint[1L] <- paste0(EC$CritName, " cor(sim, obs, \"pearson\") =")
SubCritNames <- NULL
SubCritValues <- NULL Numer <- sum((EC$VarObs[!EC$TS_ignore] - meanVarObs) * (EC$VarSim[!EC$TS_ignore] - meanVarSim))
Deno1 <- sqrt(sum((EC$VarObs[!EC$TS_ignore] - meanVarObs) ^ 2))
Deno2 <- sqrt(sum((EC$VarSim[!EC$TS_ignore] - meanVarSim) ^ 2))
##SubErrorCrit_____KGE_rPearson__________________ if (Numer == 0) {
iCrit <- iCrit+1; if (Deno1 == 0 & Deno2 == 0) {
SubCritPrint[iCrit] <- paste(CritName," cor(sim, obs, \"pearson\") =", sep = "") Crit <- 1
SubCritValues[iCrit] <- NA; } else {
SubCritNames[iCrit] <- "r" Crit <- 0
Numer <- sum( (VarObs[!TS_ignore]-meanVarObs)*(VarSim[!TS_ignore]-meanVarSim) ); }
Deno1 <- sqrt( sum((VarObs[!TS_ignore]-meanVarObs)^2) ); } else {
Deno2 <- sqrt( sum((VarSim[!TS_ignore]-meanVarSim)^2) ); Crit <- Numer / (Deno1 * Deno2)
if(Numer==0){ if(Deno1==0 & Deno2==0){ Crit <- 1; } else { Crit <- 0; } }
} else { Crit <- Numer/(Deno1*Deno2); } if (is.numeric(Crit) & is.finite(Crit)) {
if(is.numeric(Crit) & is.finite(Crit)){ SubCritValues[iCrit] <- Crit; } SubCritValues[1L] <- Crit
}
##SubErrorCrit_____KGE_alpha_____________________ ## SubErrorCrit KGE alpha
iCrit <- iCrit+1; SubCritPrint[2L] <- paste0(EC$CritName, " sd(sim)/sd(obs) =")
SubCritPrint[iCrit] <- paste(CritName," sd(sim)/sd(obs) =", sep = "")
SubCritValues[iCrit] <- NA; Numer <- sd(EC$VarSim[!EC$TS_ignore])
SubCritNames[iCrit] <- "alpha" Denom <- sd(EC$VarObs[!EC$TS_ignore])
Numer <- sd(VarSim[!TS_ignore]);
Denom <- sd(VarObs[!TS_ignore]); if (Numer == 0 & Denom == 0) {
if(Numer==0 & Denom==0){ Crit <- 1; } else { Crit <- Numer/Denom ; } Crit <- 1
if(is.numeric(Crit) & is.finite(Crit)){ SubCritValues[iCrit] <- Crit; } } else {
Crit <- Numer / Denom
}
##SubErrorCrit_____KGE_beta______________________ if (is.numeric(Crit) & is.finite(Crit)) {
iCrit <- iCrit+1; SubCritValues[2L] <- Crit
SubCritPrint[iCrit] <- paste(CritName," mean(sim)/mean(obs) =", sep = "") }
SubCritValues[iCrit] <- NA;
SubCritNames[iCrit] <- "beta" ## SubErrorCrit KGE beta
if(meanVarSim==0 & meanVarObs==0){ Crit <- 1; } else { Crit <- meanVarSim/meanVarObs ; } SubCritPrint[3L] <- paste0(EC$CritName, " mean(sim)/mean(obs) =")
if(is.numeric(Crit) & is.finite(Crit)){ SubCritValues[iCrit] <- Crit; }
if (meanVarSim == 0 & meanVarObs == 0) {
Crit <- 1
##ErrorCrit______________________________________ } else {
if(sum(is.na(SubCritValues))==0){ Crit <- meanVarSim / meanVarObs
CritValue <- ( 1 - sqrt( (SubCritValues[1]-1)^2 + (SubCritValues[2]-1)^2 + (SubCritValues[3]-1)^2 ) ); }
if (is.numeric(Crit) & is.finite(Crit)) {
SubCritValues[3L] <- Crit
}
## ErrorCrit
if (sum(is.na(SubCritValues)) == 0) {
CritValue <- (1 - sqrt((SubCritValues[1L] - 1)^2 + (SubCritValues[2L] - 1)^2 + (SubCritValues[3L] - 1)^2))
}
## Verbose
if (verbose) {
message(sprintf("Crit. %s = %.4f", EC$CritName, CritValue))
message(paste("\tSubCrit.", SubCritPrint, sprintf("%.4f", SubCritValues), "\n", sep = " "))
}
} }
##Verbose______________________________________ ## Output
if(verbose) { OutputsCrit <- list(CritValue = CritValue,
message("Crit. ", CritName, " = ", sprintf("%.4f", CritValue)) CritName = EC$CritName,
message(paste("\tSubCrit.", SubCritPrint, sprintf("%.4f", SubCritValues), "\n", sep = " ")) SubCritValues = SubCritValues,
} SubCritNames = SubCritNames,
CritBestValue = EC$CritBestValue,
Multiplier = EC$Multiplier,
Ind_notcomputed = EC$Ind_TS_ignore)
##Output_________________________________________ class(OutputsCrit) <- c("KGE", "ErrorCrit")
OutputsCrit <- list(CritValue = CritValue, CritName = CritName,
SubCritValues = SubCritValues, SubCritNames = SubCritNames, CritBestValue = CritBestValue,
Multiplier = Multiplier, Ind_notcomputed = Ind_TS_ignore)
return(OutputsCrit) return(OutputsCrit)
} }
class(ErrorCrit_KGE) <- c("FUN_CRIT", class(ErrorCrit_KGE))
R/ErrorCrit_KGE2.R
View file @ 4bccfee2
ErrorCrit_KGE2 <- function(InputsCrit,OutputsModel, warnings = TRUE, verbose = TRUE){ ErrorCrit_KGE2 <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE) {
## Arguments check
##Arguments_check________________________________ if (!inherits(OutputsModel, "OutputsModel")) {
if(inherits(InputsCrit,"InputsCrit")==FALSE){ stop("InputsCrit must be of class 'InputsCrit' \n"); return(NULL); } stop("'OutputsModel' must be of class 'OutputsModel'")
if(inherits(OutputsModel,"OutputsModel")==FALSE){ stop("OutputsModel must be of class 'OutputsModel' \n"); return(NULL); }
##Initialisation_________________________________
CritName <- NA;
if(InputsCrit$transfo=="" ){ CritName <- "KGE'[Q]" ; }
if(InputsCrit$transfo=="sqrt"){ CritName <- "KGE'[sqrt(Q)]"; }
if(InputsCrit$transfo=="log" ){ CritName <- "KGE'[log(Q)]" ; }
if(InputsCrit$transfo=="inv" ){ CritName <- "KGE'[1/Q]" ; }
if(InputsCrit$transfo=="sort"){ CritName <- "KGE'[sort(Q)]"; }
CritValue <- NA;
CritBestValue <- +1;
Multiplier <- -1; ### must be equal to -1 or +1 only
##Data_preparation_______________________________
VarObs <- InputsCrit$Qobs ; VarObs[!InputsCrit$BoolCrit] <- NA;
VarSim <- OutputsModel$Qsim; VarSim[!InputsCrit$BoolCrit] <- NA;
##Data_transformation
if("Ind_zeroes" %in% names(InputsCrit) & "epsilon" %in% names(InputsCrit)){ if(length(InputsCrit$Ind_zeroes)>0){
VarObs <- VarObs + InputsCrit$epsilon;
VarSim <- VarSim + InputsCrit$epsilon;
} }
if(InputsCrit$transfo=="sqrt"){ VarObs <- sqrt(VarObs); VarSim <- sqrt(VarSim); }
if(InputsCrit$transfo=="log" ){ VarObs <- log(VarObs) ; VarSim <- log(VarSim) ; VarSim[VarSim < -1E100] <- NA; }
if(InputsCrit$transfo=="inv" ){ VarObs <- 1/VarObs ; VarSim <- 1/VarSim ; VarSim[abs(VarSim) > 1E+100] <- NA; }
if(InputsCrit$transfo=="sort"){
VarSim[is.na(VarObs)] <- NA
VarSim <- sort(VarSim, na.last = TRUE)
VarObs <- sort(VarObs, na.last = TRUE)
InputsCrit$BoolCrit <- sort(InputsCrit$BoolCrit, decreasing = TRUE)
}
##TS_ignore
TS_ignore <- !is.finite(VarObs) | !is.finite(VarSim) | !InputsCrit$BoolCrit ;
Ind_TS_ignore <- which(TS_ignore); if(length(Ind_TS_ignore)==0){ Ind_TS_ignore <- NULL; }
if(sum(!TS_ignore)==0){ OutputsCrit <- list(NA); names(OutputsCrit) <- c("CritValue"); return(OutputsCrit); }
if(sum(!TS_ignore)==1){ OutputsCrit <- list(NA); names(OutputsCrit) <- c("CritValue"); return(OutputsCrit); } ### to avoid a problem in standard deviation computation
if(inherits(OutputsModel,"hourly" )){ WarningTS <- 365; }
if(inherits(OutputsModel,"daily" )){ WarningTS <- 365; }
if(inherits(OutputsModel,"monthly")){ WarningTS <- 12; }
if(inherits(OutputsModel,"yearly" )){ WarningTS <- 3; }
if(sum(!TS_ignore)<WarningTS & warnings){ warning("\t criterion computed on less than ", WarningTS, " time-steps") }
##Other_variables_preparation
meanVarObs <- mean(VarObs[!TS_ignore]);
meanVarSim <- mean(VarSim[!TS_ignore]);
iCrit <- 0
SubCritPrint <- NULL
SubCritNames <- NULL
SubCritValues <- NULL
##SubErrorCrit_____KGE_rPearson__________________
iCrit <- iCrit+1;
SubCritPrint[iCrit] <- paste(CritName," cor(sim, obs, \"pearson\") =", sep = "")
SubCritNames[iCrit] <- "r"
SubCritValues[iCrit] <- NA;
Numer <- sum( (VarObs[!TS_ignore]-meanVarObs)*(VarSim[!TS_ignore]-meanVarSim) );
Deno1 <- sqrt( sum((VarObs[!TS_ignore]-meanVarObs)^2) );
Deno2 <- sqrt( sum((VarSim[!TS_ignore]-meanVarSim)^2) );
if(Numer==0){ if(Deno1==0 & Deno2==0){ Crit <- 1; } else { Crit <- 0; }
} else { Crit <- Numer/(Deno1*Deno2); }
if(is.numeric(Crit) & is.finite(Crit)){ SubCritValues[iCrit] <- Crit; }
##SubErrorCrit_____KGE_gamma______________________
iCrit <- iCrit+1;
SubCritPrint[iCrit] <- paste(CritName," cv(sim)/cv(obs) =", sep = "")
SubCritNames[iCrit] <- "gamma"
SubCritValues[iCrit] <- NA;
if(meanVarSim==0){ if(sd(VarSim[!TS_ignore])==0){ CVsim <- 1; } else { CVsim <- 99999; } } else { CVsim <- sd(VarSim[!TS_ignore])/meanVarSim; }
if(meanVarObs==0){ if(sd(VarObs[!TS_ignore])==0){ CVobs <- 1; } else { CVobs <- 99999; } } else { CVobs <- sd(VarObs[!TS_ignore])/meanVarObs; }
if(CVsim==0 & CVobs==0){ Crit <- 1; } else { Crit <- CVsim/CVobs ; }
if(is.numeric(Crit) & is.finite(Crit)){ SubCritValues[iCrit] <- Crit; }
##SubErrorCrit_____KGE_beta______________________
iCrit <- iCrit+1;
SubCritPrint[iCrit] <- paste(CritName," mean(sim)/mean(obs) =", sep = "")
SubCritNames[iCrit] <- "beta"
SubCritValues[iCrit] <- NA;
if(meanVarSim==0 & meanVarObs==0){ Crit <- 1; } else { Crit <- meanVarSim/meanVarObs ; }
if(is.numeric(Crit) & is.finite(Crit)){ SubCritValues[iCrit] <- Crit; }
##ErrorCrit______________________________________
if(sum(is.na(SubCritValues))==0){
CritValue <- ( 1 - sqrt( (SubCritValues[1]-1)^2 + (SubCritValues[2]-1)^2 + (SubCritValues[3]-1)^2 ) );
} }
EC <- .ErrorCrit(InputsCrit = InputsCrit, Crit = "KGE2", OutputsModel = OutputsModel, warnings = warnings)
##Verbose______________________________________
if(verbose) { CritValue <- NA
message("Crit. ", CritName, " = ", sprintf("%.4f", CritValue)) SubCritValues <- rep(NA, 3)
message(paste("\tSubCrit.", SubCritPrint, sprintf("%.4f", SubCritValues), "\n", sep = " ")) SubCritNames <- c("r", "gamma", "beta")
SubCritPrint <- rep(NA, 3)
if (EC$CritCompute) {
## Other variables preparation
meanVarObs <- mean(EC$VarObs[!EC$TS_ignore])
meanVarSim <- mean(EC$VarSim[!EC$TS_ignore])
## SubErrorCrit KGE rPearson
SubCritPrint[1L] <- paste0(EC$CritName, " cor(sim, obs, \"pearson\") =")
Numer <- sum((EC$VarObs[!EC$TS_ignore] - meanVarObs) * (EC$VarSim[!EC$TS_ignore] - meanVarSim))
Deno1 <- sqrt(sum((EC$VarObs[!EC$TS_ignore] - meanVarObs)^2))
Deno2 <- sqrt(sum((EC$VarSim[!EC$TS_ignore] - meanVarSim)^2))
if (Numer == 0) {
if (Deno1 == 0 & Deno2 == 0) {
Crit <- 1
} else {
Crit <- 0
}
} else {
Crit <- Numer / (Deno1 * Deno2)
}
if (is.numeric(Crit) & is.finite(Crit)) {
SubCritValues[1L] <- Crit
}
## SubErrorCrit KGE gamma
SubCritPrint[2L] <- paste0(EC$CritName, " cv(sim)/cv(obs) =")
if (meanVarSim == 0) {
if (sd(EC$VarSim[!EC$TS_ignore]) == 0) {
CVsim <- 1
} else {
CVsim <- 99999
}
} else {
CVsim <- sd(EC$VarSim[!EC$TS_ignore]) / meanVarSim
}
if (meanVarObs == 0) {
if (sd(EC$VarObs[!EC$TS_ignore]) == 0) {
CVobs <- 1
} else {
CVobs <- 99999
}
} else {
CVobs <- sd(EC$VarObs[!EC$TS_ignore]) / meanVarObs
}
if (CVsim == 0 &
CVobs == 0) {
Crit <- 1
} else {
Crit <- CVsim / CVobs
}
if (is.numeric(Crit) & is.finite(Crit)) {
SubCritValues[2L] <- Crit
}
## SubErrorCrit KGE beta
SubCritPrint[3L] <- paste0(EC$CritName, " mean(sim)/mean(obs) =")
if (meanVarSim == 0 & meanVarObs == 0) {
Crit <- 1
} else {
Crit <- meanVarSim / meanVarObs
}
if (is.numeric(Crit) & is.finite(Crit)) {
SubCritValues[3L] <- Crit
}
## ErrorCrit
if (sum(is.na(SubCritValues)) == 0) {
CritValue <- (1 - sqrt((SubCritValues[1L] - 1)^2 + (SubCritValues[2L] - 1)^2 + (SubCritValues[3L] - 1)^2))
}
## Verbose
if (verbose) {
message(sprintf("Crit. %s = %.4f", EC$CritName, CritValue))
message(paste("\tSubCrit.", SubCritPrint, sprintf("%.4f", SubCritValues), "\n", sep = " "))
}
} }
##Output_________________________________________
OutputsCrit <- list(CritValue = CritValue, CritName = CritName, ## Output
SubCritValues = SubCritValues, SubCritNames = SubCritNames, CritBestValue = CritBestValue, OutputsCrit <- list(CritValue = CritValue,
Multiplier = Multiplier, Ind_notcomputed = Ind_TS_ignore) CritName = EC$CritName,
SubCritValues = SubCritValues,
SubCritNames = SubCritNames,
CritBestValue = EC$CritBestValue,
Multiplier = EC$Multiplier,
Ind_notcomputed = EC$Ind_TS_ignore)
class(OutputsCrit) <- c("KGE2", "ErrorCrit")
return(OutputsCrit) return(OutputsCrit)
} }
class(ErrorCrit_KGE2) <- c("FUN_CRIT", class(ErrorCrit_KGE2))
R/ErrorCrit_NSE.R
View file @ 4bccfee2
ErrorCrit_NSE <- function(InputsCrit,OutputsModel, warnings = TRUE, verbose = TRUE){ ErrorCrit_NSE <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE) {
## Arguments check
if (!inherits(OutputsModel, "OutputsModel")) {
stop("'OutputsModel' must be of class 'OutputsModel'")
}
##Arguments_check________________________________ EC <- .ErrorCrit(InputsCrit = InputsCrit, Crit = "NSE", OutputsModel = OutputsModel, warnings = warnings)
if(inherits(InputsCrit,"InputsCrit")==FALSE){ stop("InputsCrit must be of class 'InputsCrit' \n"); return(NULL); }
if(inherits(OutputsModel,"OutputsModel")==FALSE){ stop("OutputsModel must be of class 'OutputsModel' \n"); return(NULL); }
CritValue <- NA
##Initialisation_________________________________ if (EC$CritCompute) {
CritName <- NA; ## ErrorCrit
if(InputsCrit$transfo=="" ){ CritName <- "NSE[Q]" ; } Emod <- sum((EC$VarSim[!EC$TS_ignore] - EC$VarObs[!EC$TS_ignore])^2)
if(InputsCrit$transfo=="sqrt"){ CritName <- "NSE[sqrt(Q)]"; } Eref <- sum((EC$VarObs[!EC$TS_ignore] - mean(EC$VarObs[!EC$TS_ignore]))^2)
if(InputsCrit$transfo=="log" ){ CritName <- "NSE[log(Q)]" ; }
if(InputsCrit$transfo=="inv" ){ CritName <- "NSE[1/Q]" ; }
if(InputsCrit$transfo=="sort"){ CritName <- "NSE[sort(Q)]"; }
CritValue <- NA;
CritBestValue <- +1;
Multiplier <- -1; ### must be equal to -1 or +1 only
if (Emod == 0 & Eref == 0) {
Crit <- 0
} else {
Crit <- (1 - Emod / Eref)
}
if (is.numeric(Crit) & is.finite(Crit)) {
CritValue <- Crit
}
##Data_preparation_______________________________ ## Verbose
VarObs <- InputsCrit$Qobs ; VarObs[!InputsCrit$BoolCrit] <- NA; if (verbose) {
VarSim <- OutputsModel$Qsim; VarSim[!InputsCrit$BoolCrit] <- NA; message(sprintf("Crit. %s = %.4f", EC$CritName, CritValue))
##Data_transformation }
if("Ind_zeroes" %in% names(InputsCrit) & "epsilon" %in% names(InputsCrit)){ if(length(InputsCrit$Ind_zeroes)>0){
VarObs <- VarObs + InputsCrit$epsilon;
VarSim <- VarSim + InputsCrit$epsilon;
} }
if(InputsCrit$transfo=="sqrt"){ VarObs <- sqrt(VarObs); VarSim <- sqrt(VarSim); }
if(InputsCrit$transfo=="log" ){ VarObs <- log(VarObs) ; VarSim <- log(VarSim) ; VarSim[VarSim < -1E100] <- NA; }
if(InputsCrit$transfo=="inv" ){ VarObs <- 1/VarObs ; VarSim <- 1/VarSim ; VarSim[abs(VarSim) > 1E+100] <- NA; }
if(InputsCrit$transfo=="sort"){
VarSim[is.na(VarObs)] <- NA
VarSim <- sort(VarSim, na.last = TRUE)
VarObs <- sort(VarObs, na.last = TRUE)
InputsCrit$BoolCrit <- sort(InputsCrit$BoolCrit, decreasing = TRUE)
} }
##TS_ignore
TS_ignore <- !is.finite(VarObs) | !is.finite(VarSim) | !InputsCrit$BoolCrit ;
Ind_TS_ignore <- which(TS_ignore); if(length(Ind_TS_ignore)==0){ Ind_TS_ignore <- NULL; }
if(sum(!TS_ignore)==0){ OutputsCrit <- list(NA); names(OutputsCrit) <- c("CritValue"); return(OutputsCrit); }
if(inherits(OutputsModel,"hourly" )){ WarningTS <- 365; }
if(inherits(OutputsModel,"daily" )){ WarningTS <- 365; }
if(inherits(OutputsModel,"monthly")){ WarningTS <- 12; }
if(inherits(OutputsModel,"yearly" )){ WarningTS <- 3; }
if(sum(!TS_ignore)<WarningTS & warnings){ warning("\t criterion computed on less than ", WarningTS, " time-steps") }
##Other_variables_preparation
meanVarObs <- mean(VarObs[!TS_ignore]);
meanVarSim <- mean(VarSim[!TS_ignore]);
##ErrorCrit______________________________________
Emod <- sum((VarSim[!TS_ignore]-VarObs[!TS_ignore])^2);
Eref <- sum((VarObs[!TS_ignore]-mean(VarObs[!TS_ignore]))^2);
if(Emod==0 & Eref==0){ Crit <- 0; } else { Crit <- (1-Emod/Eref); }
if(is.numeric(Crit) & is.finite(Crit)){ CritValue <- Crit; }
##Verbose______________________________________ ## Output
if(verbose) { OutputsCrit <- list(CritValue = CritValue,
message("Crit. ", CritName, " = ", sprintf("%.4f", CritValue)) CritName = EC$CritName,
} CritBestValue = EC$CritBestValue,
Multiplier = EC$Multiplier,
Ind_notcomputed = EC$Ind_TS_ignore)
##Output_________________________________________
OutputsCrit <- list(CritValue = CritValue, CritName = CritName, class(OutputsCrit) <- c("NSE", "ErrorCrit")
CritBestValue = CritBestValue,
Multiplier = Multiplier, Ind_notcomputed = Ind_TS_ignore)
return(OutputsCrit) return(OutputsCrit)
} }
class(ErrorCrit_NSE) <- c("FUN_CRIT", class(ErrorCrit_NSE))
R/ErrorCrit_RMSE.R
View file @ 4bccfee2
ErrorCrit_RMSE <- function(InputsCrit,OutputsModel, warnings = TRUE, verbose = TRUE){ ErrorCrit_RMSE <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE) {
## Arguments check
if (!inherits(OutputsModel, "OutputsModel")) {
stop("'OutputsModel' must be of class 'OutputsModel'")
}
##Arguments_check________________________________ EC <- .ErrorCrit(InputsCrit = InputsCrit, Crit = "RMSE", OutputsModel = OutputsModel, warnings = warnings)
if(inherits(InputsCrit,"InputsCrit")==FALSE){ stop("InputsCrit must be of class 'InputsCrit' \n"); return(NULL); }
if(inherits(OutputsModel,"OutputsModel")==FALSE){ stop("OutputsModel must be of class 'OutputsModel' \n"); return(NULL); }
##Initialisation_________________________________ CritValue <- NA
CritName <- NA;
if(InputsCrit$transfo=="" ){ CritName <- "RMSE[Q]" ; }
if(InputsCrit$transfo=="sqrt"){ CritName <- "RMSE[sqrt(Q)]"; }
if(InputsCrit$transfo=="log" ){ CritName <- "RMSE[log(Q)]" ; }
if(InputsCrit$transfo=="inv" ){ CritName <- "RMSE[1/Q]" ; }
if(InputsCrit$transfo=="sort"){ CritName <- "RMSE[sort(Q)]"; }
CritValue <- NA; if (EC$CritCompute) {
CritBestValue <- +1; ## ErrorCrit
Multiplier <- +1; ### must be equal to -1 or +1 only Numer <- sum((EC$VarSim - EC$VarObs)^2, na.rm = TRUE)
Denom <- sum(!is.na(EC$VarObs))
if (Numer == 0) {
Crit <- 0
} else {
Crit <- sqrt(Numer / Denom)
}
if (is.numeric(Crit) & is.finite(Crit)) {
CritValue <- Crit
}
##Data_preparation_______________________________ ## Verbose
VarObs <- InputsCrit$Qobs ; VarObs[!InputsCrit$BoolCrit] <- NA; if (verbose) {
VarSim <- OutputsModel$Qsim; VarSim[!InputsCrit$BoolCrit] <- NA; message(sprintf("Crit. %s = %.4f", EC$CritName, CritValue))
##Data_transformation }
if("Ind_zeroes" %in% names(InputsCrit) & "epsilon" %in% names(InputsCrit)){ if(length(InputsCrit$Ind_zeroes)>0){
VarObs <- VarObs + InputsCrit$epsilon;
VarSim <- VarSim + InputsCrit$epsilon;
} }
if(InputsCrit$transfo=="sqrt"){ VarObs <- sqrt(VarObs); VarSim <- sqrt(VarSim); }
if(InputsCrit$transfo=="log" ){ VarObs <- log(VarObs) ; VarSim <- log(VarSim) ; VarSim[VarSim < -1E100] <- NA; }
if(InputsCrit$transfo=="inv" ){ VarObs <- 1/VarObs ; VarSim <- 1/VarSim ; VarSim[abs(VarSim) > 1E+100] <- NA; }
if(InputsCrit$transfo=="sort"){
VarSim[is.na(VarObs)] <- NA
VarSim <- sort(VarSim, na.last = TRUE)
VarObs <- sort(VarObs, na.last = TRUE)
InputsCrit$BoolCrit <- sort(InputsCrit$BoolCrit, decreasing = TRUE)
} }
##TS_ignore
TS_ignore <- !is.finite(VarObs) | !is.finite(VarSim) | !InputsCrit$BoolCrit ;
Ind_TS_ignore <- which(TS_ignore); if(length(Ind_TS_ignore)==0){ Ind_TS_ignore <- NULL; }
if(sum(!TS_ignore)==0){ OutputsCrit <- list(NA); names(OutputsCrit) <- c("CritValue"); return(OutputsCrit); }
if(inherits(OutputsModel,"hourly" )){ WarningTS <- 365; }
if(inherits(OutputsModel,"daily" )){ WarningTS <- 365; }
if(inherits(OutputsModel,"monthly")){ WarningTS <- 12; }
if(inherits(OutputsModel,"yearly" )){ WarningTS <- 3; }
if(sum(!TS_ignore)<WarningTS & warnings){ warning("\t criterion computed on less than ", WarningTS, " time-steps") }
##ErrorCrit______________________________________ ## Output
Numer <- sum((VarSim-VarObs)^2,na.rm=TRUE); OutputsCrit <- list(CritValue = CritValue,
Denom <- sum(!is.na(VarObs)); CritName = EC$CritName,
if(Numer==0){ Crit <- 0; } else { Crit <- sqrt(Numer/Denom); } CritBestValue = EC$CritBestValue,
if(is.numeric(Crit) & is.finite(Crit)){ CritValue <- Crit; } Multiplier = EC$Multiplier,
Ind_notcomputed = EC$Ind_TS_ignore)
##Verbose______________________________________ class(OutputsCrit) <- c("RMSE", "ErrorCrit")
if(verbose) {
message("Crit. ", CritName, " = ", sprintf("%.4f", CritValue))
}
##Output_________________________________________
OutputsCrit <- list(CritValue = CritValue, CritName = CritName, CritBestValue = CritBestValue,
Multiplier = Multiplier, Ind_notcomputed = Ind_TS_ignore)
return(OutputsCrit) return(OutputsCrit)
} }
class(ErrorCrit_RMSE) <- c("FUN_CRIT", class(ErrorCrit_RMSE))
R/Imax.R 0 → 100644
View file @ 4bccfee2
Imax <- function(InputsModel,
IndPeriod_Run,
TestedValues = seq(from = 0.1, to = 3, by = 0.1)) {
## ---------- check arguments
## InputsModel
if (!inherits(InputsModel, "InputsModel")) {
stop("'InputsModel' must be of class 'InputsModel'")
}
if (!inherits(InputsModel, "hourly")) {
stop("'InputsModel' must be of class 'hourly'")
}
## IndPeriod_Run
if (!is.vector(IndPeriod_Run)) {
stop("'IndPeriod_Run' must be a vector of numeric values")
}
if (!inherits(IndPeriod_Run, "integer")) {
stop("'IndPeriod_Run' must be of type integer")
}
if (!identical(as.integer(IndPeriod_Run), IndPeriod_Run[1]:IndPeriod_Run[length(IndPeriod_Run)])) {
stop("'IndPeriod_Run' must be a continuous sequence of integers")
}
## TestedValues
if (!(is.numeric(TestedValues))) {
stop("'TestedValues' must be 'numeric'")
}
## ---------- hourly inputs aggregation
## aggregate data at the daily time step
daily_data <- SeriesAggreg(InputsModel[IndPeriod_Run], Format = "%Y%m%d")
## ---------- calculate interception
## calculate total interception of daily GR models on the period
cum_daily <- sum(pmin(daily_data$Precip, daily_data$PotEvap))
if (anyNA(cum_daily)) {
stop("'IndPeriod_Run' must be set to select 24 hours by day")
}
## calculate total interception of the GR5H interception store on the period
## and compute difference with daily values
differences <- array(NA, c(length(TestedValues)))
for (Imax in TestedValues) {
C0 <- 0
cum_hourly <- 0
for (i in IndPeriod_Run) {
Ec <- min(InputsModel$PotEvap[i], InputsModel$Precip[i] + C0)
Pth <- max(0, InputsModel$Precip[i] - (Imax-C0) - Ec)
C0 <- C0 + InputsModel$Precip[i] - Ec - Pth
cum_hourly <- cum_hourly + Ec
}
differences[which(Imax == TestedValues)] <- abs(cum_hourly - cum_daily)
}
## return the Imax value that minimises the difference
return(TestedValues[which.min(differences)])
}
R/PE_Oudin.R 0 → 100644
View file @ 4bccfee2
PE_Oudin <- function(JD, Temp,
Lat, LatUnit = c("rad", "deg"),
TimeStepIn = "daily", TimeStepOut = "daily",
RunFortran = FALSE) {
## ---------- check arguments
if (!(inherits(JD, "numeric") | inherits(JD, "integer"))) {
stop("'JD' must be of class 'numeric'")
}
if (!(inherits(Temp, "numeric") | inherits(Temp, "integer"))) {
stop("'Temp' must be of class 'numeric'")
}
if (length(JD) != length(Temp)) {
stop("'JD' and 'Temp' must have the same length")
}
if (!RunFortran & (!inherits(Lat, "numeric") | length(Lat) != 1)) {
stop("'Lat' must be a 'numeric' of length one")
}
if (RunFortran & (!inherits(Lat, "numeric") | (!length(Lat) %in% c(1, length(Temp))))) {
stop("'Lat' must be a 'numeric' of length one or of the same length as 'Temp'")
}
LatUnit <- match.arg(LatUnit, choices = c("rad", "deg"))
if (LatUnit[1L] == "rad" & (all(Lat >= pi/2) | all(Lat <= -pi/2))) {
stop("'Lat' must be comprised between -pi/2 and +pi/2 degrees")
}
if (LatUnit[1L] == "deg" & (all(Lat >= 90) | all(Lat <= -90))) {
stop("'Lat' must be comprised between -90 and +90 degrees")
}
if (LatUnit[1L] == "rad") {
FI <- Lat
}
if (LatUnit[1L] == "deg") {
FI <- Lat / (180 / pi)
}
if (any(JD < 0) | any(JD > 366)) {
stop("'JD' must only contain integers from 1 to 366")
}
TimeStep <- c("daily", "hourly")
TimeStepIn <- match.arg(TimeStepIn , choices = TimeStep)
TimeStepOut <- match.arg(TimeStepOut, choices = TimeStep)
rleJD <- rle(JD)
msgDaliy <- "each day should have only one identical value of Julian days. The time series is not sorted, or contains duplicate or missing dates"
msgHourly <- "each day must have 24 identical values of Julian days (one for each hour). The time series is not sorted, or contains duplicate or missing dates"
if (TimeStepIn == "daily" & any(rleJD$lengths != 1)) {
warning(msgDaliy)
}
## ---------- hourly inputs aggregation
if (TimeStepIn == "hourly") {
JD <- rleJD$values
idJD <- rep(seq_along(JD), each = rleJD$lengths[1L])
if (length(Temp) != length(idJD)) {
stop(msgHourly)
} else {
Temp <- as.vector(tapply(X = Temp, INDEX = idJD, FUN = mean))
}
}
if (TimeStepIn == "hourly" & any(rleJD$lengths != 24)) {
warning(msgHourly)
}
## ---------- Oudin's formula
if (RunFortran) {
LInputs = as.integer(length(Temp))
if (length(FI) == 1) {
FI <- rep(FI, LInputs)
}
RESULTS <- .Fortran("frun_pe_oudin", PACKAGE = "airGR",
##inputs
LInputs = LInputs,
InputsLAT = as.double(FI),
InputsTT = as.double(Temp),
InputsJJ = as.double(JD),
##outputs
PE_Oudin_D = rep(as.double(-99e9), LInputs)
)
PE_Oudin_D = RESULTS$PE_Oudin_D
} else {
PE_Oudin_D <- rep(NA, length(Temp))
COSFI <- cos(FI)
for (k in seq_along(Temp)) {
TETA <- 0.4093 * sin(JD[k] / 58.1 - 1.405)
COSTETA <- cos(TETA)
COSGZ <- max(0.001, cos(FI - TETA))
GZ <- acos(COSGZ)
COSOM <- 1 - COSGZ / COSFI / COSTETA
if (COSOM < -1) {
COSOM <- -1
}
if (COSOM > 1) {
COSOM <- 1
}
COSOM2 <- COSOM * COSOM
if (COSOM2 >= 1) {
SINOM <- 0
} else {
SINOM <- sqrt(1 - COSOM2)
}
OM <- acos(COSOM)
COSPZ <- COSGZ + COSFI * COSTETA * (SINOM/OM - 1)
if (COSPZ < 0.001) {
COSPZ <- 0.001
}
ETA <- 1 + cos(JD[k] / 58.1) / 30
GE <- 446 * OM * COSPZ * ETA
if (is.na(Temp[k])) {
PE_Oudin_D[k] <- NA
} else {
if (Temp[k] >= -5.0) {
PE_Oudin_D[k] <- GE * (Temp[k] + 5) / 100 / 28.5
} else {
PE_Oudin_D[k] <- 0
}
}
}
}
## ---------- disaggregate PE from daily to hourly
if (TimeStepOut == "hourly") {
parab_D2H <- c(0.000, 0.000, 0.000, 0.000, 0.000, 0.000,
0.035, 0.062, 0.079, 0.097, 0.110, 0.117,
0.117, 0.110, 0.097, 0.079, 0.062, 0.035,
0.000, 0.000, 0.000, 0.000, 0.000, 0.000)
PE_Oudin_H <- rep(PE_Oudin_D, each = 24) * rep(parab_D2H, times = length(PE_Oudin_D))
}
## ---------- outputs warnings
if (any(is.na(Temp))) {
if (any(is.na(PE_Oudin_D))) {
warning("'Temp' time series, and therefore the returned 'PE' time series, contain missing value(s)")
} else {
warning("'Temp' time series contains missing value(s)")
}
}
if (!any(is.na(Temp)) & any(is.na(PE_Oudin_D))) {
warning("returned 'PE' time series contains missing value(s)")
}
if (TimeStepOut == "daily") {
PE_Oudin <- PE_Oudin_D
} else {
PE_Oudin <- PE_Oudin_H
}
return(PE_Oudin)
}
R/PEdaily_Oudin.R deleted 100644 → 0
View file @ 10f794b5
PEdaily_Oudin <- function(JD,Temp,LatRad){
PE_Oudin_D <- rep(NA,length(Temp));
for(k in 1:length(Temp)){
FI <- LatRad ### latitude in rad
### FI <- LatDeg/(180/pi) ### conversion from deg to rad
COSFI <- cos(FI)
AFI <- abs(LatRad/42.)
TETA <- 0.4093*sin(JD[k]/58.1-1.405)
COSTETA <- cos(TETA)
COSGZ <- max(0.001,cos(FI-TETA))
GZ <- acos(COSGZ)
COSGZ2 <- COSGZ*COSGZ
if(COSGZ2 >= 1){ SINGZ <- 0. } else { SINGZ <- sqrt(1.-COSGZ2) }
COSOM <- 1.-COSGZ/COSFI/COSTETA
if(COSOM < -1.){ COSOM <- -1. }
if(COSOM > 1.){ COSOM <- 1. }
COSOM2 <- COSOM*COSOM
if(COSOM2 >= 1.){ SINOM <- 0. } else { SINOM <- sqrt(1.-COSOM2) }
OM <- acos(COSOM)
COSPZ <- COSGZ+COSFI*COSTETA*(SINOM/OM-1.)
if(COSPZ < 0.001){ COSPZ <- 0.001 }
ETA <- 1.+cos(JD[k]/58.1)/30.
GE <- 446.*OM*COSPZ*ETA
if(Temp[k] >= -5.0) { PE_Oudin_D[k] <- GE*(Temp[k]+5.)/100./28.5 } else { PE_Oudin_D[k] <- 0 }
}
return(PE_Oudin_D);
}
R/RunModel.R
View file @ 4bccfee2
RunModel <- function(InputsModel, RunOptions, Param, FUN_MOD) { RunModel <- function(InputsModel, RunOptions, Param, FUN_MOD, ...) {
return(FUN_MOD(InputsModel, RunOptions, Param))
FUN_MOD <- match.fun(FUN_MOD)
if (inherits(InputsModel, "SD") && !identical(FUN_MOD, RunModel_Lag)) {
# Lag model take one parameter at the beginning of the vector
iFirstParamRunOffModel <- 2
RunOptions$FeatFUN_MOD$NbParam <- RunOptions$FeatFUN_MOD$NbParam - 1
} else {
# All parameters
iFirstParamRunOffModel <- 1
}
OutputsModel <- FUN_MOD(InputsModel = InputsModel, RunOptions = RunOptions,
Param = Param[iFirstParamRunOffModel:length(Param)], ...)
if (inherits(InputsModel, "SD") && !identical(FUN_MOD, RunModel_Lag)) {
OutputsModel <- RunModel_Lag(InputsModel, RunOptions, Param[1], OutputsModel)
}
return(OutputsModel)
} }
R/RunModel_CemaNeige.R
View file @ 4bccfee2
RunModel_CemaNeige <- function(InputsModel,RunOptions,Param){ RunModel_CemaNeige <- function(InputsModel, RunOptions, Param) {
NParam <- 2;
FortranOutputsCemaNeige <- c("Pliq","Psol","SnowPack","ThermalState","Gratio","PotMelt","Melt","PliqAndMelt", "Temp"); ## Initialization of variables
IsHyst <- inherits(RunOptions, "hysteresis")
##Arguments_check NParamCN <- ifelse(test = IsHyst, yes = 4L, no = 2L)
if(inherits(InputsModel,"InputsModel")==FALSE){ stop("InputsModel must be of class 'InputsModel' \n"); return(NULL); } NStates <- 4L
if(inherits(InputsModel,"daily" )==FALSE){ stop("InputsModel must be of class 'daily' \n"); return(NULL); } FortranOutputsCemaNeige <- .FortranOutputs(GR = NULL, isCN = TRUE)$CN
if(inherits(InputsModel,"CemaNeige" )==FALSE){ stop("InputsModel must be of class 'CemaNeige' \n"); return(NULL); }
if(inherits(RunOptions,"RunOptions" )==FALSE){ stop("RunOptions must be of class 'RunOptions' \n"); return(NULL); }
if(inherits(RunOptions,"CemaNeige" )==FALSE){ stop("RunOptions must be of class 'CemaNeige' \n"); return(NULL); } ## Arguments check
if(!is.vector(Param)){ stop("Param must be a vector \n"); return(NULL); } if (!inherits(InputsModel, "InputsModel")) {
if(sum(!is.na(Param))!=NParam){ stop(paste("Param must be a vector of length ",NParam," and contain no NA \n",sep="")); return(NULL); } stop("'InputsModel' must be of class 'InputsModel'")
Param <- as.double(Param); }
if (!inherits(InputsModel, "daily") & !inherits(InputsModel, "hourly")) {
##Input_data_preparation stop("'InputsModel' must be of class 'daily' or 'hourly'")
if(identical(RunOptions$IndPeriod_WarmUp,0)){ RunOptions$IndPeriod_WarmUp <- NULL; } }
IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp,RunOptions$IndPeriod_Run); if (!inherits(InputsModel, "CemaNeige")) {
IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):length(IndPeriod1); stop("'InputsModel' must be of class 'CemaNeige'")
ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim; }
ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim; if (!inherits(RunOptions, "RunOptions")) {
stop("'RunOptions' must be of class 'RunOptions'")
}
if (!inherits(RunOptions, "CemaNeige")) {
##SNOW_MODULE________________________________________________________________________________## stop("'RunOptions' must be of class 'CemaNeige'")
ParamCemaNeige <- Param; }
NLayers <- length(InputsModel$LayerPrecip); if (!is.vector(Param) | !is.numeric(Param)) {
if(sum(is.na(ParamCemaNeige))!=0){ stop("Param contains missing values \n"); return(NULL); } stop("'Param' must be a numeric vector")
if("all" %in% RunOptions$Outputs_Sim){ IndOutputsCemaNeige <- as.integer(1:length(FortranOutputsCemaNeige)); }
} else { IndOutputsCemaNeige <- which(FortranOutputsCemaNeige %in% RunOptions$Outputs_Sim); } if (sum(!is.na(Param)) != NParamCN) {
CemaNeigeLayers <- list(); CemaNeigeStateEnd <- NULL; NameCemaNeigeLayers <- "CemaNeigeLayers"; stop(sprintf("'Param' must be a vector of length %i and contain no NA", NParamCN))
}
##Call_DLL_CemaNeige_________________________ if (inherits(InputsModel, "daily")) {
for(iLayer in 1:NLayers){ time_step <- "daily"
StateStartCemaNeige <- RunOptions$IniStates[ (2*(iLayer-1)+1):(2*(iLayer-1)+2) ]; time_mult <- 1
RESULTS <- .Fortran("frun_CemaNeige",PACKAGE="airGR", }
##inputs if (inherits(InputsModel, "hourly")) {
LInputs=as.integer(length(IndPeriod1)), ### length of input and output series time_step <- "hourly"
InputsPrecip=InputsModel$LayerPrecip[[iLayer]][IndPeriod1], ### input series of total precipitation [mm/d] time_mult <- 24
InputsFracSolidPrecip=InputsModel$LayerFracSolidPrecip[[iLayer]][IndPeriod1], ### input series of fraction of solid precipitation [0-1] }
InputsTemp=InputsModel$LayerTemp[[iLayer]][IndPeriod1], ### input series of air mean temperature [degC]
MeanAnSolidPrecip=RunOptions$MeanAnSolidPrecip[iLayer], ### value of annual mean solid precip [mm/y]
NParam=as.integer(2), ### number of model parameter = 2
Param=ParamCemaNeige, ### parameter set
NStates=as.integer(2), ### number of state variables used for model initialising = 2
StateStart=StateStartCemaNeige, ### state variables used when the model run starts
NOutputs=as.integer(length(IndOutputsCemaNeige)), ### number of output series
IndOutputs=IndOutputsCemaNeige, ### indices of output series
##outputs
Outputs=matrix(as.double(-999.999),nrow=length(IndPeriod1),ncol=length(IndOutputsCemaNeige)), ### output series [mm]
StateEnd=rep(as.double(-999.999),as.integer(2)) ### state variables at the end of the model run (reservoir levels [mm] and HU)
)
RESULTS$Outputs[ round(RESULTS$Outputs ,3)==(-999.999)] <- NA;
RESULTS$StateEnd[round(RESULTS$StateEnd,3)==(-999.999)] <- NA;
##Data_storage
CemaNeigeLayers[[iLayer]] <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]);
names(CemaNeigeLayers[[iLayer]]) <- FortranOutputsCemaNeige[IndOutputsCemaNeige];
if(ExportStateEnd){ CemaNeigeStateEnd <- c(CemaNeigeStateEnd,RESULTS$StateEnd); }
rm(RESULTS);
} ###ENDFOR_iLayer
names(CemaNeigeLayers) <- paste("Layer",formatC(1:NLayers,width=2,flag="0"),sep="");
##Output_data_preparation
if(ExportDatesR==FALSE & ExportStateEnd==FALSE){
OutputsModel <- list(CemaNeigeLayers);
names(OutputsModel) <- NameCemaNeigeLayers; }
if(ExportDatesR==TRUE & ExportStateEnd==FALSE){
OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
list(CemaNeigeLayers) );
names(OutputsModel) <- c("DatesR",NameCemaNeigeLayers); }
if(ExportDatesR==FALSE & ExportStateEnd==TRUE){
OutputsModel <- c( list(CemaNeigeLayers),
list(CemaNeigeStateEnd));
names(OutputsModel) <- c(NameCemaNeigeLayers,"StateEnd"); }
if(ExportDatesR==TRUE & ExportStateEnd==TRUE){
OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
list(CemaNeigeLayers),
list(CemaNeigeStateEnd));
names(OutputsModel) <- c("DatesR",NameCemaNeigeLayers,"StateEnd"); }
##End ## Input data preparation
class(OutputsModel) <- c("OutputsModel","daily","CemaNeige"); if (identical(RunOptions$IndPeriod_WarmUp, 0L)) {
return(OutputsModel); RunOptions$IndPeriod_WarmUp <- NULL
}
IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp, RunOptions$IndPeriod_Run)
IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):length(IndPeriod1)
## Output data preparation
ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim
ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim
## CemaNeige________________________________________________________________________________
ParamCemaNeige <- Param
NLayers <- length(InputsModel$LayerPrecip)
if (sum(is.na(ParamCemaNeige)) != 0) {
stop("Param contains missing values")
}
if ("all" %in% RunOptions$Outputs_Sim) {
IndOutputsCemaNeige <- 1:length(FortranOutputsCemaNeige)
} else {
IndOutputsCemaNeige <- which(FortranOutputsCemaNeige %in% RunOptions$Outputs_Sim)
}
CemaNeigeLayers <- list()
CemaNeigeStateEnd <- NULL
NameCemaNeigeLayers <- "CemaNeigeLayers"
## Call CemaNeige Fortran_________________________
for (iLayer in 1:NLayers) {
if (!IsHyst) {
StateStartCemaNeige <- RunOptions$IniStates[(7 + 20*time_mult + 40*time_mult) + c(iLayer, iLayer+NLayers)]
} else {
StateStartCemaNeige <- RunOptions$IniStates[(7 + 20*time_mult + 40*time_mult) + c(iLayer, iLayer+NLayers, iLayer+2*NLayers, iLayer+3*NLayers)]
}
RESULTS <- .Fortran("frun_cemaneige", PACKAGE = "airGR",
## inputs
LInputs = as.integer(length(IndPeriod1)), ### length of input and output series
InputsPrecip = InputsModel$LayerPrecip[[iLayer]][IndPeriod1], ### input series of total precipitation [mm/time step]
InputsFracSolidPrecip = InputsModel$LayerFracSolidPrecip[[iLayer]][IndPeriod1], ### input series of fraction of solid precipitation [0-1]
InputsTemp = InputsModel$LayerTemp[[iLayer]][IndPeriod1], ### input series of air mean temperature [degC]
MeanAnSolidPrecip = RunOptions$MeanAnSolidPrecip[iLayer], ### value of annual mean solid precip [mm/y]
NParam = as.integer(NParamCN), ### number of model parameters = 2 or 4
Param = as.double(ParamCemaNeige), ### parameter set
NStates = as.integer(NStates), ### number of state variables used for model initialising = 4
StateStart = StateStartCemaNeige, ### state variables used when the model run starts
IsHyst = as.integer(IsHyst), ### use of hysteresis
NOutputs = as.integer(length(IndOutputsCemaNeige)), ### number of output series
IndOutputs = IndOutputsCemaNeige, ### indices of output series
## outputs
Outputs = matrix(as.double(-99e9), nrow = length(IndPeriod1), ncol = length(IndOutputsCemaNeige)), ### output series [mm, mm/time step or degC]
StateEnd = rep(as.double(-99e9), NStates) ### state variables at the end of the model run
)
RESULTS$Outputs[RESULTS$Outputs <= -99e8] <- NA
RESULTS$StateEnd[RESULTS$StateEnd <= -99e8] <- NA
## Data storage
CemaNeigeLayers[[iLayer]] <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i])
names(CemaNeigeLayers[[iLayer]]) <- FortranOutputsCemaNeige[IndOutputsCemaNeige]
if (ExportStateEnd) {
CemaNeigeStateEnd <- c(CemaNeigeStateEnd, RESULTS$StateEnd)
}
rm(RESULTS)
} ### ENDFOR iLayer
names(CemaNeigeLayers) <- sprintf("Layer%02i", seq_len(NLayers))
if (ExportStateEnd) {
idNStates <- seq_len(NStates*NLayers) %% NStates
CemaNeigeStateEnd <- CreateIniStates(FUN_MOD = RunModel_CemaNeige, InputsModel = InputsModel, IsHyst = IsHyst,
ProdStore = NULL, RoutStore = NULL, ExpStore = NULL,
UH1 = NULL, UH2 = NULL,
GCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 1]],
eTGCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 2]],
GthrCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 3]],
GlocmaxCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 0]],
verbose = FALSE)
}
## Output data preparation
## OutputsModel only
if (!ExportDatesR & !ExportStateEnd) {
OutputsModel <- list(CemaNeigeLayers)
names(OutputsModel) <- NameCemaNeigeLayers
}
## DatesR and OutputsModel only
if (ExportDatesR & !ExportStateEnd) {
OutputsModel <- c(list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
list(CemaNeigeLayers))
names(OutputsModel) <- c("DatesR", NameCemaNeigeLayers)
}
## OutputsModel and StateEnd only
if (!ExportDatesR & ExportStateEnd) {
OutputsModel <- c(list(CemaNeigeLayers),
list(CemaNeigeStateEnd))
names(OutputsModel) <- c(NameCemaNeigeLayers, "StateEnd")
}
## DatesR and OutputsModel and StateEnd
if (ExportDatesR & ExportStateEnd) {
OutputsModel <- c(list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
list(CemaNeigeLayers),
list(CemaNeigeStateEnd))
names(OutputsModel) <- c("DatesR", NameCemaNeigeLayers, "StateEnd")
}
## End
class(OutputsModel) <- c("OutputsModel", time_step, "CemaNeige")
if (IsHyst) {
class(OutputsModel) <- c(class(OutputsModel), "hysteresis")
}
return(OutputsModel)
} }
R/RunModel_CemaNeigeGR4H.R 0 → 100644
View file @ 4bccfee2
RunModel_CemaNeigeGR4H <- function(InputsModel, RunOptions, Param) {
## Initialization of variables
IsHyst <- inherits(RunOptions, "hysteresis")
NParamCN <- RunOptions$FeatFUN_MOD$NbParam - 4L
NStates <- 4L
.ArgumentsCheckGR(InputsModel, RunOptions, Param)
Param <- as.double(Param)
Param_X1X3_threshold <- 1e-2
Param_X4_threshold <- 0.5
if (Param[1L] < Param_X1X3_threshold) {
warning(sprintf("Param[1] (X1: production store capacity [mm]) < %.2f\n X1 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
Param[1L] <- Param_X1X3_threshold
}
if (Param[3L] < Param_X1X3_threshold) {
warning(sprintf("Param[3] (X3: routing store capacity [mm]) < %.2f\n X3 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
Param[3L] <- Param_X1X3_threshold
}
if (Param[4L] < Param_X4_threshold) {
warning(sprintf("Param[4] (X4: unit hydrograph time constant [h]) < %.2f\n X4 set to %.2f", Param_X4_threshold, Param_X4_threshold))
Param[4L] <- Param_X4_threshold
}
## Input data preparation
if (identical(RunOptions$IndPeriod_WarmUp, 0L)) {
RunOptions$IndPeriod_WarmUp <- NULL
}
IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp, RunOptions$IndPeriod_Run)
LInputSeries <- as.integer(length(IndPeriod1))
IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries
ParamCemaNeige <- Param[(length(Param) - 1 - 2 * as.integer(IsHyst)):length(Param)]
NParamMod <- as.integer(length(Param) - (2 + 2 * as.integer(IsHyst)))
ParamMod <- Param[1:NParamMod]
NLayers <- length(InputsModel$LayerPrecip)
NStatesMod <- as.integer(length(RunOptions$IniStates) - NStates * NLayers)
## Output data preparation
ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim
ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim
## CemaNeige________________________________________________________________________________
if (inherits(RunOptions, "CemaNeige")) {
if ("all" %in% RunOptions$Outputs_Sim) {
IndOutputsCemaNeige <- as.integer(1:length(RunOptions$FortranOutputs$CN))
} else {
IndOutputsCemaNeige <- which(RunOptions$FortranOutputs$CN %in% RunOptions$Outputs_Sim)
}
CemaNeigeLayers <- list()
CemaNeigeStateEnd <- NULL
NameCemaNeigeLayers <- "CemaNeigeLayers"
## Call CemaNeige Fortran_________________________
for (iLayer in 1:NLayers) {
if (!IsHyst) {
StateStartCemaNeige <- RunOptions$IniStates[(7 + 20*24 + 40*24) + c(iLayer, iLayer+NLayers)]
} else {
StateStartCemaNeige <- RunOptions$IniStates[(7 + 20*24 + 40*24) + c(iLayer, iLayer+NLayers, iLayer+2*NLayers, iLayer+3*NLayers)]
}
RESULTS <- .Fortran("frun_cemaneige", PACKAGE = "airGR",
## inputs
LInputs = LInputSeries, ### length of input and output series
InputsPrecip = InputsModel$LayerPrecip[[iLayer]][IndPeriod1], ### input series of total precipitation [mm/h]
InputsFracSolidPrecip = InputsModel$LayerFracSolidPrecip[[iLayer]][IndPeriod1], ### input series of fraction of solid precipitation [0-1]
InputsTemp = InputsModel$LayerTemp[[iLayer]][IndPeriod1], ### input series of air mean temperature [degC]
MeanAnSolidPrecip = RunOptions$MeanAnSolidPrecip[iLayer], ### value of annual mean solid precip [mm/y]
NParam = as.integer(NParamCN), ### number of model parameters = 2 or 4
Param = as.double(ParamCemaNeige), ### parameter set
NStates = as.integer(NStates), ### number of state variables used for model initialising = 4
StateStart = StateStartCemaNeige, ### state variables used when the model run starts
IsHyst = as.integer(IsHyst), ### use of hysteresis
NOutputs = as.integer(length(IndOutputsCemaNeige)), ### number of output series
IndOutputs = IndOutputsCemaNeige, ### indices of output series
## outputs
Outputs = matrix(as.double(-99e9), nrow = LInputSeries, ncol = length(IndOutputsCemaNeige)), ### output series [mm, mm/h or degC]
StateEnd = rep(as.double(-99e9), as.integer(NStates)) ### state variables at the end of the model run
)
RESULTS$Outputs[RESULTS$Outputs <= -99e8] <- NA
RESULTS$StateEnd[RESULTS$StateEnd <= -99e8] <- NA
## Data storage
CemaNeigeLayers[[iLayer]] <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i])
names(CemaNeigeLayers[[iLayer]]) <- RunOptions$FortranOutputs$CN[IndOutputsCemaNeige]
IndPliqAndMelt <- which(names(CemaNeigeLayers[[iLayer]]) == "PliqAndMelt")
if (iLayer == 1) {
CatchMeltAndPliq <- RESULTS$Outputs[, IndPliqAndMelt] / NLayers
}
if (iLayer > 1) {
CatchMeltAndPliq <- CatchMeltAndPliq + RESULTS$Outputs[, IndPliqAndMelt] / NLayers
}
if (ExportStateEnd) {
CemaNeigeStateEnd <- c(CemaNeigeStateEnd, RESULTS$StateEnd)
}
rm(RESULTS)
} ### ENDFOR iLayer
names(CemaNeigeLayers) <- sprintf("Layer%02i", seq_len(NLayers))
} ### ENDIF RunSnowModule
if (!inherits(RunOptions, "CemaNeige")) {
CemaNeigeLayers <- list()
CemaNeigeStateEnd <- NULL
NameCemaNeigeLayers <- NULL
CatchMeltAndPliq <- InputsModel$Precip[IndPeriod1]
}
## GR model
if ("all" %in% RunOptions$Outputs_Sim) {
IndOutputsMod <- as.integer(1:length(RunOptions$FortranOutputs$GR))
} else {
IndOutputsMod <- which(RunOptions$FortranOutputs$GR %in% RunOptions$Outputs_Sim)
}
## Use of IniResLevels
if (!is.null(RunOptions$IniResLevels)) {
RunOptions$IniStates[1] <- RunOptions$IniResLevels[1] * ParamMod[1] ### production store level (mm)
RunOptions$IniStates[2] <- RunOptions$IniResLevels[2] * ParamMod[3] ### routing store level (mm)
}
## Call GR model Fortan
RESULTS <- .Fortran("frun_gr4h", PACKAGE = "airGR",
## inputs
LInputs = LInputSeries, ### length of input and output series
InputsPrecip = CatchMeltAndPliq, ### input series of total precipitation [mm/h]
InputsPE = InputsModel$PotEvap[IndPeriod1], ### input series potential evapotranspiration [mm/h]
NParam = NParamMod, ### number of model parameter
Param = ParamMod, ### parameter set
NStates = NStatesMod, ### number of state variables used for model initialising
StateStart = RunOptions$IniStates[1:NStatesMod], ### state variables used when the model run starts
NOutputs = as.integer(length(IndOutputsMod)), ### number of output series
IndOutputs = IndOutputsMod, ### indices of output series
## outputs
Outputs = matrix(as.double(-99e9), nrow = LInputSeries, ncol = length(IndOutputsMod)), ### output series [mm or mm/h]
StateEnd = rep(as.double(-99e9), NStatesMod) ### state variables at the end of the model run
)
RESULTS$Outputs[RESULTS$Outputs <= -99e8] <- NA
RESULTS$StateEnd[RESULTS$StateEnd <= -99e8] <- NA
if (ExportStateEnd) {
RESULTS$StateEnd[-3L] <- ifelse(RESULTS$StateEnd[-3L] < 0, 0, RESULTS$StateEnd[-3L]) ### remove negative values except for the ExpStore location
idNStates <- seq_len(NStates*NLayers) %% NStates
RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_CemaNeigeGR4H, InputsModel = InputsModel, IsHyst = IsHyst,
ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL,
UH1 = RESULTS$StateEnd[(1:(20*24)) + 7],
UH2 = RESULTS$StateEnd[(1:(40*24)) + (7+20*24)],
GCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 1]],
eTGCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 2]],
GthrCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 3]],
GlocmaxCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 0]],
verbose = FALSE)
}
if (inherits(RunOptions, "CemaNeige") & "Precip" %in% RunOptions$Outputs_Sim) {
RESULTS$Outputs[, which(RunOptions$FortranOutputs$GR[IndOutputsMod] == "Precip")] <-
InputsModel$Precip[IndPeriod1]
}
## OutputsModel generation
.GetOutputsModelGR(InputsModel,
RunOptions,
RESULTS,
LInputSeries,
Param,
CemaNeigeLayers)
}
R/RunModel_CemaNeigeGR4J.R
View file @ 4bccfee2
RunModel_CemaNeigeGR4J <- function(InputsModel,RunOptions,Param){ RunModel_CemaNeigeGR4J <- function(InputsModel, RunOptions, Param) {
NParam <- 6;
FortranOutputsCemaNeige <- c("Pliq","Psol","SnowPack","ThermalState","Gratio","PotMelt","Melt","PliqAndMelt", "Temp"); ## Initialization of variables
FortranOutputsMod <- c("PotEvap","Precip","Prod","AE","Perc","PR","Q9","Q1","Rout","Exch","AExch","QR","QD","Qsim"); IsHyst <- inherits(RunOptions, "hysteresis")
NParamCN <- RunOptions$FeatFUN_MOD$NbParam - 4L
##Arguments_check NStates <- 4L
if(inherits(InputsModel,"InputsModel")==FALSE){ stop("InputsModel must be of class 'InputsModel' \n"); return(NULL); }
if(inherits(InputsModel,"daily" )==FALSE){ stop("InputsModel must be of class 'daily' \n"); return(NULL); }
if(inherits(InputsModel,"GR" )==FALSE){ stop("InputsModel must be of class 'GR' \n"); return(NULL); } .ArgumentsCheckGR(InputsModel, RunOptions, Param)
if(inherits(InputsModel,"CemaNeige" )==FALSE){ stop("InputsModel must be of class 'CemaNeige' \n"); return(NULL); }
if(inherits(RunOptions,"RunOptions" )==FALSE){ stop("RunOptions must be of class 'RunOptions' \n"); return(NULL); } Param <- as.double(Param)
if(inherits(RunOptions,"GR" )==FALSE){ stop("RunOptions must be of class 'GR' \n"); return(NULL); }
if(inherits(RunOptions,"CemaNeige" )==FALSE){ stop("RunOptions must be of class 'CemaNeige' \n"); return(NULL); }
if(!is.vector(Param)){ stop("Param must be a vector \n"); return(NULL); } Param_X1X3_threshold <- 1e-2
if(sum(!is.na(Param))!=NParam){ stop(paste("Param must be a vector of length ",NParam," and contain no NA \n",sep="")); return(NULL); } Param_X4_threshold <- 0.5
Param <- as.double(Param); if (Param[1L] < Param_X1X3_threshold) {
warning(sprintf("Param[1] (X1: production store capacity [mm]) < %.2f\n X1 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
Param_X1X3_threshold <- 1e-2 Param[1L] <- Param_X1X3_threshold
if (Param[1L] < Param_X1X3_threshold) { }
warning(sprintf("Param[1] (X1: production store capacity [mm]) < %.2f\n X1 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold)) if (Param[3L] < Param_X1X3_threshold) {
Param[1L] <- Param_X1X3_threshold warning(sprintf("Param[3] (X3: routing store capacity [mm]) < %.2f\n X3 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
Param[3L] <- Param_X1X3_threshold
}
if (Param[4L] < Param_X4_threshold) {
warning(sprintf("Param[4] (X4: unit hydrograph time constant [d]) < %.2f\n X4 set to %.2f", Param_X4_threshold, Param_X4_threshold))
Param[4L] <- Param_X4_threshold
}
## Input data preparation
if (identical(RunOptions$IndPeriod_WarmUp, 0L)) {
RunOptions$IndPeriod_WarmUp <- NULL
}
IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp, RunOptions$IndPeriod_Run)
LInputSeries <- as.integer(length(IndPeriod1))
IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries
ParamCemaNeige <- Param[(length(Param) - 1 - 2 * as.integer(IsHyst)):length(Param)]
NParamMod <- as.integer(length(Param) - (2 + 2 * as.integer(IsHyst)))
ParamMod <- Param[1:NParamMod]
NLayers <- length(InputsModel$LayerPrecip)
NStatesMod <- as.integer(length(RunOptions$IniStates) - NStates * NLayers)
## Output data preparation
ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim
ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim
## CemaNeige________________________________________________________________________________
if (inherits(RunOptions, "CemaNeige")) {
if ("all" %in% RunOptions$Outputs_Sim) {
IndOutputsCemaNeige <- as.integer(1:length(RunOptions$FortranOutputs$CN))
} else {
IndOutputsCemaNeige <- which(RunOptions$FortranOutputs$CN %in% RunOptions$Outputs_Sim)
}
CemaNeigeLayers <- list()
CemaNeigeStateEnd <- NULL
NameCemaNeigeLayers <- "CemaNeigeLayers"
## Call CemaNeige Fortran_________________________
for (iLayer in 1:NLayers) {
if (!IsHyst) {
StateStartCemaNeige <- RunOptions$IniStates[(7 + 20 + 40) + c(iLayer, iLayer+NLayers)]
} else {
StateStartCemaNeige <- RunOptions$IniStates[(7 + 20 + 40) + c(iLayer, iLayer+NLayers, iLayer+2*NLayers, iLayer+3*NLayers)]
} }
if (Param[3L] < Param_X1X3_threshold) { RESULTS <- .Fortran("frun_cemaneige", PACKAGE = "airGR",
warning(sprintf("Param[3] (X3: routing store capacity [mm]) < %.2f\n X3 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold)) ## inputs
Param[3L] <- Param_X1X3_threshold LInputs = LInputSeries, ### length of input and output series
InputsPrecip = InputsModel$LayerPrecip[[iLayer]][IndPeriod1], ### input series of total precipitation [mm/d]
InputsFracSolidPrecip = InputsModel$LayerFracSolidPrecip[[iLayer]][IndPeriod1], ### input series of fraction of solid precipitation [0-1]
InputsTemp = InputsModel$LayerTemp[[iLayer]][IndPeriod1], ### input series of air mean temperature [degC]
MeanAnSolidPrecip = RunOptions$MeanAnSolidPrecip[iLayer], ### value of annual mean solid precip [mm/y]
NParam = as.integer(NParamCN), ### number of model parameters = 2 or 4
Param = as.double(ParamCemaNeige), ### parameter set
NStates = as.integer(NStates), ### number of state variables used for model initialising = 4
StateStart = StateStartCemaNeige, ### state variables used when the model run starts
IsHyst = as.integer(IsHyst), ### use of hysteresis
NOutputs = as.integer(length(IndOutputsCemaNeige)), ### number of output series
IndOutputs = IndOutputsCemaNeige, ### indices of output series
## outputs
Outputs = matrix(as.double(-99e9), nrow = LInputSeries, ncol = length(IndOutputsCemaNeige)), ### output series [mm, mm/d or degC]
StateEnd = rep(as.double(-99e9), as.integer(NStates)) ### state variables at the end of the model run
)
RESULTS$Outputs[RESULTS$Outputs <= -99e8] <- NA
RESULTS$StateEnd[RESULTS$StateEnd <= -99e8] <- NA
## Data storage
CemaNeigeLayers[[iLayer]] <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i])
names(CemaNeigeLayers[[iLayer]]) <- RunOptions$FortranOutputs$CN[IndOutputsCemaNeige]
IndPliqAndMelt <- which(names(CemaNeigeLayers[[iLayer]]) == "PliqAndMelt")
if (iLayer == 1) {
CatchMeltAndPliq <- RESULTS$Outputs[, IndPliqAndMelt] / NLayers
} }
if (iLayer > 1) {
##Input_data_preparation CatchMeltAndPliq <- CatchMeltAndPliq + RESULTS$Outputs[, IndPliqAndMelt] / NLayers
if(identical(RunOptions$IndPeriod_WarmUp,as.integer(0))){ RunOptions$IndPeriod_WarmUp <- NULL; }
IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp,RunOptions$IndPeriod_Run);
LInputSeries <- as.integer(length(IndPeriod1))
IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries;
ParamCemaNeige <- Param[(length(Param)-1):length(Param)];
NParamMod <- as.integer(length(Param)-2);
ParamMod <- Param[1:NParamMod];
NLayers <- length(InputsModel$LayerPrecip);
NStatesMod <- as.integer(length(RunOptions$IniStates)-2*NLayers);
ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim;
ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim;
##SNOW_MODULE________________________________________________________________________________##
if(RunOptions$RunSnowModule==TRUE){
if("all" %in% RunOptions$Outputs_Sim){ IndOutputsCemaNeige <- as.integer(1:length(FortranOutputsCemaNeige));
} else { IndOutputsCemaNeige <- which(FortranOutputsCemaNeige %in% RunOptions$Outputs_Sim); }
CemaNeigeLayers <- list(); CemaNeigeStateEnd <- NULL; NameCemaNeigeLayers <- "CemaNeigeLayers";
##Call_DLL_CemaNeige_________________________
for(iLayer in 1:NLayers){
StateStartCemaNeige <- RunOptions$IniStates[ (NStatesMod+2*(iLayer-1)+1):(NStatesMod+2*(iLayer-1)+2) ];
RESULTS <- .Fortran("frun_CemaNeige",PACKAGE="airGR",
##inputs
LInputs=LInputSeries, ### length of input and output series
InputsPrecip=InputsModel$LayerPrecip[[iLayer]][IndPeriod1], ### input series of total precipitation [mm/d]
InputsFracSolidPrecip=InputsModel$LayerFracSolidPrecip[[iLayer]][IndPeriod1], ### input series of fraction of solid precipitation [0-1]
InputsTemp=InputsModel$LayerTemp[[iLayer]][IndPeriod1], ### input series of air mean temperature [degC]
MeanAnSolidPrecip=RunOptions$MeanAnSolidPrecip[iLayer], ### value of annual mean solid precip [mm/y]
NParam=as.integer(2), ### number of model parameter = 2
Param=ParamCemaNeige, ### parameter set
NStates=as.integer(2), ### number of state variables used for model initialising = 2
StateStart=StateStartCemaNeige, ### state variables used when the model run starts
NOutputs=as.integer(length(IndOutputsCemaNeige)), ### number of output series
IndOutputs=IndOutputsCemaNeige, ### indices of output series
##outputs
Outputs=matrix(as.double(-999.999),nrow=LInputSeries,ncol=length(IndOutputsCemaNeige)), ### output series [mm]
StateEnd=rep(as.double(-999.999),as.integer(2)) ### state variables at the end of the model run (reservoir levels [mm] and HU)
)
RESULTS$Outputs[ round(RESULTS$Outputs ,3)==(-999.999)] <- NA;
RESULTS$StateEnd[round(RESULTS$StateEnd,3)==(-999.999)] <- NA;
##Data_storage
CemaNeigeLayers[[iLayer]] <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]);
names(CemaNeigeLayers[[iLayer]]) <- FortranOutputsCemaNeige[IndOutputsCemaNeige];
IndPliqAndMelt <- which(names(CemaNeigeLayers[[iLayer]]) == "PliqAndMelt");
if(iLayer==1){ CatchMeltAndPliq <- RESULTS$Outputs[,IndPliqAndMelt]/NLayers; }
if(iLayer >1){ CatchMeltAndPliq <- CatchMeltAndPliq + RESULTS$Outputs[,IndPliqAndMelt]/NLayers; }
if(ExportStateEnd){ CemaNeigeStateEnd <- c(CemaNeigeStateEnd,RESULTS$StateEnd); }
rm(RESULTS);
} ###ENDFOR_iLayer
names(CemaNeigeLayers) <- paste("Layer",formatC(1:NLayers,width=2,flag="0"),sep="");
} ###ENDIF_RunSnowModule
if(RunOptions$RunSnowModule==FALSE){
CemaNeigeLayers <- list(); CemaNeigeStateEnd <- NULL; NameCemaNeigeLayers <- NULL;
CatchMeltAndPliq <- InputsModel$Precip[IndPeriod1]; }
##MODEL______________________________________________________________________________________##
if("all" %in% RunOptions$Outputs_Sim){ IndOutputsMod <- as.integer(1:length(FortranOutputsMod));
} else { IndOutputsMod <- which(FortranOutputsMod %in% RunOptions$Outputs_Sim); }
##Use_of_IniResLevels
if("IniResLevels" %in% names(RunOptions)){
RunOptions$IniStates[1] <- RunOptions$IniResLevels[1]*ParamMod[1]; ### production store level (mm)
RunOptions$IniStates[2] <- RunOptions$IniResLevels[2]*ParamMod[3]; ### routing store level (mm)
} }
if (ExportStateEnd) {
##Call_fortan CemaNeigeStateEnd <- c(CemaNeigeStateEnd, RESULTS$StateEnd)
RESULTS <- .Fortran("frun_GR4J",PACKAGE="airGR", }
##inputs rm(RESULTS)
LInputs=LInputSeries, ### length of input and output series } ### ENDFOR iLayer
InputsPrecip=CatchMeltAndPliq, ### input series of total precipitation [mm/d] names(CemaNeigeLayers) <- sprintf("Layer%02i", seq_len(NLayers))
InputsPE=InputsModel$PotEvap[IndPeriod1], ### input series potential evapotranspiration [mm/d] } ### ENDIF RunSnowModule
NParam=NParamMod, ### number of model parameter if (!inherits(RunOptions, "CemaNeige")) {
Param=ParamMod, ### parameter set CemaNeigeLayers <- list()
NStates=NStatesMod, ### number of state variables used for model initialising CemaNeigeStateEnd <- NULL
StateStart=RunOptions$IniStates[1:NStatesMod], ### state variables used when the model run starts NameCemaNeigeLayers <- NULL
NOutputs=as.integer(length(IndOutputsMod)), ### number of output series CatchMeltAndPliq <- InputsModel$Precip[IndPeriod1]
IndOutputs=IndOutputsMod, ### indices of output series }
##outputs
Outputs=matrix(as.double(-999.999),nrow=LInputSeries,ncol=length(IndOutputsMod)), ### output series [mm]
StateEnd=rep(as.double(-999.999),NStatesMod) ### state variables at the end of the model run
) ## GR model______________________________________________________________________________________
RESULTS$Outputs[ round(RESULTS$Outputs ,3)==(-999.999)] <- NA; if ("all" %in% RunOptions$Outputs_Sim) {
RESULTS$StateEnd[round(RESULTS$StateEnd,3)==(-999.999)] <- NA; IndOutputsMod <- as.integer(1:length(RunOptions$FortranOutputs$GR))
if(RunOptions$RunSnowModule & "Precip" %in% RunOptions$Outputs_Sim){ RESULTS$Outputs[,which(FortranOutputsMod[IndOutputsMod]=="Precip")] <- InputsModel$Precip[IndPeriod1]; } } else {
IndOutputsMod <- which(RunOptions$FortranOutputs$GR %in% RunOptions$Outputs_Sim)
##Output_data_preparation }
##OutputsModel_only
if(ExportDatesR==FALSE & ExportStateEnd==FALSE){ ## Use of IniResLevels
OutputsModel <- c( lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]), if (!is.null(RunOptions$IniResLevels)) {
list(CemaNeigeLayers) ); RunOptions$IniStates[1] <- RunOptions$IniResLevels[1] * ParamMod[1] ### production store level (mm)
names(OutputsModel) <- c(FortranOutputsMod[IndOutputsMod],NameCemaNeigeLayers); } RunOptions$IniStates[2] <- RunOptions$IniResLevels[2] * ParamMod[3] ### routing store level (mm)
##DatesR_and_OutputsModel_only }
if(ExportDatesR==TRUE & ExportStateEnd==FALSE){
OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]), ## Call GR model Fortan
lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]), RESULTS <- .Fortran("frun_gr4j", PACKAGE = "airGR",
list(CemaNeigeLayers) ); ## inputs
names(OutputsModel) <- c("DatesR",FortranOutputsMod[IndOutputsMod],NameCemaNeigeLayers); } LInputs = LInputSeries, ### length of input and output series
##OutputsModel_and_SateEnd_only InputsPrecip = CatchMeltAndPliq, ### input series of total precipitation [mm/d]
if(ExportDatesR==FALSE & ExportStateEnd==TRUE){ InputsPE = InputsModel$PotEvap[IndPeriod1], ### input series potential evapotranspiration [mm/d]
OutputsModel <- c( lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]), NParam = NParamMod, ### number of model parameter
list(CemaNeigeLayers), Param = ParamMod, ### parameter set
list(c(RESULTS$StateEnd,CemaNeigeStateEnd)) ); NStates = NStatesMod, ### number of state variables used for model initialising
names(OutputsModel) <- c(FortranOutputsMod[IndOutputsMod],NameCemaNeigeLayers,"StateEnd"); } StateStart = RunOptions$IniStates[1:NStatesMod], ### state variables used when the model run starts
##DatesR_and_OutputsModel_and_SateEnd NOutputs = as.integer(length(IndOutputsMod)), ### number of output series
if(ExportDatesR==TRUE & ExportStateEnd==TRUE){ IndOutputs = IndOutputsMod, ### indices of output series
OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]), ## outputs
lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]), Outputs = matrix(as.double(-99e9), nrow = LInputSeries, ncol = length(IndOutputsMod)), ### output series [mm or mm/d]
list(CemaNeigeLayers), StateEnd = rep(as.double(-99e9), NStatesMod) ### state variables at the end of the model run
list(c(RESULTS$StateEnd,CemaNeigeStateEnd)) ); )
names(OutputsModel) <- c("DatesR",FortranOutputsMod[IndOutputsMod],NameCemaNeigeLayers,"StateEnd"); } RESULTS$Outputs[RESULTS$Outputs <= -99e8] <- NA
RESULTS$StateEnd[RESULTS$StateEnd <= -99e8] <- NA
##End if (ExportStateEnd) {
rm(RESULTS); RESULTS$StateEnd[-3L] <- ifelse(RESULTS$StateEnd[-3L] < 0, 0, RESULTS$StateEnd[-3L]) ### remove negative values except for the ExpStore location
class(OutputsModel) <- c("OutputsModel","daily","GR","CemaNeige"); idNStates <- seq_len(NStates*NLayers) %% NStates
return(OutputsModel); RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_CemaNeigeGR4J, InputsModel = InputsModel, IsHyst = IsHyst,
ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL,
UH1 = RESULTS$StateEnd[(1:20) + 7],
UH2 = RESULTS$StateEnd[(1:40) + (7+20)],
GCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 1]],
eTGCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 2]],
GthrCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 3]],
GlocmaxCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 0]],
verbose = FALSE)
}
if (inherits(RunOptions, "CemaNeige") & "Precip" %in% RunOptions$Outputs_Sim) {
RESULTS$Outputs[, which(RunOptions$FortranOutputs$GR[IndOutputsMod] == "Precip")] <-
InputsModel$Precip[IndPeriod1]
}
## OutputsModel generation
.GetOutputsModelGR(InputsModel,
RunOptions,
RESULTS,
LInputSeries,
Param,
CemaNeigeLayers)
} }
R/RunModel_CemaNeigeGR5H.R 0 → 100644
View file @ 4bccfee2
RunModel_CemaNeigeGR5H <- function(InputsModel, RunOptions, Param) {
## Initialization of variables
IsHyst <- inherits(RunOptions, "hysteresis")
NParamCN <- RunOptions$FeatFUN_MOD$NbParam - 5L
NStates <- 4L
IsIntStore <- inherits(RunOptions, "interception")
if (IsIntStore) {
Imax <- RunOptions$Imax
} else {
Imax <- -99
}
.ArgumentsCheckGR(InputsModel, RunOptions, Param)
Param <- as.double(Param)
Param_X1X3_threshold <- 1e-2
Param_X4_threshold <- 0.5
if (Param[1L] < Param_X1X3_threshold) {
warning(sprintf("Param[1] (X1: production store capacity [mm]) < %.2f\n X1 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
Param[1L] <- Param_X1X3_threshold
}
if (Param[3L] < Param_X1X3_threshold) {
warning(sprintf("Param[3] (X3: routing store capacity [mm]) < %.2f\n X3 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
Param[3L] <- Param_X1X3_threshold
}
if (Param[4L] < Param_X4_threshold) {
warning(sprintf("Param[4] (X4: unit hydrograph time constant [h]) < %.2f\n X4 set to %.2f", Param_X4_threshold, Param_X4_threshold))
Param[4L] <- Param_X4_threshold
}
## Input data preparation
if (identical(RunOptions$IndPeriod_WarmUp, 0L)) {
RunOptions$IndPeriod_WarmUp <- NULL
}
IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp, RunOptions$IndPeriod_Run)
LInputSeries <- as.integer(length(IndPeriod1))
IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries
ParamCemaNeige <- Param[(length(Param) - 1 - 2 * as.integer(IsHyst)):length(Param)]
NParamMod <- as.integer(length(Param) - (2 + 2 * as.integer(IsHyst)))
ParamMod <- Param[1:NParamMod]
NLayers <- length(InputsModel$LayerPrecip)
NStatesMod <- as.integer(length(RunOptions$IniStates) - NStates * NLayers)
## Output data preparation
ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim
ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim
## CemaNeige________________________________________________________________________________
if (inherits(RunOptions, "CemaNeige")) {
if ("all" %in% RunOptions$Outputs_Sim) {
IndOutputsCemaNeige <- as.integer(1:length(RunOptions$FortranOutputs$CN))
} else {
IndOutputsCemaNeige <- which(RunOptions$FortranOutputs$CN %in% RunOptions$Outputs_Sim)
}
CemaNeigeLayers <- list()
CemaNeigeStateEnd <- NULL
NameCemaNeigeLayers <- "CemaNeigeLayers"
## Call CemaNeige Fortran_________________________
for (iLayer in 1:NLayers) {
if (!IsHyst) {
StateStartCemaNeige <- RunOptions$IniStates[(7 + 20*24 + 40*24) + c(iLayer, iLayer+NLayers)]
} else {
StateStartCemaNeige <- RunOptions$IniStates[(7 + 20*24 + 40*24) + c(iLayer, iLayer+NLayers, iLayer+2*NLayers, iLayer+3*NLayers)]
}
RESULTS <- .Fortran("frun_cemaneige", PACKAGE = "airGR",
## inputs
LInputs = LInputSeries, ### length of input and output series
InputsPrecip = InputsModel$LayerPrecip[[iLayer]][IndPeriod1], ### input series of total precipitation [mm/h]
InputsFracSolidPrecip = InputsModel$LayerFracSolidPrecip[[iLayer]][IndPeriod1], ### input series of fraction of solid precipitation [0-1]
InputsTemp = InputsModel$LayerTemp[[iLayer]][IndPeriod1], ### input series of air mean temperature [degC]
MeanAnSolidPrecip = RunOptions$MeanAnSolidPrecip[iLayer], ### value of annual mean solid precip [mm/y]
NParam = as.integer(NParamCN), ### number of model parameters = 2 or 4
Param = as.double(ParamCemaNeige), ### parameter set
NStates = as.integer(NStates), ### number of state variables used for model initialising = 4
StateStart = StateStartCemaNeige, ### state variables used when the model run starts
IsHyst = as.integer(IsHyst), ### use of hysteresis
NOutputs = as.integer(length(IndOutputsCemaNeige)), ### number of output series
IndOutputs = IndOutputsCemaNeige, ### indices of output series
## outputs
Outputs = matrix(as.double(-99e9), nrow = LInputSeries, ncol = length(IndOutputsCemaNeige)), ### output series [mm, mm/h or degC]
StateEnd = rep(as.double(-99e9), as.integer(NStates)) ### state variables at the end of the model run
)
RESULTS$Outputs[RESULTS$Outputs <= -99e8] <- NA
RESULTS$StateEnd[RESULTS$StateEnd <= -99e8] <- NA
## Data storage
CemaNeigeLayers[[iLayer]] <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i])
names(CemaNeigeLayers[[iLayer]]) <- RunOptions$FortranOutputs$CN[IndOutputsCemaNeige]
IndPliqAndMelt <- which(names(CemaNeigeLayers[[iLayer]]) == "PliqAndMelt")
if (iLayer == 1) {
CatchMeltAndPliq <- RESULTS$Outputs[, IndPliqAndMelt] / NLayers
}
if (iLayer > 1) {
CatchMeltAndPliq <- CatchMeltAndPliq + RESULTS$Outputs[, IndPliqAndMelt] / NLayers
}
if (ExportStateEnd) {
CemaNeigeStateEnd <- c(CemaNeigeStateEnd, RESULTS$StateEnd)
}
rm(RESULTS)
} ### ENDFOR iLayer
names(CemaNeigeLayers) <- sprintf("Layer%02i", seq_len(NLayers))
} ### ENDIF RunSnowModule
if (!inherits(RunOptions, "CemaNeige")) {
CemaNeigeLayers <- list()
CemaNeigeStateEnd <- NULL
NameCemaNeigeLayers <- NULL
CatchMeltAndPliq <- InputsModel$Precip[IndPeriod1]
}
## GR model
if ("all" %in% RunOptions$Outputs_Sim) {
IndOutputsMod <- as.integer(1:length(RunOptions$FortranOutputs$GR))
} else {
IndOutputsMod <- which(RunOptions$FortranOutputs$GR %in% RunOptions$Outputs_Sim)
}
## Use of IniResLevels
if (!is.null(RunOptions$IniResLevels)) {
RunOptions$IniStates[1] <- RunOptions$IniResLevels[1] * ParamMod[1] ### production store level (mm)
RunOptions$IniStates[2] <- RunOptions$IniResLevels[2] * ParamMod[3] ### routing store level (mm)
if (IsIntStore) {
RunOptions$IniStates[4] <- RunOptions$IniResLevels[4] * Imax ### interception store level (mm)
}
}
## Call GR model Fortan
RESULTS <- .Fortran("frun_gr5h", PACKAGE = "airGR",
## inputs
LInputs = LInputSeries, ### length of input and output series
InputsPrecip = CatchMeltAndPliq, ### input series of total precipitation [mm/h]
InputsPE = InputsModel$PotEvap[IndPeriod1], ### input series potential evapotranspiration [mm/h]
NParam = NParamMod, ### number of model parameter
Param = ParamMod, ### parameter set
NStates = NStatesMod, ### number of state variables used for model initialising
StateStart = RunOptions$IniStates[1:NStatesMod], ### state variables used when the model run starts
Imax = Imax, ### maximal capacity of interception store
NOutputs = as.integer(length(IndOutputsMod)), ### number of output series
IndOutputs = IndOutputsMod, ### indices of output series
## outputs
Outputs = matrix(as.double(-99e9), nrow = LInputSeries, ncol = length(IndOutputsMod)), ### output series [mm or mm/h]
StateEnd = rep(as.double(-99e9), NStatesMod) ### state variables at the end of the model run
)
RESULTS$Outputs[RESULTS$Outputs <= -99e8] <- NA
RESULTS$StateEnd[RESULTS$StateEnd <= -99e8] <- NA
if (ExportStateEnd) {
RESULTS$StateEnd[-3L] <- ifelse(RESULTS$StateEnd[-3L] < 0, 0, RESULTS$StateEnd[-3L]) ### remove negative values except for the ExpStore location
idNStates <- seq_len(NStates*NLayers) %% NStates
RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_CemaNeigeGR5H, InputsModel = InputsModel, IsHyst = IsHyst,
ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL,
IntStore = RESULTS$StateEnd[4L],
UH1 = NULL,
UH2 = RESULTS$StateEnd[(1:(40*24)) + (7+20*24)],
GCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 1]],
eTGCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 2]],
GthrCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 3]],
GlocmaxCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 0]],
verbose = FALSE)
}
if (inherits(RunOptions, "CemaNeige") & "Precip" %in% RunOptions$Outputs_Sim) {
RESULTS$Outputs[, which(RunOptions$FortranOutputs$GR[IndOutputsMod] == "Precip")] <-
InputsModel$Precip[IndPeriod1]
}
## OutputsModel generation
.GetOutputsModelGR(InputsModel,
RunOptions,
RESULTS,
LInputSeries,
Param,
CemaNeigeLayers)
}
R/RunModel_CemaNeigeGR5J.R
View file @ 4bccfee2
RunModel_CemaNeigeGR5J <- function(InputsModel,RunOptions,Param){ RunModel_CemaNeigeGR5J <- function(InputsModel, RunOptions, Param) {
NParam <- 7;
FortranOutputsCemaNeige <- c("Pliq","Psol","SnowPack","ThermalState","Gratio","PotMelt","Melt","PliqAndMelt", "Temp"); ## Initialization of variables
FortranOutputsMod <- c("PotEvap","Precip","Prod","AE","Perc","PR","Q9","Q1","Rout","Exch","AExch","QR","QD","Qsim"); IsHyst <- inherits(RunOptions, "hysteresis")
NParamCN <- RunOptions$FeatFUN_MOD$NbParam - 5L
##Arguments_check NStates <- 4L
if(inherits(InputsModel,"InputsModel")==FALSE){ stop("InputsModel must be of class 'InputsModel' \n"); return(NULL); }
if(inherits(InputsModel,"daily" )==FALSE){ stop("InputsModel must be of class 'daily' \n"); return(NULL); }
if(inherits(InputsModel,"GR" )==FALSE){ stop("InputsModel must be of class 'GR' \n"); return(NULL); } .ArgumentsCheckGR(InputsModel, RunOptions, Param)
if(inherits(InputsModel,"CemaNeige" )==FALSE){ stop("InputsModel must be of class 'CemaNeige' \n"); return(NULL); }
if(inherits(RunOptions,"RunOptions" )==FALSE){ stop("RunOptions must be of class 'RunOptions' \n"); return(NULL); } Param <- as.double(Param)
if(inherits(RunOptions,"GR" )==FALSE){ stop("RunOptions must be of class 'GR' \n"); return(NULL); }
if(inherits(RunOptions,"CemaNeige" )==FALSE){ stop("RunOptions must be of class 'CemaNeige' \n"); return(NULL); }
if(!is.vector(Param)){ stop("Param must be a vector \n"); return(NULL); } Param_X1X3_threshold <- 1e-2
if(sum(!is.na(Param))!=NParam){ stop(paste("Param must be a vector of length ",NParam," and contain no NA \n",sep="")); return(NULL); } Param_X4_threshold <- 0.5
Param <- as.double(Param); if (Param[1L] < Param_X1X3_threshold) {
warning(sprintf("Param[1] (X1: production store capacity [mm]) < %.2f\n X1 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
Param_X1X3_threshold <- 1e-2 Param[1L] <- Param_X1X3_threshold
if (Param[1L] < Param_X1X3_threshold) { }
warning(sprintf("Param[1] (X1: production store capacity [mm]) < %.2f\n X1 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold)) if (Param[3L] < Param_X1X3_threshold) {
Param[1L] <- Param_X1X3_threshold warning(sprintf("Param[3] (X3: routing store capacity [mm]) < %.2f\n X3 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
Param[3L] <- Param_X1X3_threshold
}
if (Param[4L] < Param_X4_threshold) {
warning(sprintf("Param[4] (X4: unit hydrograph time constant [d]) < %.2f\n X4 set to %.2f", Param_X4_threshold, Param_X4_threshold))
Param[4L] <- Param_X4_threshold
}
## Input data preparation
if (identical(RunOptions$IndPeriod_WarmUp, 0L)) {
RunOptions$IndPeriod_WarmUp <- NULL
}
IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp, RunOptions$IndPeriod_Run)
LInputSeries <- as.integer(length(IndPeriod1))
IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries
ParamCemaNeige <- Param[(length(Param) - 1 - 2 * as.integer(IsHyst)):length(Param)]
NParamMod <- as.integer(length(Param) - (2 + 2 * as.integer(IsHyst)))
ParamMod <- Param[1:NParamMod]
NLayers <- length(InputsModel$LayerPrecip)
NStatesMod <- as.integer(length(RunOptions$IniStates) - NStates * NLayers)
ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim
ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim
## CemaNeige________________________________________________________________________________
if (inherits(RunOptions, "CemaNeige")) {
if ("all" %in% RunOptions$Outputs_Sim) {
IndOutputsCemaNeige <- as.integer(1:length(RunOptions$FortranOutputs$CN))
} else {
IndOutputsCemaNeige <- which(RunOptions$FortranOutputs$CN %in% RunOptions$Outputs_Sim)
}
CemaNeigeLayers <- list()
CemaNeigeStateEnd <- NULL
NameCemaNeigeLayers <- "CemaNeigeLayers"
## Call CemaNeige Fortran_________________________
for (iLayer in 1:NLayers) {
if (!IsHyst) {
StateStartCemaNeige <- RunOptions$IniStates[(7 + 20 + 40) + c(iLayer, iLayer+NLayers)]
} else {
StateStartCemaNeige <- RunOptions$IniStates[(7 + 20 + 40) + c(iLayer, iLayer+NLayers, iLayer+2*NLayers, iLayer+3*NLayers)]
} }
if (Param[3L] < Param_X1X3_threshold) { RESULTS <- .Fortran("frun_cemaneige", PACKAGE = "airGR",
warning(sprintf("Param[3] (X3: routing store capacity [mm]) < %.2f\n X3 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold)) ## inputs
Param[3L] <- Param_X1X3_threshold LInputs = LInputSeries, ### length of input and output series
InputsPrecip = InputsModel$LayerPrecip[[iLayer]][IndPeriod1], ### input series of total precipitation [mm/d]
InputsFracSolidPrecip = InputsModel$LayerFracSolidPrecip[[iLayer]][IndPeriod1], ### input series of fraction of solid precipitation [0-1]
InputsTemp = InputsModel$LayerTemp[[iLayer]][IndPeriod1], ### input series of air mean temperature [degC]
MeanAnSolidPrecip = RunOptions$MeanAnSolidPrecip[iLayer], ### value of annual mean solid precip [mm/y]
NParam = as.integer(NParamCN), ### number of model parameters = 2 or 4
Param = as.double(ParamCemaNeige), ### parameter set
NStates = as.integer(NStates), ### number of state variables used for model initialising = 4
StateStart = StateStartCemaNeige, ### state variables used when the model run starts
IsHyst = as.integer(IsHyst), ### use of hysteresis
NOutputs = as.integer(length(IndOutputsCemaNeige)), ### number of output series
IndOutputs = IndOutputsCemaNeige, ### indices of output series
## outputs
Outputs = matrix(as.double(-99e9), nrow = LInputSeries, ncol = length(IndOutputsCemaNeige)), ### output series [mm, mm/d or degC]
StateEnd = rep(as.double(-99e9), as.integer(NStates)) ### state variables at the end of the model run
)
RESULTS$Outputs[RESULTS$Outputs <= -99e8] <- NA
RESULTS$StateEnd[RESULTS$StateEnd <= -99e8] <- NA
## Data storage
CemaNeigeLayers[[iLayer]] <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i])
names(CemaNeigeLayers[[iLayer]]) <- RunOptions$FortranOutputs$CN[IndOutputsCemaNeige]
IndPliqAndMelt <- which(names(CemaNeigeLayers[[iLayer]]) == "PliqAndMelt")
if (iLayer == 1) {
CatchMeltAndPliq <- RESULTS$Outputs[, IndPliqAndMelt] / NLayers
} }
if (iLayer > 1) {
##Input_data_preparation CatchMeltAndPliq <- CatchMeltAndPliq + RESULTS$Outputs[, IndPliqAndMelt] / NLayers
if(identical(RunOptions$IndPeriod_WarmUp,as.integer(0))){ RunOptions$IndPeriod_WarmUp <- NULL; }
IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp,RunOptions$IndPeriod_Run);
LInputSeries <- as.integer(length(IndPeriod1))
IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries;
ParamCemaNeige <- Param[(length(Param)-1):length(Param)];
NParamMod <- as.integer(length(Param)-2);
ParamMod <- Param[1:NParamMod];
NLayers <- length(InputsModel$LayerPrecip);
NStatesMod <- as.integer(length(RunOptions$IniStates)-2*NLayers);
ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim;
ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim;
##SNOW_MODULE________________________________________________________________________________##
if(RunOptions$RunSnowModule==TRUE){
if("all" %in% RunOptions$Outputs_Sim){ IndOutputsCemaNeige <- as.integer(1:length(FortranOutputsCemaNeige));
} else { IndOutputsCemaNeige <- which(FortranOutputsCemaNeige %in% RunOptions$Outputs_Sim); }
CemaNeigeLayers <- list(); CemaNeigeStateEnd <- NULL; NameCemaNeigeLayers <- "CemaNeigeLayers";
##Call_DLL_CemaNeige_________________________
for(iLayer in 1:NLayers){
StateStartCemaNeige <- RunOptions$IniStates[ (NStatesMod+2*(iLayer-1)+1):(NStatesMod+2*(iLayer-1)+2) ];
RESULTS <- .Fortran("frun_CemaNeige",PACKAGE="airGR",
##inputs
LInputs=LInputSeries, ### length of input and output series
InputsPrecip=InputsModel$LayerPrecip[[iLayer]][IndPeriod1], ### input series of total precipitation [mm/d]
InputsFracSolidPrecip=InputsModel$LayerFracSolidPrecip[[iLayer]][IndPeriod1], ### input series of fraction of solid precipitation [0-1]
InputsTemp=InputsModel$LayerTemp[[iLayer]][IndPeriod1], ### input series of air mean temperature [degC]
MeanAnSolidPrecip=RunOptions$MeanAnSolidPrecip[iLayer], ### value of annual mean solid precip [mm/y]
NParam=as.integer(2), ### number of model parameter = 2
Param=ParamCemaNeige, ### parameter set
NStates=as.integer(2), ### number of state variables used for model initialising = 2
StateStart=StateStartCemaNeige, ### state variables used when the model run starts
NOutputs=as.integer(length(IndOutputsCemaNeige)), ### number of output series
IndOutputs=IndOutputsCemaNeige, ### indices of output series
##outputs
Outputs=matrix(as.double(-999.999),nrow=LInputSeries,ncol=length(IndOutputsCemaNeige)), ### output series [mm]
StateEnd=rep(as.double(-999.999),as.integer(2)) ### state variables at the end of the model run (reservoir levels [mm] and HU)
)
RESULTS$Outputs[ round(RESULTS$Outputs ,3)==(-999.999)] <- NA;
RESULTS$StateEnd[round(RESULTS$StateEnd,3)==(-999.999)] <- NA;
##Data_storage
CemaNeigeLayers[[iLayer]] <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]);
names(CemaNeigeLayers[[iLayer]]) <- FortranOutputsCemaNeige[IndOutputsCemaNeige];
IndPliqAndMelt <- which(names(CemaNeigeLayers[[iLayer]]) == "PliqAndMelt");
if(iLayer==1){ CatchMeltAndPliq <- RESULTS$Outputs[,IndPliqAndMelt]/NLayers; }
if(iLayer >1){ CatchMeltAndPliq <- CatchMeltAndPliq + RESULTS$Outputs[,IndPliqAndMelt]/NLayers; }
if(ExportStateEnd){ CemaNeigeStateEnd <- c(CemaNeigeStateEnd,RESULTS$StateEnd); }
rm(RESULTS);
} ###ENDFOR_iLayer
names(CemaNeigeLayers) <- paste("Layer",formatC(1:NLayers,width=2,flag="0"),sep="");
} ###ENDIF_RunSnowModule
if(RunOptions$RunSnowModule==FALSE){
CemaNeigeLayers <- list(); CemaNeigeStateEnd <- NULL; NameCemaNeigeLayers <- NULL;
CatchMeltAndPliq <- InputsModel$Precip[IndPeriod1]; }
##MODEL______________________________________________________________________________________##
if("all" %in% RunOptions$Outputs_Sim){ IndOutputsMod <- as.integer(1:length(FortranOutputsMod));
} else { IndOutputsMod <- which(FortranOutputsMod %in% RunOptions$Outputs_Sim); }
##Use_of_IniResLevels
if("IniResLevels" %in% names(RunOptions)){
RunOptions$IniStates[1] <- RunOptions$IniResLevels[1]*ParamMod[1]; ### production store level (mm)
RunOptions$IniStates[2] <- RunOptions$IniResLevels[2]*ParamMod[3]; ### routing store level (mm)
} }
if (ExportStateEnd) {
##Call_fortan CemaNeigeStateEnd <- c(CemaNeigeStateEnd, RESULTS$StateEnd)
RESULTS <- .Fortran("frun_GR5J",PACKAGE="airGR", }
##inputs rm(RESULTS)
LInputs=LInputSeries, ### length of input and output series } ### ENDFOR iLayer
InputsPrecip=CatchMeltAndPliq, ### input series of total precipitation [mm/d] names(CemaNeigeLayers) <- sprintf("Layer%02i", seq_len(NLayers))
InputsPE=InputsModel$PotEvap[IndPeriod1], ### input series potential evapotranspiration [mm/d] } ### ENDIF RunSnowModule
NParam=NParamMod, ### number of model parameter if (!inherits(RunOptions, "CemaNeige")) {
Param=ParamMod, ### parameter set CemaNeigeLayers <- list()
NStates=NStatesMod, ### number of state variables used for model initialising CemaNeigeStateEnd <- NULL
StateStart=RunOptions$IniStates[1:NStatesMod], ### state variables used when the model run starts NameCemaNeigeLayers <- NULL
NOutputs=as.integer(length(IndOutputsMod)), ### number of output series CatchMeltAndPliq <- InputsModel$Precip[IndPeriod1]
IndOutputs=IndOutputsMod, ### indices of output series }
##outputs
Outputs=matrix(as.double(-999.999),nrow=LInputSeries,ncol=length(IndOutputsMod)), ### output series [mm]
StateEnd=rep(as.double(-999.999),NStatesMod) ### state variables at the end of the model run
) ## GR model______________________________________________________________________________________
RESULTS$Outputs[ round(RESULTS$Outputs ,3)==(-999.999)] <- NA; if ("all" %in% RunOptions$Outputs_Sim) {
RESULTS$StateEnd[round(RESULTS$StateEnd,3)==(-999.999)] <- NA; IndOutputsMod <- as.integer(1:length(RunOptions$FortranOutputs$GR))
if(RunOptions$RunSnowModule & "Precip" %in% RunOptions$Outputs_Sim){ RESULTS$Outputs[,which(FortranOutputsMod[IndOutputsMod]=="Precip")] <- InputsModel$Precip[IndPeriod1]; } } else {
IndOutputsMod <- which(RunOptions$FortranOutputs$GR %in% RunOptions$Outputs_Sim)
##Output_data_preparation }
##OutputsModel_only
if(ExportDatesR==FALSE & ExportStateEnd==FALSE){ ## Use of IniResLevels
OutputsModel <- c( lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]), if (!is.null(RunOptions$IniResLevels)) {
list(CemaNeigeLayers) ); RunOptions$IniStates[1] <- RunOptions$IniResLevels[1] * ParamMod[1] ### production store level (mm)
names(OutputsModel) <- c(FortranOutputsMod[IndOutputsMod],NameCemaNeigeLayers); } RunOptions$IniStates[2] <- RunOptions$IniResLevels[2] * ParamMod[3] ### routing store level (mm)
##DatesR_and_OutputsModel_only }
if(ExportDatesR==TRUE & ExportStateEnd==FALSE){
OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]), ## Call GR model Fortan
lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]), RESULTS <- .Fortran("frun_gr5j", PACKAGE = "airGR",
list(CemaNeigeLayers) ); ## inputs
names(OutputsModel) <- c("DatesR",FortranOutputsMod[IndOutputsMod],NameCemaNeigeLayers); } LInputs = LInputSeries, ### length of input and output series
##OutputsModel_and_SateEnd_only InputsPrecip = CatchMeltAndPliq, ### input series of total precipitation [mm/d]
if(ExportDatesR==FALSE & ExportStateEnd==TRUE){ InputsPE = InputsModel$PotEvap[IndPeriod1], ### input series potential evapotranspiration [mm/d]
OutputsModel <- c( lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]), NParam = NParamMod, ### number of model parameter
list(CemaNeigeLayers), Param = ParamMod, ### parameter set
list(c(RESULTS$StateEnd,CemaNeigeStateEnd)) ); NStates = NStatesMod, ### number of state variables used for model initialising
names(OutputsModel) <- c(FortranOutputsMod[IndOutputsMod],NameCemaNeigeLayers,"StateEnd"); } StateStart = RunOptions$IniStates[1:NStatesMod], ### state variables used when the model run starts
##DatesR_and_OutputsModel_and_SateEnd NOutputs = as.integer(length(IndOutputsMod)), ### number of output series
if(ExportDatesR==TRUE & ExportStateEnd==TRUE){ IndOutputs = IndOutputsMod, ### indices of output series
OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]), ## outputs
lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]), Outputs = matrix(as.double(-99e9), nrow = LInputSeries, ncol = length(IndOutputsMod)), ### output series [mm or mm/d]
list(CemaNeigeLayers), StateEnd = rep(as.double(-99e9), NStatesMod) ### state variables at the end of the model run
list(c(RESULTS$StateEnd,CemaNeigeStateEnd)) ); )
names(OutputsModel) <- c("DatesR",FortranOutputsMod[IndOutputsMod],NameCemaNeigeLayers,"StateEnd"); } RESULTS$Outputs[RESULTS$Outputs <= -99e8] <- NA
RESULTS$StateEnd[RESULTS$StateEnd <= -99e8] <- NA
##End if (ExportStateEnd) {
rm(RESULTS); RESULTS$StateEnd[-3L] <- ifelse(RESULTS$StateEnd[-3L] < 0, 0, RESULTS$StateEnd[-3L]) ### remove negative values except for the ExpStore location
class(OutputsModel) <- c("OutputsModel","daily","GR","CemaNeige"); idNStates <- seq_len(NStates*NLayers) %% NStates
return(OutputsModel); RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_CemaNeigeGR5J, InputsModel = InputsModel, IsHyst = IsHyst,
ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL,
UH1 = NULL,
UH2 = RESULTS$StateEnd[(1:40) + (7+20)],
GCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 1]],
eTGCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 2]],
GthrCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 3]],
GlocmaxCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 0]],
verbose = FALSE)
}
if (inherits(RunOptions, "CemaNeige") & "Precip" %in% RunOptions$Outputs_Sim) {
RESULTS$Outputs[, which(RunOptions$FortranOutputs$GR[IndOutputsMod] == "Precip")] <-
InputsModel$Precip[IndPeriod1]
}
## OutputsModel generation
.GetOutputsModelGR(InputsModel,
RunOptions,
RESULTS,
LInputSeries,
Param,
CemaNeigeLayers)
} }
R/RunModel_CemaNeigeGR6J.R
View file @ 4bccfee2
RunModel_CemaNeigeGR6J <- function(InputsModel,RunOptions,Param){ RunModel_CemaNeigeGR6J <- function(InputsModel, RunOptions, Param) {
NParam <- 8;
FortranOutputsCemaNeige <- c("Pliq","Psol","SnowPack","ThermalState","Gratio","PotMelt","Melt","PliqAndMelt", "Temp"); ## Initialization of variables
FortranOutputsMod <- c("PotEvap","Precip","Prod","AE","Perc","PR","Q9","Q1","Rout","Exch","AExch","QR","QR1","Exp","QD","Qsim"); IsHyst <- inherits(RunOptions, "hysteresis")
NParamCN <- RunOptions$FeatFUN_MOD$NbParam - 6L
##Arguments_check NStates <- 4L
if(inherits(InputsModel,"InputsModel")==FALSE){ stop("InputsModel must be of class 'InputsModel' \n"); return(NULL); }
if(inherits(InputsModel,"daily" )==FALSE){ stop("InputsModel must be of class 'daily' \n"); return(NULL); }
if(inherits(InputsModel,"GR" )==FALSE){ stop("InputsModel must be of class 'GR' \n"); return(NULL); } .ArgumentsCheckGR(InputsModel, RunOptions, Param)
if(inherits(InputsModel,"CemaNeige" )==FALSE){ stop("InputsModel must be of class 'CemaNeige' \n"); return(NULL); }
if(inherits(RunOptions,"RunOptions" )==FALSE){ stop("RunOptions must be of class 'RunOptions' \n"); return(NULL); } Param <- as.double(Param)
if(inherits(RunOptions,"GR" )==FALSE){ stop("RunOptions must be of class 'GR' \n"); return(NULL); }
if(inherits(RunOptions,"CemaNeige" )==FALSE){ stop("RunOptions must be of class 'CemaNeige' \n"); return(NULL); }
if(!is.vector(Param)){ stop("Param must be a vector \n"); return(NULL); } Param_X1X3X6_threshold <- 1e-2
if(sum(!is.na(Param))!=NParam){ stop(paste("Param must be a vector of length ",NParam," and contain no NA \n",sep="")); return(NULL); } Param_X4_threshold <- 0.5
Param <- as.double(Param); if (Param[1L] < Param_X1X3X6_threshold) {
warning(sprintf("Param[1] (X1: production store capacity [mm]) < %.2f\n X1 set to %.2f", Param_X1X3X6_threshold, Param_X1X3X6_threshold))
Param_X1X3X6_threshold <- 1e-2 Param[1L] <- Param_X1X3X6_threshold
if (Param[1L] < Param_X1X3X6_threshold) { }
warning(sprintf("Param[1] (X1: production store capacity [mm]) < %.2f\n X1 set to %.2f", Param_X1X3X6_threshold, Param_X1X3X6_threshold)) if (Param[3L] < Param_X1X3X6_threshold) {
Param[1L] <- Param_X1X3X6_threshold warning(sprintf("Param[3] (X3: routing store capacity [mm]) < %.2f\n X3 set to %.2f", Param_X1X3X6_threshold, Param_X1X3X6_threshold))
Param[3L] <- Param_X1X3X6_threshold
}
if (Param[4L] < Param_X4_threshold) {
warning(sprintf("Param[4] (X4: unit hydrograph time constant [d]) < %.2f\n X4 set to %.2f", Param_X4_threshold, Param_X4_threshold))
Param[4L] <- Param_X4_threshold
}
if (Param[6L] < Param_X1X3X6_threshold) {
warning(sprintf("Param[6] (X6: coefficient for emptying exponential store [mm]) < %.2f\n X6 set to %.2f", Param_X1X3X6_threshold, Param_X1X3X6_threshold))
Param[6L] <- Param_X1X3X6_threshold
}
## Input data preparation
if (identical(RunOptions$IndPeriod_WarmUp, 0L)) {
RunOptions$IndPeriod_WarmUp <- NULL
}
IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp, RunOptions$IndPeriod_Run)
LInputSeries <- as.integer(length(IndPeriod1))
IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries
ParamCemaNeige <- Param[(length(Param) - 1 - 2 * as.integer(IsHyst)):length(Param)]
NParamMod <- as.integer(length(Param) - (2 + 2 * as.integer(IsHyst)))
ParamMod <- Param[1:NParamMod]
NLayers <- length(InputsModel$LayerPrecip)
NStatesMod <- as.integer(length(RunOptions$IniStates) - NStates * NLayers)
ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim
ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim
## CemaNeige________________________________________________________________________________
if (inherits(RunOptions, "CemaNeige")) {
if ("all" %in% RunOptions$Outputs_Sim) {
IndOutputsCemaNeige <- as.integer(1:length(RunOptions$FortranOutputs$CN))
} else {
IndOutputsCemaNeige <- which(RunOptions$FortranOutputs$CN %in% RunOptions$Outputs_Sim)
}
CemaNeigeLayers <- list()
CemaNeigeStateEnd <- NULL
NameCemaNeigeLayers <- "CemaNeigeLayers"
## Call CemaNeige Fortran_________________________
for (iLayer in 1:NLayers) {
if (!IsHyst) {
StateStartCemaNeige <- RunOptions$IniStates[(7 + 20 + 40) + c(iLayer, iLayer+NLayers)]
} else {
StateStartCemaNeige <- RunOptions$IniStates[(7 + 20 + 40) + c(iLayer, iLayer+NLayers, iLayer+2*NLayers, iLayer+3*NLayers)]
} }
if (Param[3L] < Param_X1X3X6_threshold) { RESULTS <- .Fortran("frun_cemaneige", PACKAGE = "airGR",
warning(sprintf("Param[3] (X3: routing store capacity [mm]) < %.2f\n X3 set to %.2f", Param_X1X3X6_threshold, Param_X1X3X6_threshold)) ## inputs
Param[3L] <- Param_X1X3X6_threshold LInputs = LInputSeries, ### length of input and output series
InputsPrecip = InputsModel$LayerPrecip[[iLayer]][IndPeriod1], ### input series of total precipitation [mm/d]
InputsFracSolidPrecip = InputsModel$LayerFracSolidPrecip[[iLayer]][IndPeriod1], ### input series of fraction of solid precipitation [0-1]
InputsTemp = InputsModel$LayerTemp[[iLayer]][IndPeriod1], ### input series of air mean temperature [degC]
MeanAnSolidPrecip = RunOptions$MeanAnSolidPrecip[iLayer], ### value of annual mean solid precip [mm/y]
NParam = as.integer(NParamCN), ### number of model parameters = 2 or 4
Param = as.double(ParamCemaNeige), ### parameter set
NStates = as.integer(NStates), ### number of state variables used for model initialising = 4
StateStart = StateStartCemaNeige, ### state variables used when the model run starts
IsHyst = as.integer(IsHyst), ### use of hysteresis
NOutputs = as.integer(length(IndOutputsCemaNeige)), ### number of output series
IndOutputs = IndOutputsCemaNeige, ### indices of output series
## outputs
Outputs = matrix(as.double(-99e9), nrow = LInputSeries, ncol = length(IndOutputsCemaNeige)), ### output series [mm, mm/d or degC]
StateEnd = rep(as.double(-99e9), as.integer(NStates)) ### state variables at the end of the model run
)
RESULTS$Outputs[RESULTS$Outputs <= -99e8] <- NA
RESULTS$StateEnd[RESULTS$StateEnd <= -99e8] <- NA
## Data storage
CemaNeigeLayers[[iLayer]] <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i])
names(CemaNeigeLayers[[iLayer]]) <- RunOptions$FortranOutputs$CN[IndOutputsCemaNeige]
IndPliqAndMelt <- which(names(CemaNeigeLayers[[iLayer]]) == "PliqAndMelt")
if (iLayer == 1) {
CatchMeltAndPliq <- RESULTS$Outputs[, IndPliqAndMelt] / NLayers
} }
if (Param[6L] < Param_X1X3X6_threshold) { if (iLayer > 1) {
warning(sprintf("Param[6] (X6: coefficient for emptying exponential store [mm]) < %.2f\n X6 set to %.2f", Param_X1X3X6_threshold, Param_X1X3X6_threshold)) CatchMeltAndPliq <- CatchMeltAndPliq + RESULTS$Outputs[, IndPliqAndMelt] / NLayers
Param[6L] <- Param_X1X3X6_threshold
}
##Input_data_preparation
if(identical(RunOptions$IndPeriod_WarmUp,as.integer(0))){ RunOptions$IndPeriod_WarmUp <- NULL; }
IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp,RunOptions$IndPeriod_Run);
LInputSeries <- as.integer(length(IndPeriod1))
IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries;
ParamCemaNeige <- Param[(length(Param)-1):length(Param)];
NParamMod <- as.integer(length(Param)-2);
ParamMod <- Param[1:NParamMod];
NLayers <- length(InputsModel$LayerPrecip);
NStatesMod <- as.integer(length(RunOptions$IniStates)-2*NLayers);
ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim;
ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim;
##SNOW_MODULE________________________________________________________________________________##
if(RunOptions$RunSnowModule==TRUE){
if("all" %in% RunOptions$Outputs_Sim){ IndOutputsCemaNeige <- as.integer(1:length(FortranOutputsCemaNeige));
} else { IndOutputsCemaNeige <- which(FortranOutputsCemaNeige %in% RunOptions$Outputs_Sim); }
CemaNeigeLayers <- list(); CemaNeigeStateEnd <- NULL; NameCemaNeigeLayers <- "CemaNeigeLayers";
##Call_DLL_CemaNeige_________________________
for(iLayer in 1:NLayers){
StateStartCemaNeige <- RunOptions$IniStates[ (NStatesMod+2*(iLayer-1)+1):(NStatesMod+2*(iLayer-1)+2) ];
RESULTS <- .Fortran("frun_CemaNeige",PACKAGE="airGR",
##inputs
LInputs=LInputSeries, ### length of input and output series
InputsPrecip=InputsModel$LayerPrecip[[iLayer]][IndPeriod1], ### input series of total precipitation [mm/d]
InputsFracSolidPrecip=InputsModel$LayerFracSolidPrecip[[iLayer]][IndPeriod1], ### input series of fraction of solid precipitation [0-1]
InputsTemp=InputsModel$LayerTemp[[iLayer]][IndPeriod1], ### input series of air mean temperature [degC]
MeanAnSolidPrecip=RunOptions$MeanAnSolidPrecip[iLayer], ### value of annual mean solid precip [mm/y]
NParam=as.integer(2), ### number of model parameter = 2
Param=ParamCemaNeige, ### parameter set
NStates=as.integer(2), ### number of state variables used for model initialising = 2
StateStart=StateStartCemaNeige, ### state variables used when the model run starts
NOutputs=as.integer(length(IndOutputsCemaNeige)), ### number of output series
IndOutputs=IndOutputsCemaNeige, ### indices of output series
##outputs
Outputs=matrix(as.double(-999.999),nrow=LInputSeries,ncol=length(IndOutputsCemaNeige)), ### output series [mm]
StateEnd=rep(as.double(-999.999),as.integer(2)) ### state variables at the end of the model run (reservoir levels [mm] and HU)
)
RESULTS$Outputs[ round(RESULTS$Outputs ,3)==(-999.999)] <- NA;
RESULTS$StateEnd[round(RESULTS$StateEnd,3)==(-999.999)] <- NA;
##Data_storage
CemaNeigeLayers[[iLayer]] <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]);
names(CemaNeigeLayers[[iLayer]]) <- FortranOutputsCemaNeige[IndOutputsCemaNeige];
IndPliqAndMelt <- which(names(CemaNeigeLayers[[iLayer]]) == "PliqAndMelt");
if(iLayer==1){ CatchMeltAndPliq <- RESULTS$Outputs[,IndPliqAndMelt]/NLayers; }
if(iLayer >1){ CatchMeltAndPliq <- CatchMeltAndPliq + RESULTS$Outputs[,IndPliqAndMelt]/NLayers; }
if(ExportStateEnd){ CemaNeigeStateEnd <- c(CemaNeigeStateEnd,RESULTS$StateEnd); }
rm(RESULTS);
} ###ENDFOR_iLayer
names(CemaNeigeLayers) <- paste("Layer",formatC(1:NLayers,width=2,flag="0"),sep="");
} ###ENDIF_RunSnowModule
if(RunOptions$RunSnowModule==FALSE){
CemaNeigeLayers <- list(); CemaNeigeStateEnd <- NULL; NameCemaNeigeLayers <- NULL;
CatchMeltAndPliq <- InputsModel$Precip[IndPeriod1]; }
##MODEL______________________________________________________________________________________##
if("all" %in% RunOptions$Outputs_Sim){ IndOutputsMod <- as.integer(1:length(FortranOutputsMod));
} else { IndOutputsMod <- which(FortranOutputsMod %in% RunOptions$Outputs_Sim); }
##Use_of_IniResLevels
if("IniResLevels" %in% names(RunOptions)){
RunOptions$IniStates[1] <- RunOptions$IniResLevels[1]*ParamMod[1]; ### production store level (mm)
RunOptions$IniStates[2] <- RunOptions$IniResLevels[2]*ParamMod[3]; ### routing store level (mm)
} }
if (ExportStateEnd) {
##Call_fortan CemaNeigeStateEnd <- c(CemaNeigeStateEnd, RESULTS$StateEnd)
RESULTS <- .Fortran("frun_GR6J",PACKAGE="airGR", }
##inputs rm(RESULTS)
LInputs=LInputSeries, ### length of input and output series } ### ENDFOR iLayer
InputsPrecip=CatchMeltAndPliq, ### input series of total precipitation [mm/d] names(CemaNeigeLayers) <- sprintf("Layer%02i", seq_len(NLayers))
InputsPE=InputsModel$PotEvap[IndPeriod1], ### input series potential evapotranspiration [mm/d] } ### ENDIF RunSnowModule
NParam=NParamMod, ### number of model parameter if (!inherits(RunOptions, "CemaNeige")) {
Param=ParamMod, ### parameter set CemaNeigeLayers <- list()
NStates=NStatesMod, ### number of state variables used for model initialising CemaNeigeStateEnd <- NULL
StateStart=RunOptions$IniStates[1:NStatesMod], ### state variables used when the model run starts NameCemaNeigeLayers <- NULL
NOutputs=as.integer(length(IndOutputsMod)), ### number of output series CatchMeltAndPliq <- InputsModel$Precip[IndPeriod1]
IndOutputs=IndOutputsMod, ### indices of output series }
##outputs
Outputs=matrix(as.double(-999.999),nrow=LInputSeries,ncol=length(IndOutputsMod)), ### output series [mm]
StateEnd=rep(as.double(-999.999),NStatesMod) ### state variables at the end of the model run
) ## GR model______________________________________________________________________________________
RESULTS$Outputs[ round(RESULTS$Outputs ,3)==(-999.999)] <- NA; if ("all" %in% RunOptions$Outputs_Sim) {
RESULTS$StateEnd[round(RESULTS$StateEnd,3)==(-999.999)] <- NA; IndOutputsMod <- as.integer(1:length(RunOptions$FortranOutputs$GR))
if(RunOptions$RunSnowModule & "Precip" %in% RunOptions$Outputs_Sim){ RESULTS$Outputs[,which(FortranOutputsMod[IndOutputsMod]=="Precip")] <- InputsModel$Precip[IndPeriod1]; } } else {
IndOutputsMod <- which(RunOptions$FortranOutputs$GR %in% RunOptions$Outputs_Sim)
##Output_data_preparation }
##OutputsModel_only
if(ExportDatesR==FALSE & ExportStateEnd==FALSE){ ## Use of IniResLevels
OutputsModel <- c( lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]), if (!is.null(RunOptions$IniResLevels)) {
list(CemaNeigeLayers) ); RunOptions$IniStates[1] <- RunOptions$IniResLevels[1] * ParamMod[1] ### production store level (mm)
names(OutputsModel) <- c(FortranOutputsMod[IndOutputsMod],NameCemaNeigeLayers); } RunOptions$IniStates[2] <- RunOptions$IniResLevels[2] * ParamMod[3] ### routing store level (mm)
##DatesR_and_OutputsModel_only RunOptions$IniStates[3] <- RunOptions$IniResLevels[3] ### exponential store level (mm)
if(ExportDatesR==TRUE & ExportStateEnd==FALSE){ }
OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]), ## Call GR model Fortan
list(CemaNeigeLayers) ); RESULTS <- .Fortran("frun_gr6j", PACKAGE = "airGR",
names(OutputsModel) <- c("DatesR",FortranOutputsMod[IndOutputsMod],NameCemaNeigeLayers); } ## inputs
##OutputsModel_and_SateEnd_only LInputs = LInputSeries, ### length of input and output series
if(ExportDatesR==FALSE & ExportStateEnd==TRUE){ InputsPrecip = CatchMeltAndPliq, ### input series of total precipitation [mm/d]
OutputsModel <- c( lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]), InputsPE = InputsModel$PotEvap[IndPeriod1], ### input series potential evapotranspiration [mm/d]
list(CemaNeigeLayers), NParam = NParamMod, ### number of model parameter
list(c(RESULTS$StateEnd,CemaNeigeStateEnd)) ); Param = ParamMod, ### parameter set
names(OutputsModel) <- c(FortranOutputsMod[IndOutputsMod],NameCemaNeigeLayers,"StateEnd"); } NStates = NStatesMod, ### number of state variables used for model initialising
##DatesR_and_OutputsModel_and_SateEnd StateStart = RunOptions$IniStates[1:NStatesMod], ### state variables used when the model run starts
if(ExportDatesR==TRUE & ExportStateEnd==TRUE){ NOutputs = as.integer(length(IndOutputsMod)), ### number of output series
OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]), IndOutputs = IndOutputsMod, ### indices of output series
lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]), ## outputs
list(CemaNeigeLayers), Outputs = matrix(as.double(-99e9), nrow = LInputSeries, ncol = length(IndOutputsMod)), ### output series [mm or mm/d]
list(c(RESULTS$StateEnd,CemaNeigeStateEnd)) ); StateEnd = rep(as.double(-99e9), NStatesMod) ### state variables at the end of the model run
names(OutputsModel) <- c("DatesR",FortranOutputsMod[IndOutputsMod],NameCemaNeigeLayers,"StateEnd"); } )
RESULTS$Outputs[RESULTS$Outputs <= -99e8] <- NA
##End RESULTS$StateEnd[RESULTS$StateEnd <= -99e8] <- NA
rm(RESULTS); if (ExportStateEnd) {
class(OutputsModel) <- c("OutputsModel","daily","GR","CemaNeige"); RESULTS$StateEnd[-3L] <- ifelse(RESULTS$StateEnd[-3L] < 0, 0, RESULTS$StateEnd[-3L]) ### remove negative values except for the ExpStore location
return(OutputsModel); idNStates <- seq_len(NStates*NLayers) %% NStates
RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_CemaNeigeGR6J, InputsModel = InputsModel, IsHyst = IsHyst,
ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = RESULTS$StateEnd[3L],
UH1 = RESULTS$StateEnd[(1:20) + 7],
UH2 = RESULTS$StateEnd[(1:40) + (7+20)],
GCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 1]],
eTGCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 2]],
GthrCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 3]],
GlocmaxCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 0]],
verbose = FALSE)
}
if (inherits(RunOptions, "CemaNeige") & "Precip" %in% RunOptions$Outputs_Sim) {
RESULTS$Outputs[, which(RunOptions$FortranOutputs$GR[IndOutputsMod] == "Precip")] <-
InputsModel$Precip[IndPeriod1]
}
## OutputsModel generation
.GetOutputsModelGR(InputsModel,
RunOptions,
RESULTS,
LInputSeries,
Param,
CemaNeigeLayers)
} }
R/RunModel_GR1A.R
View file @ 4bccfee2
RunModel_GR1A <- function(InputsModel,RunOptions,Param){ RunModel_GR1A <- function(InputsModel, RunOptions, Param) {
NParam <- 1; .ArgumentsCheckGR(InputsModel, RunOptions, Param)
FortranOutputs <- c("PotEvap","Precip","Qsim");
Param <- as.double(Param)
##Arguments_check
if(inherits(InputsModel,"InputsModel")==FALSE){ stop("InputsModel must be of class 'InputsModel' \n"); return(NULL); }
if(inherits(InputsModel,"yearly" )==FALSE){ stop("InputsModel must be of class 'yearly' \n"); return(NULL); } ## Input data preparation
if(inherits(InputsModel,"GR" )==FALSE){ stop("InputsModel must be of class 'GR' \n"); return(NULL); } if (identical(RunOptions$IndPeriod_WarmUp, 0L)) {
if(inherits(RunOptions,"RunOptions" )==FALSE){ stop("RunOptions must be of class 'RunOptions' \n"); return(NULL); } RunOptions$IndPeriod_WarmUp <- NULL
if(inherits(RunOptions,"GR" )==FALSE){ stop("RunOptions must be of class 'GR' \n"); return(NULL); } }
if(!is.vector(Param)){ stop("Param must be a vector \n"); return(NULL); } IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp, RunOptions$IndPeriod_Run)
if(sum(!is.na(Param))!=NParam){ stop(paste("Param must be a vector of length ",NParam," and contain no NA \n",sep="")); return(NULL); } LInputSeries <- as.integer(length(IndPeriod1))
Param <- as.double(Param); if ("all" %in% RunOptions$Outputs_Sim) {
IndOutputs <- as.integer(1:length(RunOptions$FortranOutputs$GR))
##Input_data_preparation } else {
if(identical(RunOptions$IndPeriod_WarmUp,as.integer(0))){ RunOptions$IndPeriod_WarmUp <- NULL; } IndOutputs <- which(RunOptions$FortranOutputs$GR %in% RunOptions$Outputs_Sim)
IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp,RunOptions$IndPeriod_Run); }
LInputSeries <- as.integer(length(IndPeriod1))
if("all" %in% RunOptions$Outputs_Sim){ IndOutputs <- as.integer(1:length(FortranOutputs));
} else { IndOutputs <- which(FortranOutputs %in% RunOptions$Outputs_Sim); } ## Output data preparation
IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries
ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim
BOOL_Fortran <- FALSE; if(BOOL_Fortran){ ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim
##Call_fortan
RESULTS <- .Fortran("frun_gr1a",PACKAGE="airGR",
##inputs ## Call GR model Fortan
LInputs=LInputSeries, ### length of input and output series RESULTS <- .Fortran("frun_gr1a", PACKAGE = "airGR",
InputsPrecip=InputsModel$Precip[IndPeriod1], ### input series of total precipitation [mm/y] ## inputs
InputsPE=InputsModel$PotEvap[IndPeriod1], ### input series potential evapotranspiration [mm/y] LInputs = LInputSeries, ### length of input and output series
NParam=as.integer(length(Param)), ### number of model parameter InputsPrecip = InputsModel$Precip[IndPeriod1], ### input series of total precipitation [mm/y]
Param=Param, ### parameter set InputsPE = InputsModel$PotEvap[IndPeriod1], ### input series potential evapotranspiration [mm/y]
NStates=as.integer(length(RunOptions$IniStates)), ### number of state variables used for model initialising NParam = as.integer(length(Param)), ### number of model parameter
StateStart=RunOptions$IniStates, ### state variables used when the model run starts Param = Param, ### parameter set
NOutputs=as.integer(length(IndOutputs)), ### number of output series NStates = as.integer(length(RunOptions$IniStates)), ### number of state variables used for model initialising
IndOutputs=IndOutputs, ### indices of output series StateStart = RunOptions$IniStates, ### state variables used when the model run starts
##outputs NOutputs = as.integer(length(IndOutputs)), ### number of output series
Outputs=matrix(as.double(-999.999),nrow=LInputSeries,ncol=length(IndOutputs)), ### output series [mm] IndOutputs = IndOutputs, ### indices of output series
StateEnd=rep(as.double(-999.999),length(RunOptions$IniStates)) ### state variables at the end of the model run ## outputs
) Outputs = matrix(as.double(-99e9), nrow = LInputSeries, ncol = length(IndOutputs)), ### output series [mm/y]
RESULTS$Outputs[ round(RESULTS$Outputs ,3)==(-999.999)] <- NA; StateEnd = rep(as.double(-99e9), length(RunOptions$IniStates)) ### state variables at the end of the model run
RESULTS$StateEnd[round(RESULTS$StateEnd,3)==(-999.999)] <- NA; )
RESULTS$Outputs[RESULTS$Outputs <= -99e8] <- NA
} else { RESULTS$StateEnd[RESULTS$StateEnd <= -99e8] <- NA
##R_version
L <- length(IndPeriod1) ## OutputsModel generation
P0 <- InputsModel$Precip[ IndPeriod1][1:(L-1)] .GetOutputsModelGR(InputsModel,
P1 <- InputsModel$Precip[ IndPeriod1][2: L ] RunOptions,
E1 <- InputsModel$PotEvap[IndPeriod1][2: L ] RESULTS,
Q1 <- P1*(1.-1./(1.+((0.7*P1+0.3*P0)/Param[1]/E1)^2.0)^0.5) LInputSeries,
PEQ <- rbind(c(NA,NA,NA),cbind(P1,E1,Q1)) Param)
Outputs <- PEQ[,IndOutputs]
if(is.vector(Outputs)){ Outputs <- cbind(Outputs); }
RESULTS <- list(NOutputs=length(IndOutputs),IndOutputs=IndOutputs,Outputs=Outputs,StatesEnd=NA)
}
##Output_data_preparation
IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries;
ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim;
ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim;
##OutputsModel_only
if(ExportDatesR==FALSE & ExportStateEnd==FALSE){
OutputsModel <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]);
names(OutputsModel) <- FortranOutputs[IndOutputs]; }
##DatesR_and_OutputsModel_only
if(ExportDatesR==TRUE & ExportStateEnd==FALSE){
OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]) );
names(OutputsModel) <- c("DatesR",FortranOutputs[IndOutputs]); }
##OutputsModel_and_SateEnd_only
if(ExportDatesR==FALSE & ExportStateEnd==TRUE){
OutputsModel <- c( lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
list(RESULTS$StateEnd) );
names(OutputsModel) <- c(FortranOutputs[IndOutputs],"StateEnd"); }
##DatesR_and_OutputsModel_and_SateEnd
if((ExportDatesR==TRUE & ExportStateEnd==TRUE) | "all" %in% RunOptions$Outputs_Sim){
OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
list(RESULTS$StateEnd) );
names(OutputsModel) <- c("DatesR",FortranOutputs[IndOutputs],"StateEnd"); }
##End
rm(RESULTS);
class(OutputsModel) <- c("OutputsModel","yearly","GR");
return(OutputsModel);
} }
R/RunModel_GR2M.R
View file @ 4bccfee2
RunModel_GR2M <- function(InputsModel,RunOptions,Param){ RunModel_GR2M <- function(InputsModel, RunOptions, Param) {
NParam <- 2; .ArgumentsCheckGR(InputsModel, RunOptions, Param)
FortranOutputs <- c("PotEvap","Precip","AE","Perc","P3","Exch","Prod","Rout","Qsim");
##Arguments_check Param <- as.double(Param)
if(inherits(InputsModel,"InputsModel")==FALSE){ stop("InputsModel must be of class 'InputsModel' \n"); return(NULL); }
if(inherits(InputsModel,"monthly" )==FALSE){ stop("InputsModel must be of class 'monthly' \n"); return(NULL); }
if(inherits(InputsModel,"GR" )==FALSE){ stop("InputsModel must be of class 'GR' \n"); return(NULL); }
if(inherits(RunOptions,"RunOptions" )==FALSE){ stop("RunOptions must be of class 'RunOptions' \n"); return(NULL); }
if(inherits(RunOptions,"GR" )==FALSE){ stop("RunOptions must be of class 'GR' \n"); return(NULL); }
if(!is.vector(Param)){ stop("Param must be a vector \n"); return(NULL); }
if(sum(!is.na(Param))!=NParam){ stop(paste("Param must be a vector of length ",NParam," and contain no NA \n",sep="")); return(NULL); }
Param <- as.double(Param);
Param_X1_threshold <- 1e-2
if (Param[1L] < Param_X1_threshold) {
warning(sprintf("Param[1] (X1: production store capacity [mm]) < %.2f\n X1 set to %.2f", Param_X1_threshold, Param_X1_threshold))
Param[1L] <- Param_X1_threshold
}
##Input_data_preparation Param_X1X2_threshold <- 1e-2
if(identical(RunOptions$IndPeriod_WarmUp,as.integer(0))){ RunOptions$IndPeriod_WarmUp <- NULL; } if (Param[1L] < Param_X1X2_threshold) {
IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp,RunOptions$IndPeriod_Run); warning(sprintf("Param[1] (X1: production store capacity [mm]) < %.2f\n X1 set to %.2f", Param_X1X2_threshold, Param_X1X2_threshold))
LInputSeries <- as.integer(length(IndPeriod1)) Param[1L] <- Param_X1X2_threshold
if("all" %in% RunOptions$Outputs_Sim){ IndOutputs <- as.integer(1:length(FortranOutputs)); }
} else { IndOutputs <- which(FortranOutputs %in% RunOptions$Outputs_Sim); } if (Param[2L] < Param_X1X2_threshold) {
warning(sprintf("Param[2] (X2: routing store capacity [mm]) < %.2f\n X2 set to %.2f", Param_X1X2_threshold, Param_X1X2_threshold))
Param[2L] <- Param_X1X2_threshold
}
##Use_of_IniResLevels ## Input data preparation
if("IniResLevels" %in% names(RunOptions)){ if (identical(RunOptions$IndPeriod_WarmUp, 0L)) {
RunOptions$IniStates[1] <- RunOptions$IniResLevels[1]*Param[1]; ### production store level (mm) RunOptions$IndPeriod_WarmUp <- NULL
} }
IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp, RunOptions$IndPeriod_Run)
LInputSeries <- as.integer(length(IndPeriod1))
if ("all" %in% RunOptions$Outputs_Sim) {
IndOutputs <- as.integer(1:length(RunOptions$FortranOutputs$GR))
} else {
IndOutputs <- which(RunOptions$FortranOutputs$GR %in% RunOptions$Outputs_Sim)
}
##Call_fortan ## Output data preparation
RESULTS <- .Fortran("frun_GR2M",PACKAGE="airGR", IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries
##inputs ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim
LInputs=LInputSeries, ### length of input and output series ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim
InputsPrecip=InputsModel$Precip[IndPeriod1], ### input series of total precipitation [mm/month]
InputsPE=InputsModel$PotEvap[IndPeriod1], ### input series potential evapotranspiration [mm/month]
NParam=as.integer(length(Param)), ### number of model parameter
Param=Param, ### parameter set
NStates=as.integer(length(RunOptions$IniStates)), ### number of state variables used for model initialising
StateStart=RunOptions$IniStates, ### state variables used when the model run starts
NOutputs=as.integer(length(IndOutputs)), ### number of output series
IndOutputs=IndOutputs, ### indices of output series
##outputs
Outputs=matrix(as.double(-999.999),nrow=LInputSeries,ncol=length(IndOutputs)), ### output series [mm]
StateEnd=rep(as.double(-999.999),length(RunOptions$IniStates)) ### state variables at the end of the model run
)
RESULTS$Outputs [round(RESULTS$Outputs ,3)==(-999.999)] <- NA;
RESULTS$StateEnd[round(RESULTS$StateEnd,3)==(-999.999)] <- NA;
##Output_data_preparation ## Use of IniResLevels
IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries; if (!is.null(RunOptions$IniResLevels)) {
ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim; RunOptions$IniStates[1] <- RunOptions$IniResLevels[1] * Param[1] ### production store level (mm)
ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim; RunOptions$IniStates[2] <- RunOptions$IniResLevels[2] * Param[2] ### routing store level (mm)
##OutputsModel_only }
if(ExportDatesR==FALSE & ExportStateEnd==FALSE){
OutputsModel <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]);
names(OutputsModel) <- FortranOutputs[IndOutputs]; }
##DatesR_and_OutputsModel_only
if(ExportDatesR==TRUE & ExportStateEnd==FALSE){
OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]) );
names(OutputsModel) <- c("DatesR",FortranOutputs[IndOutputs]); }
##OutputsModel_and_SateEnd_only
if(ExportDatesR==FALSE & ExportStateEnd==TRUE){
OutputsModel <- c( lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
list(RESULTS$StateEnd) );
names(OutputsModel) <- c(FortranOutputs[IndOutputs],"StateEnd"); }
##DatesR_and_OutputsModel_and_SateEnd
if((ExportDatesR==TRUE & ExportStateEnd==TRUE) | "all" %in% RunOptions$Outputs_Sim){
OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
list(RESULTS$StateEnd) );
names(OutputsModel) <- c("DatesR",FortranOutputs[IndOutputs],"StateEnd"); }
##End ## Call GR model Fortan
rm(RESULTS); RESULTS <- .Fortran("frun_gr2m", PACKAGE = "airGR",
class(OutputsModel) <- c("OutputsModel","monthly","GR"); ## inputs
return(OutputsModel); LInputs = LInputSeries, ### length of input and output series
InputsPrecip = InputsModel$Precip[IndPeriod1], ### input series of total precipitation [mm/month]
InputsPE = InputsModel$PotEvap[IndPeriod1], ### input series potential evapotranspiration [mm/month]
NParam = as.integer(length(Param)), ### number of model parameter
Param = Param, ### parameter set
NStates = as.integer(length(RunOptions$IniStates)), ### number of state variables used for model initialising
StateStart = RunOptions$IniStates, ### state variables used when the model run starts
NOutputs = as.integer(length(IndOutputs)), ### number of output series
IndOutputs = IndOutputs, ### indices of output series
## outputs
Outputs = matrix(as.double(-99e9), nrow = LInputSeries, ncol = length(IndOutputs)), ### output series [mm or mm/month]
StateEnd = rep(as.double(-99e9), length(RunOptions$IniStates)) ### state variables at the end of the model run
)
RESULTS$Outputs[RESULTS$Outputs <= -99e8] <- NA
RESULTS$StateEnd[RESULTS$StateEnd <= -99e8] <- NA
if (ExportStateEnd) {
RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_GR2M, InputsModel = InputsModel,
ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL,
UH1 = NULL,
UH2 = NULL,
GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL,
verbose = FALSE)
}
## Output data preparation
## OutputsModel generation
.GetOutputsModelGR(InputsModel,
RunOptions,
RESULTS,
LInputSeries,
Param)
} }
R/RunModel_GR4H.R
View file @ 4bccfee2
RunModel_GR4H <- function(InputsModel,RunOptions,Param){ RunModel_GR4H <- function(InputsModel, RunOptions, Param) {
NParam <- 4; .ArgumentsCheckGR(InputsModel, RunOptions, Param)
FortranOutputs <- c("PotEvap","Precip","Prod","AE","Perc","PR","Q9","Q1","Rout","Exch","AExch","QR","QD","Qsim");
##Arguments_check Param <- as.double(Param)
if(inherits(InputsModel,"InputsModel")==FALSE){ stop("InputsModel must be of class 'InputsModel' \n"); return(NULL); }
if(inherits(InputsModel,"hourly" )==FALSE){ stop("InputsModel must be of class 'hourly' \n"); return(NULL); }
if(inherits(InputsModel,"GR" )==FALSE){ stop("InputsModel must be of class 'GR' \n"); return(NULL); }
if(inherits(RunOptions,"RunOptions" )==FALSE){ stop("RunOptions must be of class 'RunOptions' \n"); return(NULL); }
if(inherits(RunOptions,"GR" )==FALSE){ stop("RunOptions must be of class 'GR' \n"); return(NULL); }
if(!is.vector(Param)){ stop("Param must be a vector \n"); return(NULL); }
if(sum(!is.na(Param))!=NParam){ stop(paste("Param must be a vector of length ",NParam," and contain no NA \n",sep="")); return(NULL); }
Param <- as.double(Param);
Param_X1X3_threshold <- 1e-2
if (Param[1L] < Param_X1X3_threshold) {
warning(sprintf("Param[1] (X1: production store capacity [mm]) < %.2f\n X1 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
Param[1L] <- Param_X1X3_threshold
}
if (Param[3L] < Param_X1X3_threshold) {
warning(sprintf("Param[3] (X3: routing store capacity [mm]) < %.2f\n X3 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
Param[3L] <- Param_X1X3_threshold
}
##Input_data_preparation Param_X1X3_threshold <- 1e-2
if(identical(RunOptions$IndPeriod_WarmUp,as.integer(0))){ RunOptions$IndPeriod_WarmUp <- NULL; } Param_X4_threshold <- 0.5
IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp,RunOptions$IndPeriod_Run); if (Param[1L] < Param_X1X3_threshold) {
LInputSeries <- as.integer(length(IndPeriod1)) warning(sprintf("Param[1] (X1: production store capacity [mm]) < %.2f\n X1 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
if("all" %in% RunOptions$Outputs_Sim){ IndOutputs <- as.integer(1:length(FortranOutputs)); Param[1L] <- Param_X1X3_threshold
} else { IndOutputs <- which(FortranOutputs %in% RunOptions$Outputs_Sim); } }
if (Param[3L] < Param_X1X3_threshold) {
warning(sprintf("Param[3] (X3: routing store capacity [mm]) < %.2f\n X3 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
Param[3L] <- Param_X1X3_threshold
}
if (Param[4L] < Param_X4_threshold) {
warning(sprintf("Param[4] (X4: unit hydrograph time constant [h]) < %.2f\n X4 set to %.2f", Param_X4_threshold, Param_X4_threshold))
Param[4L] <- Param_X4_threshold
}
##Use_of_IniResLevels ## Input data preparation
if("IniResLevels" %in% names(RunOptions)){ if (identical(RunOptions$IndPeriod_WarmUp, 0L)) {
RunOptions$IniStates[1] <- RunOptions$IniResLevels[1]*Param[1]; ### production store level (mm) RunOptions$IndPeriod_WarmUp <- NULL
RunOptions$IniStates[2] <- RunOptions$IniResLevels[2]*Param[3]; ### routing store level (mm) }
} IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp, RunOptions$IndPeriod_Run)
LInputSeries <- as.integer(length(IndPeriod1))
if ("all" %in% RunOptions$Outputs_Sim) {
IndOutputs <- as.integer(1:length(RunOptions$FortranOutputs$GR))
} else {
IndOutputs <- which(RunOptions$FortranOutputs$GR %in% RunOptions$Outputs_Sim)
}
##Call_fortan ## Output data preparation
RESULTS <- .Fortran("frun_GR4H",PACKAGE="airGR", IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries
##inputs ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim
LInputs=LInputSeries, ### length of input and output series ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim
InputsPrecip=InputsModel$Precip[IndPeriod1], ### input series of total precipitation [mm/h]
InputsPE=InputsModel$PotEvap[IndPeriod1], ### input series potential evapotranspiration [mm/h]
NParam=as.integer(length(Param)), ### number of model parameter
Param=Param, ### parameter set
NStates=as.integer(length(RunOptions$IniStates)), ### number of state variables used for model initialising
StateStart=RunOptions$IniStates, ### state variables used when the model run starts
NOutputs=as.integer(length(IndOutputs)), ### number of output series
IndOutputs=IndOutputs, ### indices of output series
##outputs
Outputs=matrix(as.double(-999.999),nrow=LInputSeries,ncol=length(IndOutputs)), ### output series [mm]
StateEnd=rep(as.double(-999.999),length(RunOptions$IniStates)) ### state variables at the end of the model run
)
RESULTS$Outputs[ round(RESULTS$Outputs ,3)==(-999.999)] <- NA;
RESULTS$StateEnd[round(RESULTS$StateEnd,3)==(-999.999)] <- NA;
##Output_data_preparation ## Use of IniResLevels
IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries; if (!is.null(RunOptions$IniResLevels)) {
ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim; RunOptions$IniStates[1] <- RunOptions$IniResLevels[1] * Param[1] ### production store level (mm)
ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim; RunOptions$IniStates[2] <- RunOptions$IniResLevels[2] * Param[3] ### routing store level (mm)
##OutputsModel_only }
if(ExportDatesR==FALSE & ExportStateEnd==FALSE){
OutputsModel <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]);
names(OutputsModel) <- FortranOutputs[IndOutputs]; }
##DatesR_and_OutputsModel_only
if(ExportDatesR==TRUE & ExportStateEnd==FALSE){
OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]) );
names(OutputsModel) <- c("DatesR",FortranOutputs[IndOutputs]); }
##OutputsModel_and_SateEnd_only
if(ExportDatesR==FALSE & ExportStateEnd==TRUE){
OutputsModel <- c( lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
list(RESULTS$StateEnd) );
names(OutputsModel) <- c(FortranOutputs[IndOutputs],"StateEnd"); }
##DatesR_and_OutputsModel_and_SateEnd
if((ExportDatesR==TRUE & ExportStateEnd==TRUE) | "all" %in% RunOptions$Outputs_Sim){
OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
list(RESULTS$StateEnd) );
names(OutputsModel) <- c("DatesR",FortranOutputs[IndOutputs],"StateEnd"); }
##End ## Call GR model Fortan
rm(RESULTS); RESULTS <- .Fortran("frun_gr4h", PACKAGE = "airGR",
class(OutputsModel) <- c("OutputsModel","hourly","GR"); ## inputs
return(OutputsModel); LInputs = LInputSeries, ### length of input and output series
InputsPrecip = InputsModel$Precip[IndPeriod1], ### input series of total precipitation [mm/h]
InputsPE = InputsModel$PotEvap[IndPeriod1], ### input series potential evapotranspiration [mm/h]
NParam = as.integer(length(Param)), ### number of model parameter
Param = Param, ### parameter set
NStates = as.integer(length(RunOptions$IniStates)), ### number of state variables used for model initialising
StateStart = RunOptions$IniStates, ### state variables used when the model run starts
NOutputs = as.integer(length(IndOutputs)), ### number of output series
IndOutputs = IndOutputs, ### indices of output series
## outputs
Outputs = matrix(as.double(-99e9), nrow = LInputSeries, ncol = length(IndOutputs)), ### output series [mm or mm/h]
StateEnd = rep(as.double(-99e9), length(RunOptions$IniStates)) ### state variables at the end of the model run
)
RESULTS$Outputs[RESULTS$Outputs <= -99e8] <- NA
RESULTS$StateEnd[RESULTS$StateEnd <= -99e8] <- NA
if (ExportStateEnd) {
RESULTS$StateEnd[-3L] <- ifelse(RESULTS$StateEnd[-3L] < 0, 0, RESULTS$StateEnd[-3L]) ### remove negative values except for the ExpStore location
RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_GR4H, InputsModel = InputsModel,
ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL,
UH1 = RESULTS$StateEnd[(1:(20*24)) + 7],
UH2 = RESULTS$StateEnd[(1:(40*24)) + (7+20*24)],
GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL,
verbose = FALSE)
}
## OutputsModel generation
.GetOutputsModelGR(InputsModel,
RunOptions,
RESULTS,
LInputSeries,
Param)
} }
R/RunModel_GR4J.R
View file @ 4bccfee2
RunModel_GR4J <- function(InputsModel,RunOptions,Param){ RunModel_GR4J <- function(InputsModel, RunOptions, Param) {
NParam <- 4; .ArgumentsCheckGR(InputsModel, RunOptions, Param)
FortranOutputs <- c("PotEvap","Precip","Prod","AE","Perc","PR","Q9","Q1","Rout","Exch","AExch","QR","QD","Qsim");
##Arguments_check Param <- as.double(Param)
if(inherits(InputsModel,"InputsModel")==FALSE){ stop("InputsModel must be of class 'InputsModel' \n"); return(NULL); }
if(inherits(InputsModel,"daily" )==FALSE){ stop("InputsModel must be of class 'daily' \n"); return(NULL); }
if(inherits(InputsModel,"GR" )==FALSE){ stop("InputsModel must be of class 'GR' \n"); return(NULL); }
if(inherits(RunOptions,"RunOptions" )==FALSE){ stop("RunOptions must be of class 'RunOptions' \n"); return(NULL); }
if(inherits(RunOptions,"GR" )==FALSE){ stop("RunOptions must be of class 'GR' \n"); return(NULL); }
if(!is.vector(Param)){ stop("Param must be a vector \n"); return(NULL); }
if(sum(!is.na(Param))!=NParam){ stop(paste("Param must be a vector of length ",NParam," and contain no NA \n",sep="")); return(NULL); }
Param <- as.double(Param);
Param_X1X3_threshold <- 1e-2
if (Param[1L] < Param_X1X3_threshold) {
warning(sprintf("Param[1] (X1: production store capacity [mm]) < %.2f\n X1 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
Param[1L] <- Param_X1X3_threshold
}
if (Param[3L] < Param_X1X3_threshold) {
warning(sprintf("Param[3] (X3: routing store capacity [mm]) < %.2f\n X3 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
Param[3L] <- Param_X1X3_threshold
}
##Input_data_preparation Param_X1X3_threshold <- 1e-2
if(identical(RunOptions$IndPeriod_WarmUp,as.integer(0))){ RunOptions$IndPeriod_WarmUp <- NULL; } Param_X4_threshold <- 0.5
IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp,RunOptions$IndPeriod_Run); if (Param[1L] < Param_X1X3_threshold) {
LInputSeries <- as.integer(length(IndPeriod1)) warning(sprintf("Param[1] (X1: production store capacity [mm]) < %.2f\n X1 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
if("all" %in% RunOptions$Outputs_Sim){ IndOutputs <- as.integer(1:length(FortranOutputs)); Param[1L] <- Param_X1X3_threshold
} else { IndOutputs <- which(FortranOutputs %in% RunOptions$Outputs_Sim); } }
if (Param[3L] < Param_X1X3_threshold) {
warning(sprintf("Param[3] (X3: routing store capacity [mm]) < %.2f\n X3 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
Param[3L] <- Param_X1X3_threshold
}
if (Param[4L] < Param_X4_threshold) {
warning(sprintf("Param[4] (X4: unit hydrograph time constant [d]) < %.2f\n X4 set to %.2f", Param_X4_threshold, Param_X4_threshold))
Param[4L] <- Param_X4_threshold
}
##Use_of_IniResLevels ## Input data preparation
if("IniResLevels" %in% names(RunOptions)){ if (identical(RunOptions$IndPeriod_WarmUp, 0L)) {
RunOptions$IniStates[1] <- RunOptions$IniResLevels[1]*Param[1]; ### production store level (mm) RunOptions$IndPeriod_WarmUp <- NULL
RunOptions$IniStates[2] <- RunOptions$IniResLevels[2]*Param[3]; ### routing store level (mm) }
} IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp, RunOptions$IndPeriod_Run)
LInputSeries <- as.integer(length(IndPeriod1))
if ("all" %in% RunOptions$Outputs_Sim) {
IndOutputs <- as.integer(1:length(RunOptions$FortranOutputs$GR))
} else {
IndOutputs <- which(RunOptions$FortranOutputs$GR %in% RunOptions$Outputs_Sim)
}
##Call_fortan ## Use of IniResLevels
RESULTS <- .Fortran("frun_GR4J",PACKAGE="airGR", if (!is.null(RunOptions$IniResLevels)) {
##inputs RunOptions$IniStates[1] <- RunOptions$IniResLevels[1] * Param[1] ### production store level (mm)
LInputs=LInputSeries, ### length of input and output series RunOptions$IniStates[2] <- RunOptions$IniResLevels[2] * Param[3] ### routing store level (mm)
InputsPrecip=InputsModel$Precip[IndPeriod1], ### input series of total precipitation [mm/d] }
InputsPE=InputsModel$PotEvap[IndPeriod1], ### input series potential evapotranspiration [mm/d]
NParam=as.integer(length(Param)), ### number of model parameter
Param=Param, ### parameter set
NStates=as.integer(length(RunOptions$IniStates)), ### number of state variables used for model initialising
StateStart=RunOptions$IniStates, ### state variables used when the model run starts
NOutputs=as.integer(length(IndOutputs)), ### number of output series
IndOutputs=IndOutputs, ### indices of output series
##outputs
Outputs=matrix(as.double(-999.999),nrow=LInputSeries,ncol=length(IndOutputs)), ### output series [mm]
StateEnd=rep(as.double(-999.999),length(RunOptions$IniStates)) ### state variables at the end of the model run
)
RESULTS$Outputs[ round(RESULTS$Outputs ,3)==(-999.999)] <- NA;
RESULTS$StateEnd[round(RESULTS$StateEnd,3)==(-999.999)] <- NA;
##Output_data_preparation ## Call GR model Fortan
IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries; RESULTS <- .Fortran("frun_gr4j", PACKAGE = "airGR",
ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim; ## inputs
ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim; LInputs = LInputSeries, ### length of input and output series
##OutputsModel_only InputsPrecip = InputsModel$Precip[IndPeriod1], ### input series of total precipitation [mm/d]
if(ExportDatesR==FALSE & ExportStateEnd==FALSE){ InputsPE = InputsModel$PotEvap[IndPeriod1], ### input series potential evapotranspiration [mm/d]
OutputsModel <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]); NParam = as.integer(length(Param)), ### number of model parameter
names(OutputsModel) <- FortranOutputs[IndOutputs]; } Param = Param, ### parameter set
##DatesR_and_OutputsModel_only NStates = as.integer(length(RunOptions$IniStates)), ### number of state variables used for model initialising
if(ExportDatesR==TRUE & ExportStateEnd==FALSE){ StateStart = RunOptions$IniStates, ### state variables used when the model run starts
OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]), NOutputs = as.integer(length(IndOutputs)), ### number of output series
lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]) ); IndOutputs = IndOutputs, ### indices of output series
names(OutputsModel) <- c("DatesR",FortranOutputs[IndOutputs]); } ## outputs
##OutputsModel_and_SateEnd_only Outputs = matrix(as.double(-99e9), nrow = LInputSeries, ncol = length(IndOutputs)), ### output series [mm or mm/d]
if(ExportDatesR==FALSE & ExportStateEnd==TRUE){ StateEnd = rep(as.double(-99e9), length(RunOptions$IniStates)) ### state variables at the end of the model run
OutputsModel <- c( lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]), )
list(RESULTS$StateEnd) ); RESULTS$Outputs[RESULTS$Outputs <= -99e8] <- NA
names(OutputsModel) <- c(FortranOutputs[IndOutputs],"StateEnd"); } RESULTS$StateEnd[RESULTS$StateEnd <= -99e8] <- NA
##DatesR_and_OutputsModel_and_SateEnd if ("StateEnd" %in% RunOptions$Outputs_Sim) {
if((ExportDatesR==TRUE & ExportStateEnd==TRUE) | "all" %in% RunOptions$Outputs_Sim){ RESULTS$StateEnd[-3L] <- ifelse(RESULTS$StateEnd[-3L] < 0, 0, RESULTS$StateEnd[-3L]) ### remove negative values except for the ExpStore location
OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]), RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_GR4J, InputsModel = InputsModel,
lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]), ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL,
list(RESULTS$StateEnd) ); UH1 = RESULTS$StateEnd[(1:20) + 7],
names(OutputsModel) <- c("DatesR",FortranOutputs[IndOutputs],"StateEnd"); } UH2 = RESULTS$StateEnd[(1:40) + (7+20)],
GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL,
##End verbose = FALSE)
rm(RESULTS); }
class(OutputsModel) <- c("OutputsModel","daily","GR");
return(OutputsModel);
## OutputsModel generation
.GetOutputsModelGR(InputsModel,
RunOptions,
RESULTS,
LInputSeries,
Param)
} }