Add basic operators

Project.toml

@@ -19,6 +19,7 @@ OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 PoissonRandom = "e409e4f3-bfea-5376-8464-e040bb5c01ab"
 PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"

src/Forcing.jl

@@ -58,8 +58,8 @@ function update(self::ForcingBase{ForcingStandard}, GMV::GridMeanVariables)
     if self.apply_subsidence
         @inbounds for k in real_center_indicies(self.Gr)
             # Apply large-scale subsidence tendencies
-            GMV.H.tendencies[k] -= (GMV.H.values[k + 1] - GMV.H.values[k]) * self.Gr.dzi * self.subsidence[k]
-            GMV.QT.tendencies[k] -= (GMV.QT.values[k + 1] - GMV.QT.values[k]) * self.Gr.dzi * self.subsidence[k]
+            GMV.H.tendencies[k] -= ∇_upwind(GMV.H.values, self.Gr, k) * self.subsidence[k]
+            GMV.QT.tendencies[k] -= ∇_upwind(GMV.QT.values, self.Gr, k) * self.subsidence[k]
         end
     end
     if self.apply_coriolis
@@ -85,8 +85,8 @@ function update(self::ForcingBase{ForcingDYCOMS_RF01}, GMV::GridMeanVariables)
     @inbounds for k in real_center_indicies(self.Gr)
         GMV.QT.tendencies[k] += self.dqtdt[k]
         # Apply large-scale subsidence tendencies
-        GMV.H.tendencies[k] -= (GMV.H.values[k + 1] - GMV.H.values[k]) * self.Gr.dzi * self.subsidence[k]
-        GMV.QT.tendencies[k] -= (GMV.QT.values[k + 1] - GMV.QT.values[k]) * self.Gr.dzi * self.subsidence[k]
+        GMV.H.tendencies[k] -= ∇_upwind(GMV.H.values, self.Gr, k) * self.subsidence[k]
+        GMV.QT.tendencies[k] -= ∇_upwind(GMV.QT.values, self.Gr, k) * self.subsidence[k]
     end
 
     if self.apply_coriolis

src/Operators.jl

@@ -0,0 +1,35 @@
+
+
+function ∇_collocated(f, grid::Grid, k::Int)
+    return (f[k+1] - f[k-1]) / 2 * grid.dzi
+end
+
+function ∇_collocated(f::SVector, grid::Grid)
+    @assert length(f) == 3
+    return (∇_staggered(f[2:3], grid) - ∇_staggered(f[1:2], grid)) * grid.dzi
+end
+
+function ∇_staggered(f, grid::Grid)
+    @assert length(f) == 2
+    return (f[2] - f[1]) * grid.dzi
+end
+
+function ∇f2c(f, grid::Grid, k::Int)
+    return (f[k] - f[k-1]) * grid.dzi
+end
+
+function ∇c2f(f, grid::Grid, k::Int)
+    return (f[k+1] - f[k]) * grid.dzi
+end
+
+function ∇_upwind(f, grid::Grid, k::Int)
+    return (f[k+1] - f[k]) * grid.dzi
+end
+
+function interpc2f(f, grid::Grid, k::Int)
+    return 0.5*(f[k+1] + f[k])
+end
+
+function cut(f, k)
+    return SVector(f[k-1], f[k], f[k+1])
+end

src/Turbulence.jl

@@ -27,7 +27,7 @@ end
 # Update the diagnosis of the inversion height, using the maximum temperature gradient method
 function update_inversion(self::ParameterizationBase, GMV::GridMeanVariables, option)
     theta_rho = center_field(self.Gr)
-    maxgrad = 0.0
+    ∇θ_liq_max = 0.0
     k_fi = first_center(self.Gr)
 
     @inbounds for k in real_center_indicies(self.Gr)
@@ -46,9 +46,9 @@ function update_inversion(self::ParameterizationBase, GMV::GridMeanVariables, op
     elseif option == "thetal_maxgrad"
 
         @inbounds for k in real_center_indicies(self.Gr)
-            grad = (GMV.THL.values[k + 1] - GMV.THL.values[k]) * self.Gr.dzi
-            if grad > maxgrad
-                maxgrad = grad
+            ∇θ_liq = ∇_upwind(GMV.THL.values, self.Gr, k)
+            if ∇θ_liq > ∇θ_liq_max
+                ∇θ_liq_max = ∇θ_liq
                 self.zi = self.Gr.z[k]
             end
         end

src/TurbulenceConvection.jl

@@ -2,6 +2,7 @@ module TurbulenceConvection
 
 import Dierckx
 using StatsBase
+using StaticArrays
 import Statistics
 using PoissonRandom: pois_rand
 import LambertW
@@ -33,6 +34,7 @@ include("python_primitives.jl")
 include("parameters.jl")
 include("Grid.jl")
 include("Fields.jl")
+include("Operators.jl")
 include("NetCDFIO.jl")
 include("ReferenceState.jl")
 include("types.jl")

src/Turbulence_PrognosticTKE.jl

@@ -411,9 +411,9 @@ function compute_mixing_length(self, obukhov_length, ustar, GMV::GridMeanVariabl
 
         # Buoyancy-shear-subdomain exchange-dissipation TKE equilibrium scale
         shear2 =
-            pow((GMV.U.values[k + 1] - GMV.U.values[k - 1]) * 0.5 * dzi, 2) +
-            pow((GMV.V.values[k + 1] - GMV.V.values[k - 1]) * 0.5 * dzi, 2) +
-            pow((self.EnvVar.W.values[k] - self.EnvVar.W.values[k - 1]) * dzi, 2)
+            pow(∇_collocated(GMV.U.values, grid, k), 2) +
+            pow(∇_collocated(GMV.V.values, grid, k), 2) +
+            pow(∇f2c(self.EnvVar.W.values, grid, k), 2)
         qt_dry = self.EnvThermo.qt_dry[k]
         th_dry = self.EnvThermo.th_dry[k]
         t_cloudy = self.EnvThermo.t_cloudy[k]
@@ -428,6 +428,8 @@ function compute_mixing_length(self, obukhov_length, ustar, GMV::GridMeanVariabl
         grad_qt_plus = (self.EnvVar.QT.values[k + 1] - self.EnvVar.QT.values[k]) * dzi
         grad_thl = interp2pt(grad_thl_low, grad_thl_plus)
         grad_qt = interp2pt(grad_qt_low, grad_qt_plus)
+        # grad_thl = ∇_collocated(self.EnvVar.THL.values, grid, k)
+        # grad_qt = ∇_collocated(self.EnvVar.QT.values, grid, k)
         # g/theta_ref
         prefactor = g * (Rd / ref_state.alpha0_half[k] / ref_state.p0_half[k]) * exner_c(ref_state.p0_half[k])
         d_buoy_thetal_dry = prefactor * (1.0 + (eps_vi - 1.0) * qt_dry)