Performance improvements

DDDFunctions/Big2SmallLambda.jl

@@ -12,14 +12,11 @@ function Big2SmallLambda(GshInt,GscInt)
 # using LsqFit
 # using Distributions
 # using Statistics
-
- arg = zeros(Float64,10)
- avglam = zeros(Float64,10) # number of levels  
- arg[1:9] = [0.1:0.1:0.9;]
- arg[10] = 0.99
+
+arg = vcat([.1:.1:.9;],[.99])
 
 g = Gamma(GshInt,GscInt)
-avglam[1:10] = quantile.(g,arg)
+avglam = quantile.(g,arg) # number of levels
 
 trlam = zeros(Float64,10)  # Lowercase lambda for level
 trlam[1] = avglam[1]

DDDFunctions/BogLayerUpdate.jl

@@ -7,20 +7,15 @@
 #     Revised: 16.12.2019
 #--------------------------------------------------------------------------
 
-function BogLayerUpdate(outbog, BogLayers, UHBog, nodaysvector)
-
-    qlayer = zeros(nodaysvector)
-    qlayer .= outbog*UHBog    #finds response  in mm!!!!  for bog, a vector     
-    if(nodaysvector > 1)    
-      BogLayers[1:(nodaysvector-1)] .= BogLayers[2:nodaysvector] .+ qlayer[2:nodaysvector]# shifts the level of the matrix one timestep ahead
-      BogLayers[nodaysvector] = 0.0
-    end 
-    if(nodaysvector == 1)
-	  BogLayers[1:nodaysvector] .= 0.0 
-      #BogLayers[1:nodaysvector] .= qlayer
+function BogLayerUpdate!(outbog, BogLayers, UHBog, nodaysvector)
+    qlayer = outbog * UHBog    #finds response  in mm!!!!  for bog, a vector     
+    for t in 2:nodaysvector
+        BogLayers[t-1] = BogLayers[t] + qlayer[t] # shifts the level of the matrix one timestep ahead
+    end
+    if nodaysvector == 1
+        BogLayers .= 0.
+#        BogLayers .= qlayer
     end
-
-return BogLayers
 end            
 
 

DDDFunctions/DDDAllTerrain22012024.jl

@@ -44,7 +44,6 @@ include("LayerEstimation.jl")
 include("PyrAreas.jl")
 include("GrWPoint.jl")
 include("RiverPoint.jl")
-include("TemperatureVector.jl")
 # EB and Snow Routines
 include("NedbEBGlac_debug04072022.jl")
 include("SnowpackTemp.jl")
@@ -171,8 +170,6 @@ snofritt = zeros(Lty)
 wcs = zeros(Lty)
 GIsoil = zeros(Lty)  # mean glacial melt 
 scaobx = zeros(10)
-hprecip = zeros(10)
-htemp = zeros(10)
 hfelt = zeros(10)
 
 ########################Reading parameters#########################################################
@@ -291,9 +288,7 @@ CFR = 2.5*(Timeresinsec/86400)*0.0833   # Fixed as 1/12 of estimate of CX= 2.5 f
 len = Int(5*(86400/Timeresinsec))       # number of timestepes to spin up the model. Recommended to use timesteps equal to a minimum of 5 days to estimate the snowpack temperature.)
 
 startsim = startsim + len               # taking into account estimating snowpack temperature
-tempstart = zeros(len,10)               # matrix for storing temperatures when starting from states
-tempstart = ptqinn[(startsim-len+1):startsim,16:25]  # assigning temperature values
-STempvec = zeros(len)              # Temperature vector for estimating snowpack temperature
+tempstart = Matrix{Float64}(ptqinn[(startsim-len+1):startsim,16:25])  # matrix for storing temperatures when starting from states
 
 Pa[1:10] = 101.3*((293 .- 0.0065.*hfelt[1:10])/293.0).^5.26     # Air pressure as a function of height Zhu et al. 2012 and Stoy et al, 2018
 MPa = mean(Pa)
@@ -507,19 +502,18 @@ for i in startsim:days
   DN = Dates.dayofyear(dato)           #daynumber        
 
   #Reads Precipitation and temperature for each elevation zone   
-  htemp = ptqinn[i,16:25]
-  hprecip = ptqinn[i,6:15] 
+  htemp = Vector{Float64}(ptqinn[i,16:25])
+  hprecip = Vector{Float64}(ptqinn[i,6:15])
 
   meanprecip = mean(hprecip)
   meantemp =  mean(htemp)
-  tempstart = TempstartUpdate(tempstart,htemp, len) #Updating the tempstart with this timesteps temperature 
+  TempstartUpdate!(tempstart,htemp, len) #Updating the tempstart with this timesteps temperature 
 
   for Lst in 1:Lty     # landscape types, one snow regime for each landscape type P and IP. The other Lty have no snow
 
     for idim in 1:hson # elevation zones  
 
-      #STempvec = TemperatureVector(ptqinn[(i-len+1):i,16:25],idim, len) # i is always greater than len
-      STempvec = TemperatureVector(tempstart, idim, len)
+      @views STempvec = reverse(tempstart[:,idim]) # Temperature vector for estimating snowpack temperatur
       CGLAC = 0.0 # dummy, has no role in this version, energy balance is used
       CX = 0.0    # dummy, has no role in this version, energy balance is used
       TS = 0.0    # dummy, has no role in this version, energy balance is used
@@ -776,14 +770,14 @@ for i in startsim:days
   #Updating the saturation Layers
     for Lst in 1:Lty
       if(Lst==1)
-         LayersP = LayerUpdate(ddist[Lst,1:NoL],outx[Lst], LayersP, layerUH_P, nodaysvector[Lst,1:NoL], NoL)
+         LayerUpdate!(ddist[Lst,1:NoL],outx[Lst], LayersP, layerUH_P, nodaysvector[Lst,1:NoL], NoL)
     end       
       if(Lst==2)
-        LayersIP = LayerUpdate(ddist[Lst,1:NoL],outx[Lst], LayersIP, layerUH_IP, nodaysvector[Lst,1:NoL], NoL)       
+        LayerUpdate!(ddist[Lst,1:NoL],outx[Lst], LayersIP, layerUH_IP, nodaysvector[Lst,1:NoL], NoL)       
       end
     end  
 
-  BogLayers = BogLayerUpdate(outbog, BogLayers, UHbog, antBogsteps)#
+  BogLayerUpdate!(outbog, BogLayers, UHbog, antBogsteps)#
 
   #summing up groundwater states
     for Lst in 1:Lty
@@ -858,20 +852,20 @@ for i in startsim:days
   #Qmm = (QRD*Timeresinsec*1000/totarea)  #QRD in mm/day
   #Routing the contributions in the Lake
   Qm3s = LakeLayers[1] + QRD*UHLake[1]            #Lakewater to be discharged from outlet + this timesteps contribution, i.e. catchment response in m3/s
-  LakeLayers = BogLayerUpdate(QRD, LakeLayers, UHLake, nodaysLake) # Same routine updating as for Bogs
+  BogLayerUpdate!(QRD, LakeLayers, UHLake, nodaysLake) # Same routine updating as for Bogs
   Qmm = (Qm3s*Timeresinsec*1000/totarea)  #GDT_Lake in mm/Timeresinsec
 
   P_Qm3s = P_LakeLayers[1] + QRDP*UHLake[1]            #Lakewater to be discharged from outlet + this timesteps contribution, i.e. catchemnt response 
-  P_LakeLayers = BogLayerUpdate(QRDP, P_LakeLayers, UHLake, nodaysLake) # Same routine updating as for Bogs
+  BogLayerUpdate!(QRDP, P_LakeLayers, UHLake, nodaysLake) # Same routine updating as for Bogs
 
   IP_Qm3s = IP_LakeLayers[1] + QRDIP*UHLake[1]            #Lakewater to be discharged from outlet + this timesteps contribution, i.e. catchemnt response 
-  IP_LakeLayers = BogLayerUpdate(QRDIP, IP_LakeLayers, UHLake, nodaysLake) # Same routine updating as for Bogs
+  BogLayerUpdate!(QRDIP, IP_LakeLayers, UHLake, nodaysLake) # Same routine updating as for Bogs
 
   Bog_Qm3s = Bog_LakeLayers[1] + QRDBog*UHLake[1]            #Lakewater to be discharged from outlet + this timesteps contribution, i.e. catchemnt response 
-  Bog_LakeLayers = BogLayerUpdate(QRDBog, Bog_LakeLayers, UHLake, nodaysLake) # Same routine updating as for Bogs
+  BogLayerUpdate!(QRDBog, Bog_LakeLayers, UHLake, nodaysLake) # Same routine updating as for Bogs
 
   OF_Qm3s = OF_LakeLayers[1] + QRDOF*UHLake[1]            #not a contribution, just sttratifying the runoff 
-  OF_LakeLayers = BogLayerUpdate(QRDOF, OF_LakeLayers, UHLake, nodaysLake) # Same routine updating as for Bogs
+  BogLayerUpdate!(QRDOF, OF_LakeLayers, UHLake, nodaysLake) # Same routine updating as for Bogs
 
    #WBP = (SPinn + RinnP) - (lyrs[1] + sum(QRivxP)*(Timeresinsec*1000/area[1]))    #mm  last term inludes discharge and storage in rivernetwork                                         
    #println("WBP = ", WBP, " Total inn P= ",SPinn)                                             
@@ -909,7 +903,7 @@ for i in startsim:days
 #---------------------------------------------------------------------------------------------------
 
 # Updating the  temperature matrix for estimating snowpack temperature
- tempstart = ptqinn[(i-len+1):i,16:25]  # assigning temperature values, i is always larger than len
+ tempstart = Matrix{Float64}(ptqinn[(i-len+1):i,16:25])  # assigning temperature values, i is always larger than len
 
 # Saving states, SWE, sm, Layers, etc in a JLD2 file
 if(i == (38165 + len) && savestate == 1) #  This number is ONE timestep less than startsim, i.e. the statefile is for the timestep before startsim 

DDDFunctions/DDDEventFloodDesign.jl

@@ -390,10 +390,10 @@ for i in 1:days2 #Length of precip timeseries +5 (recession)
   #Updating the saturation Layers. The "boxes are shifted on step ahead"
     for Lst in 1:Lty
       if(Lst==1)
-         LayersP = LayerUpdate(ddist[Lst,1:NoL],outx[Lst], LayersP, layerUH_P, nodaysvector[Lst,1:NoL], NoL)
+         LayerUpdate!(ddist[Lst,1:NoL],outx[Lst], LayersP, layerUH_P, nodaysvector[Lst,1:NoL], NoL)
     end       
       if(Lst==2)
-        LayersIP = LayerUpdate(ddist[Lst,1:NoL],outx[Lst], LayersIP, layerUH_IP, nodaysvector[Lst,1:NoL], NoL)       
+         LayerUpdate!(ddist[Lst,1:NoL],outx[Lst], LayersIP, layerUH_IP, nodaysvector[Lst,1:NoL], NoL)       
       end
     end  
 

DDDFunctions/DDDEventFloodDesignFixedPDynSM.jl

@@ -454,10 +454,10 @@ for i in 1:days2 #Length of precip timeseries +5 (recession)
   #Updating the saturation Layers. The "boxes are shifted on step ahead"
     for Lst in 1:Lty
       if(Lst==1)
-         LayersP = LayerUpdate(ddist[Lst,1:NoL],outx[Lst], LayersP, layerUH_P, nodaysvector[Lst,1:NoL], NoL)
+         LayerUpdate!(ddist[Lst,1:NoL],outx[Lst], LayersP, layerUH_P, nodaysvector[Lst,1:NoL], NoL)
     end       
       if(Lst==2)
-        LayersIP = LayerUpdate(ddist[Lst,1:NoL],outx[Lst], LayersIP, layerUH_IP, nodaysvector[Lst,1:NoL], NoL)       
+        LayerUpdate!(ddist[Lst,1:NoL],outx[Lst], LayersIP, layerUH_IP, nodaysvector[Lst,1:NoL], NoL)       
       end
     end  
 

DDDFunctions/DDDModulFunc-MAD-LAKE.jl

@@ -495,8 +495,8 @@ FromSnow = SnowGamma(PR[idim],PS[idim],MW[idim],sca[idim],spd[idim],wcd[idim],pr
  GDT_Bog = BogLayers[1] + outbog*UHbog[1]            #Bogwater to be discharged into the river network + this timesteps contribution
 
   #Updating the saturation Layers
-  Layers = LayerUpdate(ddist, outx, Layers, layerUH, nodaysvector, NoL) # mm
-  BogLayers = BogLayerUpdate(outbog, BogLayers, UHbog, antBogsteps)     # mm
+  LayerUpdate!(ddist, outx, Layers, layerUH, nodaysvector, NoL) # mm
+  BogLayerUpdate!(outbog, BogLayers, UHbog, antBogsteps)     # mm
 
   lyrs = sum(Layers)# gives the sum of layers after todays runoff has taken place
   boglyrs = sum(BogLayers)
@@ -524,7 +524,7 @@ FromSnow = SnowGamma(PR[idim],PS[idim],MW[idim],sca[idim],spd[idim],wcd[idim],pr
   qmm_state = (sum(QRivx) - QRD)*(Timeresinsec*1000/totarea)       # This is also a reservoir [mm], water from todays event is stored for future runoff in the RN, relevant for WB.  
 
   GDT_Lake = LakeLayers[1] + QRD*UHLake[1]            #Lakewater to be discharged from outlet + this timesteps contribution, i.e. catchemnt response 
-  LakeLayers = BogLayerUpdate(QRD, LakeLayers, UHLake, nodaysLake)#
+  BogLayerUpdate!(QRD, LakeLayers, UHLake, nodaysLake)#
   Qm3s = GDT_Lake # runoff in m3/s
   Qmm = (GDT_Lake*Timeresinsec*1000/totarea)  #GDT_Lake in mm/Timeresinsec
 

DDDFunctions/DDDUrbanFunc10012024_InfWB.jl

@@ -574,14 +574,14 @@ for i in startsim:days
   #Updating the saturation Layers
     for Lst in 1:Lty
       if(Lst==1)
-         LayersP = LayerUpdate(ddist[Lst,1:NoL],outx[Lst], LayersP, layerUH_P, nodaysvector[Lst,1:NoL], NoL)
+         LayerUpdate!(ddist[Lst,1:NoL],outx[Lst], LayersP, layerUH_P, nodaysvector[Lst,1:NoL], NoL)
     end       
       if(Lst==2)
-        LayersIP = LayerUpdate(ddist[Lst,1:NoL],outx[Lst], LayersIP, layerUH_IP, nodaysvector[Lst,1:NoL], NoL)       
+         LayerUpdate!(ddist[Lst,1:NoL],outx[Lst], LayersIP, layerUH_IP, nodaysvector[Lst,1:NoL], NoL)       
       end
     end  
 
-  BogLayers = BogLayerUpdate(outbog, BogLayers, UHbog, antBogsteps)#
+  BogLayerUpdate!(outbog, BogLayers, UHbog, antBogsteps)#
 
   #summing up groundwater states
     for Lst in 1:Lty

DDDFunctions/DDDUrbanFunc19062023_InfWB.jl

@@ -574,14 +574,14 @@ for i in startsim:days
   #Updating the saturation Layers
     for Lst in 1:Lty
       if(Lst==1)
-         LayersP = LayerUpdate(ddist[Lst,1:NoL],outx[Lst], LayersP, layerUH_P, nodaysvector[Lst,1:NoL], NoL)
+         LayerUpdate!(ddist[Lst,1:NoL],outx[Lst], LayersP, layerUH_P, nodaysvector[Lst,1:NoL], NoL)
     end       
       if(Lst==2)
-        LayersIP = LayerUpdate(ddist[Lst,1:NoL],outx[Lst], LayersIP, layerUH_IP, nodaysvector[Lst,1:NoL], NoL)       
+         LayerUpdate!(ddist[Lst,1:NoL],outx[Lst], LayersIP, layerUH_IP, nodaysvector[Lst,1:NoL], NoL)       
       end
     end  
 
-  BogLayers = BogLayerUpdate(outbog, BogLayers, UHbog, antBogsteps)#
+  BogLayerUpdate!(outbog, BogLayers, UHbog, antBogsteps)#
 
   #summing up groundwater states
     for Lst in 1:Lty

DDDFunctions/LayerCapacityUpdate.jl

@@ -9,16 +9,15 @@
 
 function LayerCapacityUpdate(Layers, nodaysvector, Magkap, NoL)
 
-ddistx = zeros(Float64,NoL)
-aktMag = zeros(Float64,NoL)
+ddistx = Vector{Float64}(undef, NoL)
 
 #Below are the states (in mm) for each saturation level
 for j in reverse(1:NoL)
                                       #state after this timesteps' water is gone. amount of water  in mm, minus current timestep
-  aktMag[j] = sum(Layers[j,2:nodaysvector[j]])
+  aktMag = sum(Layers[j,2:nodaysvector[j]])
 
-  if (aktMag[j] < Magkap[j])
-   ddistx[j] = Magkap[j] - aktMag[j]
+  if (aktMag < Magkap[j])
+   ddistx[j] = Magkap[j] - aktMag
   end
 
 end

DDDFunctions/LayerEstimation.jl

@@ -26,18 +26,14 @@ mLam = GshInt*GscInt
 varLam = GshInt*(GscInt)^2                           #Yevjevich p.145
 meanIntk = mLam*midDL/Timeresinsec                   #mean celerity estimated through Integrated Celerity
 antBox = Int(trunc(maxDl/(meanIntk*Timeresinsec)))+1 #Temporal length UH_MAD
-UH_MAD = zeros(Float64,antBox)
-sRes = zeros(Float64,antBox) # saturation sum
+sRes = Vector{Float64}(undef, antBox) # saturation sum
 
 #Unit hydrograph for MAD
 UH_MAD = SingleUH(meanIntk,Timeresinsec, midDL, maxDl, 0)
 
 StSt = (1000*MAD*Timeresinsec)/(area2)         # Steady state Input eq. output in mm
-sRes[1] = 0
-sRes[2:antBox] .= StSt.*UH_MAD[2:antBox]
-
-for i in 3: antBox
-  sRes[i:antBox] .= sRes[i:antBox] + StSt.*UH_MAD[i:antBox]
+for i in 1:antBox
+  sRes[i] = (i - 1) * StSt * UH_MAD[i]
 end
 
 mRes = sum(sRes)

DDDFunctions/LayerEvap.jl

@@ -15,84 +15,75 @@ function LayerEvap(Layers, nodaysvector, ea_S, layerUH, NoL)
   #ea_S = 4.5
   #layerUH = [0.6 0.4 0.0 0.0 0.0; 0.4 0.35 0.25 0.0 0.0; 0.35 0.25 0.25 0.15 0.0; 0.3 0.25 0.2 0.15 0.1]
 
-  LayersLast = zeros(NoL,nodaysvector[NoL])
-  for i in 1: NoL
-    LayersLast[i,1:nodaysvector[i]].= Layers[i,1:nodaysvector[i]]
-  end
+  total_layers_last = sum(Layers)
+
   #Below are the states (in mm) for each saturation level
   redea = ea_S
-   #println(sum(LayersLast))
+   #println(total_layers_last)
    #println(sum(Layers))
 
 
   for j in 1 : NoL                     # 1 is the top(fastest) Layer NoL is the bottom layer 
 
-    newLayer = zeros(nodaysvector[j])
-
     if(redea > 0.0)
 
       if(sum(Layers) > 0.0)             
-        newLayer .= Layers[j,1:nodaysvector[j]]
+        newLayer = Layers[j,1:nodaysvector[j]]
         aktMag = sum(Layers[j,1:nodaysvector[j]]) # this is correct because ea is a nonintegrated with a continuum variable as opposed to discharge
         differ = aktMag-redea 
         # println(aktMag)
         # println(differ)
         # println(newLayer)
-        # println(sum(LayersLast))
+        # println(total_layers_last)
         # println(sum(Layers))
 
 
         if (differ > 0.0) # the Layer has more water than is to be evaporated > ea_S                
             ea_excess = 0.0
-            evapUH = redea .* layerUH[j,1:nodaysvector[j]]
-            newLayer .= Layers[j,1:nodaysvector[j]] .- evapUH[1:nodaysvector[j]]  
+            @views evapUH = redea .* layerUH[j,1:nodaysvector[j]]
+            @views newLayer = Layers[j,1:nodaysvector[j]] .- evapUH[1:nodaysvector[j]]  
 
             tull = findall(newLayer .< 0.0)
             if(length(tull) > 0)
               x = tull    # locate which boxes have not enough water to evaporate
               ea_excess = sum(evapUH[x]) # the amount which is not evaporated
               evapUH[x] .= 0.0            # the identified boxes have zero instead of negative values
-              newLayer .= Layers[j,1:nodaysvector[j]]- evapUH[1:nodaysvector[j]] # updates so that everthing is positive
+              @views newLayer = Layers[j,1:nodaysvector[j]] .- evapUH[1:nodaysvector[j]] # updates so that everthing is positive
             end   
 
-            if (round(sum(Layers[j,1:nodaysvector[j]])-sum(newLayer)- redea+ea_excess; digits= 8)!= 0.0)        
-              avvik = round(sum(Layers[j,1:nodaysvector[j]]) -sum(newLayer)- redea+ea_excess; digits= 8)
+            if (round(aktMag - sum(newLayer) - redea + ea_excess; digits= 8) != 0.0)        
+              avvik = round(aktMag -sum(newLayer)- redea+ea_excess; digits= 8)
               println("Hei, feil i fordampning", avvik)
             end  
             redea = 0.0
-        end #if differ >0.0
-
-        if (differ <= 0.0)# corresponds to that the Layer looses all water and redea is reduced, layer content < ea
-
-
-          newLayer[1:nodaysvector[j]] .=  0.0
+        else
+          newLayer = zeros(nodaysvector[j])
           redea = redea - aktMag                    
           if(j == NoL)           
               redea = 0.0   # 
           end
           #println(newLayer)
-
         end
-        #println(sum(LayersLast))       
-     Layers[j,1:nodaysvector[j]] .= newLayer[1:nodaysvector[j]]  # updating the actual layer
+        #println(total_layers_last)       
+        Layers[j,1:nodaysvector[j]] .= newLayer  # updating the actual layer
         #println(sum(Layers))
-        #println(sum(LayersLast))    
+        #println(total_layers_last)    
       else                           # sumLayers == 0    
         ea = 0.0
       end                             # sumLayers ==0      
     end                   # redea > 0.0 
   end                          # for loop
 
-  ea = sum(LayersLast)-sum(Layers)
-
+  total_layers_current = sum(Layers)
+  ea = total_layers_last - total_layers_current
 
   if (ea > 0.0)
-     if (round(sum(LayersLast)-(ea+sum(Layers)); digits=8) != 0.0)
+     if (round(total_layers_last - (ea + total_layers_current); digits=8) != 0.0)
        println("Layers ut not OK")
      end
   end
 #println(ea)
-#println(sum(LayersLast))
+#println(total_layers_last)
 #println(sum(Layers))                                   
   return Layers, ea
 end