DataAltiExtrapolation_Valery.R 5.72 KB
Newer Older
1
2
3
4
5
DataAltiExtrapolation_Valery <- function(DatesR,
                                         Precip,  PrecipScale = TRUE,
                                         TempMean, TempMin = NULL, TempMax = NULL,
                                         ZInputs,  HypsoData, NLayers,
                                         verbose = TRUE) {
6
7
8

  ##Altitudinal_gradient_functions_______________________________________________________________
  ##unique_gradient_for_precipitation
9
  GradP_Valery2010 <- 0.00041  ### value from Valery PhD thesis page 126
10
11
12



13
  ##Format_______________________________________________________________________________________
14
15
16
17
18
  HypsoData <- as.double(HypsoData)
  ZInputs   <- as.double(ZInputs)



19
  ##ElevationLayers_Creation_____________________________________________________________________
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
  ZLayers   <- as.double(rep(NA, NLayers))

  if (!identical(HypsoData, as.double(rep(NA, 101)))) {
    nmoy   <- 100 %/% NLayers
    nreste <- 100 %% NLayers
    ncont  <- 0

    for (iLayer in 1:NLayers) {
      if (nreste > 0) {
        nn <- nmoy + 1
        nreste <- nreste - 1
      } else {
        nn <- nmoy
      }
      if (nn == 1) {
        ZLayers[iLayer] <- HypsoData[ncont + 1]
36
      }
37
38
39
40
41
42
43
      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
44
    }
45
46
47
  }


48
49
  ##Precipitation_extrapolation__________________________________________________________________
  ##Initialisation
50
51
52
53
  if (identical(ZInputs, HypsoData[51]) & NLayers == 1) {
    LayerPrecip <- list(as.double(Precip))
  } else {
    ##Elevation_gradients_for_daily_mean_precipitation
54
    GradP    <- GradP_Valery2010 ### single value
55
56
57
58
59
60
61
62
63
64
65
66
67
68
    TabGradP <- rep(GradP, length(Precip))
    ##Extrapolation
    ##Thresold_of_inputs_median_elevation
    Zthreshold <- 4000
    LayerPrecip_mat <- sapply(1:NLayers, function(iLayer) {
      ##If_layer_elevation_smaller_than_Zthreshold
      if (ZLayers[iLayer] <= Zthreshold) {
        prcp <- as.double(Precip * exp(TabGradP * (ZLayers[iLayer] - ZInputs)))
        ##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
69
        } else {
70
          prcp <- as.double(Precip)
Delaigue Olivier's avatar
Delaigue Olivier committed
71
72
        }
      }
73
74
75
76
77
      return(prcp)
    })
    if (PrecipScale) {
      LayerPrecip_mat <- LayerPrecip_mat / rowMeans(LayerPrecip_mat) * Precip
      LayerPrecip_mat[is.nan(LayerPrecip_mat)] <- 0
78
    }
79
80
81
82
83
    LayerPrecip <- as.list(as.data.frame(LayerPrecip_mat))
  }



84
85
  ##Temperature_extrapolation____________________________________________________________________
  ##Initialisation
86
87
88
89
90
91
92
93
94
95
96
97
98
  LayerTempMean <- list()
  LayerTempMin <- list()
  LayerTempMax <- list()

  if (identical(ZInputs, HypsoData[51]) & NLayers == 1) {
    LayerTempMean[[1]] <- as.double(TempMean)

    if (!is.null(TempMin) & !is.null(TempMax)) {
      LayerTempMin[[1]] <- as.double(TempMin)
      LayerTempMax[[1]] <- as.double(TempMax)
    }
  } else {
    ##Elevation_gradients_for_daily_mean_min_and_max_temperature
99
    GradT <- .GradT_Valery2010
100
101
    iday <- match(format(DatesR, format = "%d%m"),
                  sprintf("%02i%02i", GradT[, "day"], GradT[, "month"]))
102
    TabGradT <- GradT[iday, c("grad_Tmean", "grad_Tmin", "grad_Tmax")]
103
104
105
106
    ##Extrapolation
    ##On_each_elevation_layer...
    for (iLayer in 1:NLayers) {
      LayerTempMean[[iLayer]] <- as.double(TempMean + (ZInputs - ZLayers[iLayer]) * abs(TabGradT[, "grad_Tmean"]) /  100)
107
      if (!is.null(TempMin) & !is.null(TempMax)) {
108
109
        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)
Delaigue Olivier's avatar
Delaigue Olivier committed
110
      }
111
    }
112
113
114
115
  }



116
  ##Solid_Fraction_for_each_elevation_layer______________________________________________________
117
118
119
120
121
122
123
124
125
126
  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"
127
    }
128
  }
129

130
131
132
133
134
135
136
137
138
139
  ##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
140
    }
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
    ##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



163
  ##END__________________________________________________________________________________________
164
165
166
167
168
169
170
171
  return(list(LayerPrecip          = LayerPrecip,
              LayerTempMean        = LayerTempMean,
              LayerTempMin         = LayerTempMin,
              LayerTempMax         = LayerTempMax,
              LayerFracSolidPrecip = LayerFracSolidPrecip,
              ZLayers              = ZLayers))


Delaigue Olivier's avatar
Delaigue Olivier committed
172
173
}