Merge #406

406: Make updraft w, theta, and q_tot fields r=charleskawczynski a=charleskawczynski

This PR also changes the names to abide by our naming conventions, and use unicode.

Co-authored-by: Charles Kawczynski <[email protected]>

integration_tests/utils/main.jl

@@ -67,15 +67,15 @@ face_diagnostic_vars(FT, n_up) = (; face_diagnostic_vars_gm(FT)..., face_diagnos
 # Center only
 cent_prognostic_vars(FT, n_up) = (; cent_prognostic_vars_gm(FT)..., cent_prognostic_vars_edmf(FT, n_up)...)
 cent_prognostic_vars_gm(FT) = (; U = FT(0), V = FT(0), H = FT(0), QT = FT(0))
-cent_prognostic_vars_up(FT) = (; area = FT(0), H = FT(0), QT = FT(0))
+cent_prognostic_vars_up(FT) = (; area = FT(0), θ_liq_ice = FT(0), q_tot = FT(0))
 cent_prognostic_vars_en(FT) = (; TKE = FT(0), Hvar = FT(0), QTvar = FT(0), HQTcov = FT(0))
 cent_prognostic_vars_edmf(FT, n_up) =
     (; turbconv = (; en = cent_prognostic_vars_en(FT), up = ntuple(i -> cent_prognostic_vars_up(FT), n_up)))
 # cent_prognostic_vars_edmf(FT, n_up) = (;) # could also use this for empty model
 
 # Face only
 face_prognostic_vars(FT, n_up) = (; face_prognostic_vars_edmf(FT, n_up)...)
-face_prognostic_vars_up(FT) = (; W = FT(0))
+face_prognostic_vars_up(FT) = (; w = FT(0))
 face_prognostic_vars_edmf(FT, n_up) = (; turbconv = (; up = ntuple(i -> face_prognostic_vars_up(FT), n_up)))
 # face_prognostic_vars_edmf(FT, n_up) = (;) # could also use this for empty model
 

src/EDMF_Updrafts.jl

@@ -2,9 +2,13 @@ function initialize(edmf, grid, state, up::UpdraftVariables, gm::GridMeanVariabl
     kc_surf = kc_surface(grid)
 
     prog_up = center_prog_updrafts(state)
-    up.W.values .= 0
+    prog_up_f = face_prog_updrafts(state)
     up.B.values .= 0
     @inbounds for i in 1:(up.n_updrafts)
+        @inbounds for k in real_face_indices(grid)
+            prog_up_f[i].w[k] = 0
+        end
+
         @inbounds for k in real_center_indices(grid)
             # Simple treatment for now, revise when multiple updraft closures
             # become more well defined
@@ -13,9 +17,9 @@ function initialize(edmf, grid, state, up::UpdraftVariables, gm::GridMeanVariabl
             else
                 prog_up[i].area[k] = up.updraft_fraction / up.n_updrafts
             end
-            up.QT.values[i, k] = gm.QT.values[k]
+            prog_up[i].q_tot[k] = gm.QT.values[k]
+            prog_up[i].θ_liq_ice[k] = gm.H.values[k]
             up.QL.values[i, k] = gm.QL.values[k]
-            up.H.values[i, k] = gm.H.values[k]
             up.T.values[i, k] = gm.T.values[k]
         end
 
@@ -95,28 +99,29 @@ function initialize_DryBubble(edmf, grid, state, up::UpdraftVariables, gm::GridM
     #! format: on
 
     prog_up = center_prog_updrafts(state)
+    prog_up_f = face_prog_updrafts(state)
     Area_in = pyinterp(grid.zc, z_in, Area_in)
     θ_liq_in = pyinterp(grid.zc, z_in, θ_liq_in)
     T_in = pyinterp(grid.zc, z_in, T_in)
     @inbounds for i in 1:(up.n_updrafts)
         @inbounds for k in real_face_indices(grid)
             if minimum(z_in) <= grid.zf[k] <= maximum(z_in)
-                up.W.values[i, k] = 0.0
+                prog_up_f[i].w[k] = 0.0
             end
         end
 
         @inbounds for k in real_center_indices(grid)
             if minimum(z_in) <= grid.zc[k] <= maximum(z_in)
                 prog_up[i].area[k] = Area_in[k] #up.updraft_fraction/up.n_updrafts
-                up.H.values[i, k] = θ_liq_in[k]
-                up.QT.values[i, k] = 0.0
+                prog_up[i].θ_liq_ice[k] = θ_liq_in[k]
+                prog_up[i].q_tot[k] = 0.0
                 up.QL.values[i, k] = 0.0
 
                 # for now temperature is provided as diagnostics from LES
                 up.T.values[i, k] = T_in[k]
             else
                 prog_up[i].area[k] = 0.0 #up.updraft_fraction/up.n_updrafts
-                up.H.values[i, k] = gm.H.values[k]
+                prog_up[i].θ_liq_ice[k] = gm.H.values[k]
                 up.T.values[i, k] = gm.T.values[k]
             end
         end
@@ -144,15 +149,16 @@ end
 # quick utility to set "new" arrays with values in the "values" arrays
 function set_new_with_values(up::UpdraftVariables, grid, state)
     prog_up = center_prog_updrafts(state)
+    prog_up_f = face_prog_updrafts(state)
     @inbounds for i in 1:(up.n_updrafts)
         @inbounds for k in real_face_indices(grid)
-            up.W.new[i, k] = up.W.values[i, k]
+            up.W.new[i, k] = prog_up_f[i].w[k]
         end
 
         @inbounds for k in real_center_indices(grid)
             up.Area.new[i, k] = prog_up[i].area[k]
-            up.QT.new[i, k] = up.QT.values[i, k]
-            up.H.new[i, k] = up.H.values[i, k]
+            up.QT.new[i, k] = prog_up[i].q_tot[k]
+            up.H.new[i, k] = prog_up[i].θ_liq_ice[k]
         end
     end
     return
@@ -161,16 +167,16 @@ end
 # quick utility to set "tmp" arrays with values in the "new" arrays
 function set_values_with_new(up::UpdraftVariables, grid, state)
     prog_up = center_prog_updrafts(state)
+    prog_up_f = face_prog_updrafts(state)
     @inbounds for i in 1:(up.n_updrafts)
         @inbounds for k in real_face_indices(grid)
-            up.W.values[i, k] = up.W.new[i, k]
+            prog_up_f[i].w[k] = up.W.new[i, k]
         end
 
         @inbounds for k in real_center_indices(grid)
-            up.W.values[i, k] = up.W.new[i, k]
             prog_up[i].area[k] = up.Area.new[i, k]
-            up.QT.values[i, k] = up.QT.new[i, k]
-            up.H.values[i, k] = up.H.new[i, k]
+            prog_up[i].q_tot[k] = up.QT.new[i, k]
+            prog_up[i].θ_liq_ice[k] = up.H.new[i, k]
         end
     end
     return
@@ -244,7 +250,7 @@ function compute_rain_formation_tendencies(
     @inbounds for i in 1:(up.n_updrafts)
         @inbounds for k in real_center_indices(grid)
             T_up = up.T.values[i, k]
-            q_tot_up = up.QT.values[i, k]
+            q_tot_up = prog_up[i].q_tot[k]
             ts_up = TD.PhaseEquil_pTq(param_set, p0_c[k], T_up, q_tot_up)
 
             # autoconversion and accretion

src/Operators.jl

@@ -325,7 +325,7 @@ Used when
      - traversing cell centers
      - grabbing centered center stencils
 """
-function ccut(f, grid, k)
+function ccut(f, grid, k::Cent)
     if is_surface_center(grid, k)
         return SA.SVector(f[k], f[k + 1])
     elseif is_toa_center(grid, k)
@@ -334,7 +334,7 @@ function ccut(f, grid, k)
         return SA.SVector(f[k - 1], f[k], f[k + 1])
     end
 end
-function ccut(f, grid, k, i_up::Int)
+function ccut(f, grid, k::Cent, i_up::Int)
     if is_surface_center(grid, k)
         return SA.SVector(f[i_up, k], f[i_up, k + 1])
     elseif is_toa_center(grid, k)
@@ -351,7 +351,7 @@ Used when
      - traversing cell faces
      - grabbing centered face stencils
 """
-function fcut(f, grid, k)
+function fcut(f, grid, k::CCO.PlusHalf)
     if is_surface_face(grid, k)
         return SA.SVector(f[k], f[k + 1])
     elseif is_toa_face(grid, k)
@@ -360,7 +360,7 @@ function fcut(f, grid, k)
         return SA.SVector(f[k - 1], f[k], f[k + 1])
     end
 end
-function fcut(f, grid, k, i_up::Int)
+function fcut(f, grid, k::CCO.PlusHalf, i_up::Int)
     if is_surface_face(grid, k)
         return SA.SVector(f[i_up, k], f[i_up, k + 1])
     elseif is_toa_face(grid, k)
@@ -379,14 +379,14 @@ Used when
 
 This is needed for "upwinding" rain, which travels _down_ (hence the direction change).
 """
-function ccut_downwind(f, grid, k)
+function ccut_downwind(f, grid, k::Cent)
     if is_toa_center(grid, k)
         return SA.SVector(f[k])
     else
         return SA.SVector(f[k], f[k + 1])
     end
 end
-function ccut_downwind(f, grid, k, i_up::Int)
+function ccut_downwind(f, grid, k::Cent, i_up::Int)
     if is_toa_center(grid, k)
         return SA.SVector(f[i_up, k])
     else
@@ -401,14 +401,14 @@ Used when
      - traversing cell centers
      - grabbing one-sided (upwind) stencil of cell center `k` and cell center `k-1`
 """
-function ccut_upwind(f, grid, k)
+function ccut_upwind(f, grid, k::Cent)
     if is_surface_center(grid, k)
         return SA.SVector(f[k])
     else
         return SA.SVector(f[k - 1], f[k])
     end
 end
-function ccut_upwind(f, grid, k, i_up::Int)
+function ccut_upwind(f, grid, k::Cent, i_up::Int)
     if is_surface_center(grid, k)
         return SA.SVector(f[i_up, k])
     else
@@ -477,18 +477,20 @@ Used when
      - traversing cell centers
      - grabbing _interpolated_ one-sided (upwind) stencil of cell center `k` and cell center `k-1`
 """
-function daul_f2c_upwind(f, grid, k)
+function daul_f2c_upwind(f, grid, k::Cent)
+    kf = CCO.PlusHalf(k.i)
     if is_surface_center(grid, k)
-        return SA.SVector((f[k] + f[k + 1]) / 2)
+        return SA.SVector((f[kf] + f[kf + 1]) / 2)
     else
-        return SA.SVector((f[k - 1] + f[k]) / 2, (f[k] + f[k + 1]) / 2)
+        return SA.SVector((f[kf - 1] + f[kf]) / 2, (f[kf] + f[kf + 1]) / 2)
     end
 end
-function daul_f2c_upwind(f, grid, k, i_up::Int)
+function daul_f2c_upwind(f, grid, k::Cent, i_up::Int)
+    kf = CCO.PlusHalf(k.i)
     if is_surface_center(grid, k)
-        return SA.SVector((f[i_up, k] + f[i_up, k + 1]) / 2)
+        return SA.SVector((f[i_up, kf] + f[i_up, kf + 1]) / 2)
     else
-        return SA.SVector((f[i_up, k - 1] + f[i_up, k]) / 2, (f[i_up, k] + f[i_up, k + 1]) / 2)
+        return SA.SVector((f[i_up, kf - 1] + f[i_up, kf]) / 2, (f[i_up, kf] + f[i_up, kf + 1]) / 2)
     end
 end
 
@@ -499,8 +501,14 @@ Used when
      - traversing cell centers
      - grabbing stencil of 2 neighboring cell faces
 """
-dual_faces(f, grid, k) = SA.SVector(f[k], f[k + 1])
-dual_faces(f, grid, k, i_up::Int) = SA.SVector(f[i_up, k], f[i_up, k + 1])
+function dual_faces(f, grid, k::Cent)
+    kf = CCO.PlusHalf(k.i)
+    SA.SVector(f[kf], f[kf + 1])
+end
+function dual_faces(f, grid, k::Cent, i_up::Int)
+    kf = CCO.PlusHalf(k.i)
+    SA.SVector(f[i_up, kf], f[i_up, kf + 1])
+end
 
 """
     dual_centers
@@ -613,8 +621,6 @@ function construct_tridiag_diffusion_en(
     TS,
     KM,
     KH,
-    prog_up,
-    w_up,
     w_en,
     tke_en,
     n_updrafts::Int,
@@ -637,6 +643,8 @@ function construct_tridiag_diffusion_en(
     ρ0_c = center_ref_state(state).ρ0
     ρ0_f = face_ref_state(state).ρ0
     aux_tc = center_aux_tc(state)
+    prog_up = center_prog_updrafts(state)
+    prog_up_f = face_prog_updrafts(state)
 
     ae = vec(1 .- aux_tc.bulk.area)
     rho_ae_K_m = face_field(grid)
@@ -661,7 +669,7 @@ function construct_tridiag_diffusion_en(
             if prog_up[i].area[k] > minimum_area
                 turb_entr = frac_turb_entr[i, k]
                 R_up = pressure_plume_spacing[i]
-                w_up_c = interpf2c(w_up, grid, k, i)
+                w_up_c = interpf2c(prog_up_f[i].w, grid, k)
                 D_env += ρ0_c[k] * prog_up[i].area[k] * w_up_c * (entr_sc[i, k] + turb_entr)
             else
                 D_env = 0.0