diff --git a/DDDFunctions/Big2SmallLambda.jl b/DDDFunctions/Big2SmallLambda.jl
index 67006b4..28dd01e 100644
--- a/DDDFunctions/Big2SmallLambda.jl
+++ b/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]
diff --git a/DDDFunctions/BogLayerUpdate.jl b/DDDFunctions/BogLayerUpdate.jl
index c1aaf85..2b33557 100644
--- a/DDDFunctions/BogLayerUpdate.jl
+++ b/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            
 
 
diff --git a/DDDFunctions/DDDAllTerrain22012024.jl b/DDDFunctions/DDDAllTerrain22012024.jl
index 65c8753..96abbbe 100644
--- a/DDDFunctions/DDDAllTerrain22012024.jl
+++ b/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 
diff --git a/DDDFunctions/DDDEventFloodDesign.jl b/DDDFunctions/DDDEventFloodDesign.jl
index 785d748..934b89c 100644
--- a/DDDFunctions/DDDEventFloodDesign.jl
+++ b/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  
    
diff --git a/DDDFunctions/DDDEventFloodDesignFixedPDynSM.jl b/DDDFunctions/DDDEventFloodDesignFixedPDynSM.jl
index 6f782db..6e3b8e0 100644
--- a/DDDFunctions/DDDEventFloodDesignFixedPDynSM.jl
+++ b/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  
    
diff --git a/DDDFunctions/DDDModulFunc-MAD-LAKE.jl b/DDDFunctions/DDDModulFunc-MAD-LAKE.jl
index d2efe6f..182b6bc 100644
--- a/DDDFunctions/DDDModulFunc-MAD-LAKE.jl
+++ b/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
   
diff --git a/DDDFunctions/DDDUrbanFunc10012024_InfWB.jl b/DDDFunctions/DDDUrbanFunc10012024_InfWB.jl
index c49bd82..2ed820d 100644
--- a/DDDFunctions/DDDUrbanFunc10012024_InfWB.jl
+++ b/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
diff --git a/DDDFunctions/DDDUrbanFunc19062023_InfWB.jl b/DDDFunctions/DDDUrbanFunc19062023_InfWB.jl
index b62703c..2e23b13 100644
--- a/DDDFunctions/DDDUrbanFunc19062023_InfWB.jl
+++ b/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
diff --git a/DDDFunctions/LayerCapacityUpdate.jl b/DDDFunctions/LayerCapacityUpdate.jl
index 4b6c462..58a8425 100644
--- a/DDDFunctions/LayerCapacityUpdate.jl
+++ b/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
diff --git a/DDDFunctions/LayerEstimation.jl b/DDDFunctions/LayerEstimation.jl
index 9909e80..9a0fd44 100644
--- a/DDDFunctions/LayerEstimation.jl
+++ b/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)
diff --git a/DDDFunctions/LayerEvap.jl b/DDDFunctions/LayerEvap.jl
index e35275d..8991b2b 100644
--- a/DDDFunctions/LayerEvap.jl
+++ b/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
diff --git a/DDDFunctions/LayerUpdate.jl b/DDDFunctions/LayerUpdate.jl
index 0172d79..c91afff 100644
--- a/DDDFunctions/LayerUpdate.jl
+++ b/DDDFunctions/LayerUpdate.jl
@@ -7,28 +7,13 @@
 #     Revised: 17.12.2019
 #--------------------------------------------------------------------------
 
-function LayerUpdate(ddist, outx, Layers, layerUH, nodaysvector, NoL)
+function LayerUpdate!(ddist::Vector{Float64}, outx::Float64, Layers::Matrix{Float64}, layerUH::Matrix{Float64}, nodaysvector::Vector{Int}, NoL::Int)
 
-# ddist: 1 dim array float
-# Layers: 2 dim array float
-# layersUH: 2 dim array float
-# NoL: scalar, integer    
-# nodaysvector: 1 dim array integer
-# outx: scalar, float
-        
-  for j in 1 : NoL
-    qlayer = zeros(nodaysvector[j])    
-    qlayer .= ddist[j]*outx .* layerUH[j,1:nodaysvector[j]]    #finds response  in mm!!!!  for the actual layer, en vektor     
-        
-    if(nodaysvector[j] > 1)
-      Layers[j,(1:(nodaysvector[j]-1))] .= Layers[j,2:nodaysvector[j]] .+ qlayer[2:nodaysvector[j]]# flytter the level of the matrix one timestep ahead
-      Layers[j,nodaysvector[j]] = 0.0       
-    end 
-
-    if(nodaysvector[j] == 1)
-      Layers[j,1:nodaysvector[j]] .= qlayer[1:nodaysvector[j]]
+  for j in 1:NoL
+    multiplier = ddist[j] * outx
+    Layers[j,1] = multiplier * layerUH[j,1]
+    for h in 2:nodaysvector[j]
+      Layers[j,h-1] = Layers[j,h] + multiplier * layerUH[j,h]
     end
   end
-
-return Layers
 end                
diff --git a/DDDFunctions/OverlandFlowDynamicDD.jl b/DDDFunctions/OverlandFlowDynamicDD.jl
index 7aa514c..72c461d 100644
--- a/DDDFunctions/OverlandFlowDynamicDD.jl
+++ b/DDDFunctions/OverlandFlowDynamicDD.jl
@@ -58,7 +58,7 @@ function OverlandFlowDynamicDD(k,ddist,outx, layerUH, nodaysvector,NoL, midDL, C
 
       OFUH = SingleUH(celerityOF,Timeresinsec, dmean, maxdistOF, zOF)
       if (length(OFUH) < nodaysvector[1])             # ensures DNR is more dense than original RN
-         layerUH[1,1:length(OFUH)].= OFUH[1:length(OFUH)]
+         layerUH[1,1:length(OFUH)] .= OFUH
          layerUH[1,((length(OFUH)+1):nodaysvector[NoL])] .= 0.0
       end
   end
diff --git a/DDDFunctions/PyrAreas.jl b/DDDFunctions/PyrAreas.jl
index e0e9c19..04586eb 100644
--- a/DDDFunctions/PyrAreas.jl
+++ b/DDDFunctions/PyrAreas.jl
@@ -32,8 +32,7 @@ function PyrAreas(NoL,totarea,maxDl,nodaysvector, layerUH, antHorlag)
        Areas[i,j] = totarea.*layerUH[i,j] 
     end
   end
-  delta_d = zeros(NoL)                            # in meters (height as in the pyramid plots) pr. time-step box
-  delta_d[1:NoL] = maxDl ./nodaysvector[1:NoL]  
+  delta_d = maxDl ./ nodaysvector  # in meters (height as in the pyramid plots) pr. time-step box
   
 
 #To be used in GRW_point subroutne: Layers__mm <- Layers*totarea/areas ######################
diff --git a/DDDFunctions/SingleUH.jl b/DDDFunctions/SingleUH.jl
index 69754c8..0ed759c 100644
--- a/DDDFunctions/SingleUH.jl
+++ b/DDDFunctions/SingleUH.jl
@@ -10,20 +10,22 @@ function SingleUH(k,Timeresinsec,meanD,maxD, nugget)
 #Timeresinsec: scalar Integer = 86400
 #k: scalar float, NB enters as an array float in other subroutines
     
-ant = trunc(maxD/(k*Timeresinsec))+1
-ant= Int(ant)    
-UHvec = zeros(Float64,ant)
+ant = Int(trunc(maxD/(k*Timeresinsec))) + 1
+UHvec = Vector{Float64}(undef, ant)
 
 escl =  (meanD/k)/Timeresinsec 
 Uhexp = Exponential(escl)    
-UHvec[1] = (nugget+(1-nugget)*cdf(Uhexp,1)) #1/escl is parameter of distribution, escl is the mean of timesteps needed to drain.
+new = (1 - nugget)*cdf(Uhexp, 1) # 1/escl is parameter of distribution, escl is the mean of timesteps needed to drain.
+UHvec[1] = nugget + new
 if(ant > 1)
    for i in 2:ant
-    UHvec[i] = (nugget+(1-nugget)*cdf(Uhexp,i))-(nugget+(1-nugget)*cdf(Uhexp,(i-1))) #Makes exponential pdf.  
+      old = new
+      new = (1 - nugget)*cdf(Uhexp, i)
+      UHvec[i] = new - old #Makes exponential pdf.  
                                            # UH is 1/e1ascl*exp(-(1/e1ascl)*t)/sum_for_all_t(1/e1ascl*exp(-(1/e1ascl)*t))
    end
 end
-UHvec[1:ant] = UHvec[1:ant]/sum(UHvec)                       #Must normalise so that the sum equals 1
+UHvec /= sum(UHvec)                       #Must normalise so that the sum equals 1
 
 return UHvec                                                 #Unit hydrograph, i.e. weights distributing runoff in time
 end
diff --git a/DDDFunctions/SnowpackTemp.jl b/DDDFunctions/SnowpackTemp.jl
index 4f444c7..776cab0 100644
--- a/DDDFunctions/SnowpackTemp.jl
+++ b/DDDFunctions/SnowpackTemp.jl
@@ -1,27 +1,13 @@
 function SnowpackTemp(Temp)
-
 #This function calculates a weighted (linearly decreasing) mean of the air temperature as an estimate of the skin tmeperature
 # of the snowpack. It is assumed that the temperature at ground is zero such that the mean temperature is estT/2
-    
 # Temp: : 1 dim array float
-    
+# The temperature vector is organzed such that  the first element is todays temp
     len = length(Temp)
-
-    wgt = zeros(Float64,len)
-                               #The temperature vector is organzed such that  the first element is todays temp
-    for i in 1:len
-      wgt[i] = (len-i+1)/len
-    end
-
-    summen = sum(wgt)
-
-    for i in 1:len 
-        wgt[i] = wgt[i]/summen
-    end
-     snittT = sum(wgt.*Temp)          # This is the snowpack temperature. Tss is snittT *2
+    wgt = collect(len:-1:1) / (len*(len+1)/2)
+    snittT = sum(wgt.*Temp)          # This is the snowpack temperature. Tss is snittT *2
     if snittT > 0.0 
         snittT = 0.0 # Snowpack temperature cannot be more than zero
-    end      
-    
+    end
 return snittT
 end
diff --git a/DDDFunctions/TempstartUpdate.jl b/DDDFunctions/TempstartUpdate.jl
index 7d35c73..8db076f 100644
--- a/DDDFunctions/TempstartUpdate.jl
+++ b/DDDFunctions/TempstartUpdate.jl
@@ -7,16 +7,9 @@
 #     Revised: 19.02.2024
 #--------------------------------------------------------------------------
 
-function TempstartUpdate(tempmatrix, temps, len)
-
+function TempstartUpdate!(tempmatrix::Matrix{Float64}, temps::Vector{Float64}, len::Int)
 # tempmatrix: timeseries of length 5 days 8length dep on temporal resolution) for 10 elevation zones
 # temps: this timesteps temperatures
-
-  dummy = tempmatrix      
-
-      dummy[1:(len-1), 1:10] .= dummy[2:len, 1:10]# flytter the level of the matrix one timestep ahead
-      dummy[len, 1:10] = temps[1:10]
-      tempmatrix = dummy 
-
-return tempmatrix
+      @views tempmatrix[1:len-1,:] .= tempmatrix[2:len,:] # move the level of the matrix one timestep ahead
+      tempmatrix[len,:] .= temps
 end                
diff --git a/DDDFunctions/Varc.jl b/DDDFunctions/Varc.jl
index 5119e32..cb6e1ec 100644
--- a/DDDFunctions/Varc.jl
+++ b/DDDFunctions/Varc.jl
@@ -90,7 +90,7 @@ end #for u > 0.0 END ABLATION
        nudyn = (meantot^2/vartot)
        alphadyn = meantot/vartot  # Conditional distribution
       end
-      if(vartot =="NaN")
+      if isnan(vartot)
        println("vartot=NaN (ny_varc) nudyn= ",nudyn," alphadyn= ",alphadyn, " meantot= ", meantot," vartot= ",vartot," nnn= ",nnn," n= ",n," u= ",u)
        sleep(5) #stops for 5 seconds
       end
diff --git a/Manifest.toml b/Manifest.toml
index 9ba9315..af51e21 100644
--- a/Manifest.toml
+++ b/Manifest.toml
@@ -2,7 +2,7 @@
 
 julia_version = "1.10.10"
 manifest_format = "2.0"
-project_hash = "2e9983a1e4d4bad7f0d6f34959fca9cc5d523b5a"
+project_hash = "bc40fddbefe16d7163d0816bfb3ac92f722eafa3"
 
 [[deps.ADTypes]]
 git-tree-sha1 = "f7304359109c768cf32dc5fa2d371565bb63b68a"
@@ -19,6 +19,11 @@ version = "1.21.0"
     ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
     EnzymeCore = "f151be2c-9106-41f4-ab19-57ee4f262869"
 
+[[deps.AbstractTrees]]
+git-tree-sha1 = "2d9c9a55f9c93e8887ad391fbae72f8ef55e1177"
+uuid = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
+version = "0.4.5"
+
 [[deps.Adapt]]
 deps = ["LinearAlgebra", "Requires"]
 git-tree-sha1 = "7e35fca2bdfba44d797c53dfe63a51fabf39bfc0"
@@ -110,6 +115,11 @@ git-tree-sha1 = "1b96ea4a01afe0ea4090c5c8039690672dd13f2e"
 uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0"
 version = "1.0.9+0"
 
+[[deps.CEnum]]
+git-tree-sha1 = "389ad5c84de1ae7cf0e28e381131c98ea87d54fc"
+uuid = "fa961155-64e5-5f13-b03f-caf6b980ea82"
+version = "0.5.0"
+
 [[deps.CPUTime]]
 git-tree-sha1 = "2dcc50ea6a0a1ef6440d6eecd0fe3813e5671f45"
 uuid = "a9c8d775-2e2e-55fc-8582-045d282d599e"
@@ -121,6 +131,12 @@ git-tree-sha1 = "deddd8725e5e1cc49ee205a1964256043720a6c3"
 uuid = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 version = "0.10.15"
 
+[[deps.Cairo]]
+deps = ["Cairo_jll", "Colors", "Glib_jll", "Graphics", "Libdl", "Pango_jll"]
+git-tree-sha1 = "71aa551c5c33f1a4415867fe06b7844faadb0ae9"
+uuid = "159f3aea-2a34-519c-b102-8c37f9878175"
+version = "1.1.1"
+
 [[deps.Cairo_jll]]
 deps = ["Artifacts", "Bzip2_jll", "CompilerSupportLibraries_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "JLLWrappers", "LZO_jll", "Libdl", "Pixman_jll", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Zlib_jll", "libpng_jll"]
 git-tree-sha1 = "fde3bf89aead2e723284a8ff9cdf5b551ed700e8"
@@ -248,6 +264,18 @@ git-tree-sha1 = "249fe38abf76d48563e2f4556bebd215aa317e15"
 uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
 version = "4.1.1"
 
+[[deps.Cthulhu]]
+deps = ["CodeTracking", "FoldingTrees", "InteractiveUtils", "JuliaSyntax", "PrecompileTools", "Preferences", "REPL", "TypedSyntax", "UUIDs", "Unicode", "WidthLimitedIO"]
+git-tree-sha1 = "af124ef442564dd4d87dcbf67fb1d2f9fc8f7297"
+uuid = "f68482b8-f384-11e8-15f7-abe071a5a75f"
+version = "2.16.7"
+
+    [deps.Cthulhu.extensions]
+    CthulhuCompilerExt = "Compiler"
+
+    [deps.Cthulhu.weakdeps]
+    Compiler = "807dbc54-b67e-4c79-8afb-eafe4df6f2e1"
+
 [[deps.DataAPI]]
 git-tree-sha1 = "abe83f3a2f1b857aac70ef8b269080af17764bbe"
 uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
@@ -479,6 +507,18 @@ git-tree-sha1 = "05882d6995ae5c12bb5f36dd2ed3f61c98cbb172"
 uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
 version = "0.8.5"
 
+[[deps.FlameGraphs]]
+deps = ["AbstractTrees", "Colors", "FileIO", "FixedPointNumbers", "IndirectArrays", "LeftChildRightSiblingTrees", "Profile"]
+git-tree-sha1 = "0166baf81babb91cf78bfcc771d8e87c43d568df"
+uuid = "08572546-2f56-4bcf-ba4e-bab62c3a3f89"
+version = "1.1.0"
+
+[[deps.FoldingTrees]]
+deps = ["AbstractTrees", "REPL"]
+git-tree-sha1 = "c1b0164369256b26f71d9830df9000a9c39757fc"
+uuid = "1eca21be-9b9b-4ed8-839a-6d8ae26b1781"
+version = "1.2.2"
+
 [[deps.Fontconfig_jll]]
 deps = ["Artifacts", "Bzip2_jll", "Expat_jll", "FreeType2_jll", "JLLWrappers", "Libdl", "Libuuid_jll", "Zlib_jll"]
 git-tree-sha1 = "f85dac9a96a01087df6e3a749840015a0ca3817d"
@@ -542,6 +582,12 @@ git-tree-sha1 = "7dd7173f7129a1b6f84e0f03e0890cd1189b0659"
 uuid = "d2c73de3-f751-5644-a686-071e5b155ba9"
 version = "0.73.22+0"
 
+[[deps.GTK4_jll]]
+deps = ["Artifacts", "Cairo_jll", "Fontconfig_jll", "FreeType2_jll", "FriBidi_jll", "GettextRuntime_jll", "Glib_jll", "Graphene_jll", "HarfBuzz_jll", "JLLWrappers", "Libdl", "Libepoxy_jll", "Libtiff_jll", "PCRE2_jll", "Pango_jll", "Wayland_jll", "Wayland_protocols_jll", "Xorg_libX11_jll", "Xorg_libXcursor_jll", "Xorg_libXdamage_jll", "Xorg_libXext_jll", "Xorg_libXfixes_jll", "Xorg_libXi_jll", "Xorg_libXinerama_jll", "Xorg_libXrandr_jll", "Xorg_libXrender_jll", "gdk_pixbuf_jll", "iso_codes_jll", "xkbcommon_jll"]
+git-tree-sha1 = "41c6d2ab0a034c35c6ba2c7d509d3a03da64fc69"
+uuid = "6ebb71f1-8434-552f-b6b1-dc18babcca63"
+version = "4.18.6+0"
+
 [[deps.GettextRuntime_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Libiconv_jll"]
 git-tree-sha1 = "45288942190db7c5f760f59c04495064eedf9340"
@@ -560,6 +606,18 @@ git-tree-sha1 = "24f6def62397474a297bfcec22384101609142ed"
 uuid = "7746bdde-850d-59dc-9ae8-88ece973131d"
 version = "2.86.3+0"
 
+[[deps.Graphene_jll]]
+deps = ["Artifacts", "Glib_jll", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "47a6897c1ce13a484b5f62643ed8f84e02babb05"
+uuid = "75302f13-0b7e-5bab-a6d1-23fa92e4c2ea"
+version = "1.10.8+0"
+
+[[deps.Graphics]]
+deps = ["Colors", "LinearAlgebra", "NaNMath"]
+git-tree-sha1 = "a641238db938fff9b2f60d08ed9030387daf428c"
+uuid = "a2bd30eb-e257-5431-a919-1863eab51364"
+version = "1.1.3"
+
 [[deps.Graphite2_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
 git-tree-sha1 = "8a6dbda1fd736d60cc477d99f2e7a042acfa46e8"
@@ -571,6 +629,18 @@ git-tree-sha1 = "53bb909d1151e57e2484c3d1b53e19552b887fb2"
 uuid = "42e2da0e-8278-4e71-bc24-59509adca0fe"
 version = "1.0.2"
 
+[[deps.Gtk4]]
+deps = ["BitFlags", "CEnum", "Cairo", "Cairo_jll", "ColorTypes", "Dates", "FixedPointNumbers", "GTK4_jll", "Glib_jll", "Graphene_jll", "Graphics", "JLLWrappers", "Libdl", "Librsvg_jll", "Pango_jll", "Preferences", "Reexport", "Scratch", "Xorg_xkeyboard_config_jll", "adwaita_icon_theme_jll", "gdk_pixbuf_jll", "hicolor_icon_theme_jll", "libpng_jll"]
+git-tree-sha1 = "10578b975a68cbe862390e7219b1f6da4b468766"
+uuid = "9db2cae5-386f-4011-9d63-a5602296539b"
+version = "0.7.12"
+
+[[deps.GtkObservables]]
+deps = ["Cairo", "Colors", "Dates", "FixedPointNumbers", "Graphics", "Gtk4", "IntervalSets", "LinearAlgebra", "Observables", "PrecompileTools", "Reexport", "RoundingIntegers"]
+git-tree-sha1 = "7cd1331aece756ab5b76c9e401996faad42ea654"
+uuid = "8710efd8-4ad6-11eb-33ea-2d5ceb25a41c"
+version = "2.2.2"
+
 [[deps.HTTP]]
 deps = ["Base64", "CodecZlib", "ConcurrentUtilities", "Dates", "ExceptionUnwrapping", "Logging", "LoggingExtras", "MbedTLS", "NetworkOptions", "OpenSSL", "PrecompileTools", "Random", "SimpleBufferStream", "Sockets", "URIs", "UUIDs"]
 git-tree-sha1 = "5e6fe50ae7f23d171f44e311c2960294aaa0beb5"
@@ -594,6 +664,11 @@ git-tree-sha1 = "68c173f4f449de5b438ee67ed0c9c748dc31a2ec"
 uuid = "34004b35-14d8-5ef3-9330-4cdb6864b03a"
 version = "0.3.28"
 
+[[deps.IndirectArrays]]
+git-tree-sha1 = "012e604e1c7458645cb8b436f8fba789a51b257f"
+uuid = "9b13fd28-a010-5f03-acff-a1bbcff69959"
+version = "1.0.0"
+
 [[deps.Infiltrator]]
 deps = ["InteractiveUtils", "Markdown", "REPL", "UUIDs"]
 git-tree-sha1 = "d07c6718c1ce7324341947f8f608f5d69795dde5"
@@ -617,6 +692,17 @@ version = "1.4.5"
 deps = ["Markdown"]
 uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 
+[[deps.IntervalSets]]
+git-tree-sha1 = "d966f85b3b7a8e49d034d27a189e9a4874b4391a"
+uuid = "8197267c-284f-5f27-9208-e0e47529a953"
+version = "0.7.13"
+weakdeps = ["Random", "RecipesBase", "Statistics"]
+
+    [deps.IntervalSets.extensions]
+    IntervalSetsRandomExt = "Random"
+    IntervalSetsRecipesBaseExt = "RecipesBase"
+    IntervalSetsStatisticsExt = "Statistics"
+
 [[deps.InvertedIndices]]
 git-tree-sha1 = "6da3c4316095de0f5ee2ebd875df8721e7e0bdbe"
 uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f"
@@ -740,6 +826,12 @@ version = "0.16.10"
     SymEngine = "123dc426-2d89-5057-bbad-38513e3affd8"
     tectonic_jll = "d7dd28d6-a5e6-559c-9131-7eb760cdacc5"
 
+[[deps.LeftChildRightSiblingTrees]]
+deps = ["AbstractTrees"]
+git-tree-sha1 = "95ba48564903b43b2462318aa243ee79d81135ff"
+uuid = "1d6d02ad-be62-4b6b-8a6d-2f90e265016e"
+version = "0.2.1"
+
 [[deps.LibCURL]]
 deps = ["LibCURL_jll", "MozillaCACerts_jll"]
 uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"
@@ -767,6 +859,12 @@ version = "1.11.0+1"
 [[deps.Libdl]]
 uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
 
+[[deps.Libepoxy_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Libglvnd_jll", "Xorg_libX11_jll"]
+git-tree-sha1 = "94d068f57b4241dd090693b6aba63416892298de"
+uuid = "42c93a91-0102-5b3f-8f9d-e41de60ac950"
+version = "1.5.11+0"
+
 [[deps.Libffi_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
 git-tree-sha1 = "c8da7e6a91781c41a863611c7e966098d783c57a"
@@ -791,6 +889,12 @@ git-tree-sha1 = "97bbca976196f2a1eb9607131cb108c69ec3f8a6"
 uuid = "4b2f31a3-9ecc-558c-b454-b3730dcb73e9"
 version = "2.41.3+0"
 
+[[deps.Librsvg_jll]]
+deps = ["Artifacts", "Cairo_jll", "FreeType2_jll", "Glib_jll", "JLLWrappers", "Libdl", "Pango_jll", "XML2_jll", "gdk_pixbuf_jll"]
+git-tree-sha1 = "e6ab5dda9916d7041356371c53cdc00b39841c31"
+uuid = "925c91fb-5dd6-59dd-8e8c-345e74382d89"
+version = "2.54.7+0"
+
 [[deps.Libtiff_jll]]
 deps = ["Artifacts", "JLLWrappers", "JpegTurbo_jll", "LERC_jll", "Libdl", "XZ_jll", "Zlib_jll", "Zstd_jll"]
 git-tree-sha1 = "f04133fe05eff1667d2054c53d59f9122383fe05"
@@ -869,6 +973,12 @@ git-tree-sha1 = "b513cedd20d9c914783d8ad83d08120702bf2c77"
 uuid = "442fdcdd-2543-5da2-b0f3-8c86c306513e"
 version = "0.3.3"
 
+[[deps.MethodAnalysis]]
+deps = ["AbstractTrees"]
+git-tree-sha1 = "39277c4f6a771f45613e552a9ae826be1bf873fe"
+uuid = "85b6ec6f-f7df-4429-9514-a64bcd9ee824"
+version = "1.0.0"
+
 [[deps.Missings]]
 deps = ["DataAPI"]
 git-tree-sha1 = "ec4f7fbeab05d7747bdf98eb74d130a2a2ed298d"
@@ -898,6 +1008,11 @@ version = "1.1.3"
 uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
 version = "1.2.0"
 
+[[deps.Observables]]
+git-tree-sha1 = "7438a59546cf62428fc9d1bc94729146d37a7225"
+uuid = "510215fc-4207-5dde-b226-833fc4488ee2"
+version = "0.5.5"
+
 [[deps.Ogg_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
 git-tree-sha1 = "b6aa4566bb7ae78498a5e68943863fa8b5231b59"
@@ -943,6 +1058,12 @@ git-tree-sha1 = "05868e21324cede2207c6f0f466b4bfef6d5e7ee"
 uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 version = "1.8.1"
 
+[[deps.OwnTime]]
+deps = ["Printf", "Profile"]
+git-tree-sha1 = "00d9140789be6f702ef5846fc9b9e34a62cb8c8b"
+uuid = "18732c20-e27e-497f-aa49-3bf01a8fc721"
+version = "0.1.0"
+
 [[deps.PCRE2_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "efcefdf7-47ab-520b-bdef-62a2eaa19f15"
@@ -1045,6 +1166,16 @@ uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 deps = ["Printf"]
 uuid = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79"
 
+[[deps.ProfileView]]
+deps = ["Cairo", "Colors", "Dates", "FileIO", "FlameGraphs", "Graphics", "Gtk4", "GtkObservables", "InteractiveUtils", "IntervalSets", "MethodAnalysis", "PrecompileTools", "Preferences", "Profile", "UUIDs"]
+git-tree-sha1 = "47125aaefe988290d6384bdad2c060fd9ca3f10e"
+uuid = "c46f51b8-102a-5cf2-8d2c-8597cb0e0da7"
+version = "1.10.3"
+weakdeps = ["Cthulhu"]
+
+    [deps.ProfileView.extensions]
+    ProfileViewCthulhuExt = "Cthulhu"
+
 [[deps.PtrArrays]]
 git-tree-sha1 = "1d36ef11a9aaf1e8b74dacc6a731dd1de8fd493d"
 uuid = "43287f4e-b6f4-7ad1-bb20-aadabca52c3d"
@@ -1135,6 +1266,11 @@ git-tree-sha1 = "58cdd8fb2201a6267e1db87ff148dd6c1dbd8ad8"
 uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f"
 version = "0.5.1+0"
 
+[[deps.RoundingIntegers]]
+git-tree-sha1 = "99acd97f396ea71a5be06ba6de5c9defe188a778"
+uuid = "d5f540fe-1c90-5db3-b776-2e2f362d9394"
+version = "1.1.0"
+
 [[deps.SHA]]
 uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
 version = "0.7.0"
@@ -1302,6 +1438,12 @@ git-tree-sha1 = "0c45878dcfdcfa8480052b6ab162cdd138781742"
 uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
 version = "0.11.3"
 
+[[deps.TypedSyntax]]
+deps = ["CodeTracking", "JuliaSyntax"]
+git-tree-sha1 = "9b3d96dc9de6e356e56d6b80f46187fbdd6ee23b"
+uuid = "d265eb64-f81a-44ad-a842-4247ee1503de"
+version = "1.4.4"
+
 [[deps.URIs]]
 git-tree-sha1 = "bef26fb046d031353ef97a82e3fdb6afe7f21b1a"
 uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4"
@@ -1337,17 +1479,35 @@ git-tree-sha1 = "96478df35bbc2f3e1e791bc7a3d0eeee559e60e9"
 uuid = "a2964d1f-97da-50d4-b82a-358c7fce9d89"
 version = "1.24.0+0"
 
+[[deps.Wayland_protocols_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "54b8a029ac145ebe8299463447fd1590b2b1d92f"
+uuid = "2381bf8a-dfd0-557d-9999-79630e7b1b91"
+version = "1.44.0+0"
+
 [[deps.WeakRefStrings]]
 deps = ["DataAPI", "InlineStrings", "Parsers"]
 git-tree-sha1 = "b1be2855ed9ed8eac54e5caff2afcdb442d52c23"
 uuid = "ea10d353-3f73-51f8-a26c-33c1cb351aa5"
 version = "1.4.2"
 
+[[deps.WidthLimitedIO]]
+deps = ["Unicode"]
+git-tree-sha1 = "71142739e695823729a335e9bc124ef41ec1433a"
+uuid = "b8c1c048-cf81-46c6-9da0-18c1d99e41f2"
+version = "1.0.1"
+
 [[deps.WorkerUtilities]]
 git-tree-sha1 = "cd1659ba0d57b71a464a29e64dbc67cfe83d54e7"
 uuid = "76eceee3-57b5-4d4a-8e66-0e911cebbf60"
 version = "1.6.1"
 
+[[deps.XML2_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Zlib_jll"]
+git-tree-sha1 = "80d3930c6347cfce7ccf96bd3bafdf079d9c0390"
+uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a"
+version = "2.13.9+0"
+
 [[deps.XZ_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
 git-tree-sha1 = "9cce64c0fdd1960b597ba7ecda2950b5ed957438"
@@ -1384,6 +1544,12 @@ git-tree-sha1 = "6c74ca84bbabc18c4547014765d194ff0b4dc9da"
 uuid = "935fb764-8cf2-53bf-bb30-45bb1f8bf724"
 version = "1.2.4+0"
 
+[[deps.Xorg_libXdamage_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libXfixes_jll"]
+git-tree-sha1 = "45917154defafcbff585d6c04f323b5797da1ecf"
+uuid = "0aeada51-83db-5f97-b67e-184615cfc6f6"
+version = "1.1.7+0"
+
 [[deps.Xorg_libXdmcp_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
 git-tree-sha1 = "52858d64353db33a56e13c341d7bf44cd0d7b309"
@@ -1503,6 +1669,12 @@ git-tree-sha1 = "446b23e73536f84e8037f5dce465e92275f6a308"
 uuid = "3161d3a3-bdf6-5164-811a-617609db77b4"
 version = "1.5.7+1"
 
+[[deps.adwaita_icon_theme_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "hicolor_icon_theme_jll"]
+git-tree-sha1 = "28401767f30e5743ef5e3b0be71417bc911d3952"
+uuid = "b437f822-2cd6-5e08-a15c-8bac984d38ee"
+version = "43.0.1+0"
+
 [[deps.eudev_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
 git-tree-sha1 = "c3b0e6196d50eab0c5ed34021aaa0bb463489510"
@@ -1515,6 +1687,24 @@ git-tree-sha1 = "b6a34e0e0960190ac2a4363a1bd003504772d631"
 uuid = "214eeab7-80f7-51ab-84ad-2988db7cef09"
 version = "0.61.1+0"
 
+[[deps.gdk_pixbuf_jll]]
+deps = ["Artifacts", "Glib_jll", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Libtiff_jll", "Xorg_libX11_jll", "libpng_jll"]
+git-tree-sha1 = "895f21b699121d1a57ecac57e65a852caf569254"
+uuid = "da03df04-f53b-5353-a52f-6a8b0620ced0"
+version = "2.42.13+0"
+
+[[deps.hicolor_icon_theme_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "8f7bcbb45ea61c1ab3d5508413ff5d69565e3ff1"
+uuid = "059c91fe-1bad-52ad-bddd-f7b78713c282"
+version = "0.18.0+0"
+
+[[deps.iso_codes_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "9200917315ff176d445ccc759945ea3fdb7c49c3"
+uuid = "bf975903-5238-5d20-8243-bc370bc1e7e5"
+version = "4.18.0+0"
+
 [[deps.libaom_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
 git-tree-sha1 = "371cc681c00a3ccc3fbc5c0fb91f58ba9bec1ecf"
diff --git a/Project.toml b/Project.toml
index 1a071ba..4f4e101 100644
--- a/Project.toml
+++ b/Project.toml
@@ -3,6 +3,7 @@ BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 BlackBoxOptim = "a134a8b2-14d6-55f6-9291-3336d3ab0209"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 Configurations = "5218b696-f38b-4ac9-8b61-a12ec717816d"
+Cthulhu = "f68482b8-f384-11e8-15f7-abe071a5a75f"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
@@ -10,6 +11,9 @@ Infiltrator = "5903a43b-9cc3-4c30-8d17-598619ec4e9b"
 JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 LsqFit = "2fda8390-95c7-5789-9bda-21331edee243"
+OwnTime = "18732c20-e27e-497f-aa49-3bf01a8fc721"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+Profile = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79"
+ProfileView = "c46f51b8-102a-5cf2-8d2c-8597cb0e0da7"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
diff --git a/calibration/run_single_catchment.jl b/calibration/run_single_catchment.jl
index 9c793df..eaa9624 100644
--- a/calibration/run_single_catchment.jl
+++ b/calibration/run_single_catchment.jl
@@ -1,5 +1,5 @@
 # To run this script:
-# julia --project=. run_single_catchment.jl settings/calibration_251120.toml 2.11 calibration
+# julia --project=.. run_single_catchment.jl settings/calibration_251120.toml 2.11 calibration
 include("calibration.jl")
 
 @assert length(ARGS) == 3
@@ -22,3 +22,10 @@ using BenchmarkTools
 #using ProfileView
 #using Cthulhu
 #@profview runSingleCatchment("settings/calibration_251120.toml", "2.11", "calibration")
+
+# REPL
+#using Profile, OwnTime
+#include("calibration.jl")
+#Profile.clear()
+#@profile runSingleCatchment("settings/calibration_251120.toml", "2.11", "calibration")
+#owntime()