Vectorization of snow model; snow related SGS fix

kuadrat · Aug 21, 2024 · 5eec56f · 5eec56f
1 parent ad04bf6
commit 5eec56f
Show file tree

Hide file tree

Showing 4 changed files with 29 additions and 22 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: growR
 Type: Package
-Version: 1.3.0.9004
+Version: 1.3.0.9005
 Date: 2024-05-23
 Authors@R: person(
     given = "Kevin",

diff --git a/NEWS.md b/NEWS.md
@@ -13,10 +13,14 @@
 * autocut accepts 'medium' as alias for 'middle' and applies a single cut 
   according to phenology.
 * Performance improvements to ModvegeSite$un() in the order of 25%.
+* Small performance improvement of `WeatherData` initialization through 
+  vectorization of snow model.
 
 ## Fixed
 
 * `ModvegeParameters$set_parameters()` can now handle empty lists.
+* Fix in determination of start of growing season for cases with no snow in 
+  first half-year.
 
 ## Removed
 

diff --git a/R/modvegesite.R b/R/modvegesite.R
@@ -665,11 +665,13 @@ ModvegeSite = R6Class(
       W = self$get_weather()
       # Find the last day of the first half of the year that still had snow 
       # cover. Assume that growth can only start after that day.
+      j_snow = 1
       if (consider_snow) {
         # Go one day beyond the last day with too much snow
-        j_snow = max(which(W[["snow"]][1:180] > critical_snow)) + 1
-      } else {
-        j_snow = 1
+        snow_days = which(W[["snow"]][1:180] > critical_snow)
+        if (length(snow_days) > 0) {
+          j_snow = max(snow_days) + 1
+        }
       }
 
       if (private$SGS_method == "MTD") {

diff --git a/R/weather.R b/R/weather.R
@@ -157,27 +157,28 @@ WeatherData = R6Class(
       liquidP_vec = (1./(1. + exp(-1.5*(Ta_vec - 2.)))) * PP_vec
       solidP_vec  = PP_vec - liquidP_vec
 
-      melt_vec   = numeric(vec_size)
-      freeze_vec = numeric(vec_size)
-      snow_vec   = numeric(vec_size)
-      melt_vec[1]   = 0.
-      freeze_vec[1] = 0.
-      snow_vec[1]   = 0.
+      # Calculate as much as possible vectorially.
+      dT_melt = Ta_vec - T_melt
+      T_above_melting = dT_melt >= 0
+      freeze_vec = pmin(liquidP_vec, -C_freeze*dT_melt)
+      freeze_vec[T_above_melting] = 0
+      melt0 = C_melt * dT_melt
+      dsnow0 = solidP_vec + freeze_vec
 
+      snow_vec = numeric(vec_size)
+      melt_vec = numeric(vec_size)
+      snow_vec[1] = 0.
+      melt_vec[1] = 0.
       for (j in 2:vec_size) {
-        # :TODO: Is time_step needed here?
-        time_step = 1
-        melt_vec[j]   = ifelse(snow_vec[j - 1] > 0. & Ta_vec[j] >= T_melt,
-                               min(snow_vec[j - 1]/time_step, 
-                                   C_melt * (Ta_vec[j] - T_melt)), 
-                               0.)
-        freeze_vec[j] = ifelse(Ta_vec[j] < T_melt, 
-                               min(liquidP_vec[j], 
-                                   C_freeze * (T_melt - Ta_vec[j])), 
-                               0.)
-        dsnow_dt     = solidP_vec[j] + freeze_vec[j] - melt_vec[j]
-        snow_vec[j]   = max(0., snow_vec[j - 1] + dsnow_dt * time_step)
+        if (snow_vec[j-1] > 0. & T_above_melting[j]) {
+          melt_vec[j] = min(snow_vec[j-1], melt0[j])
+        } else {
+          melt_vec[j] = 0
+        }
+        dsnow = dsnow0[j] - melt_vec[j]
+        snow_vec[j] = snow_vec[j-1] + dsnow
       }
+
       self[["liquidP_vec"]] = liquidP_vec
       self[["melt_vec"]] = melt_vec
       self[["snow_vec"]] = snow_vec