switch to conservative variables in the updrafts

integration_tests/utils/initial_conditions.jl

@@ -33,33 +33,32 @@ end
 
 function initialize_updrafts(edmf, grid, state, up::TC.UpdraftVariables, gm::TC.GridMeanVariables)
     kc_surf = TC.kc_surface(grid)
-
-    prog_up = TC.center_prog_updrafts(state)
+    aux_up = TC.center_aux_updrafts(state)
     prog_gm = TC.center_prog_grid_mean(state)
     aux_up = TC.center_aux_updrafts(state)
-    prog_up_f = TC.face_prog_updrafts(state)
+    aux_up_f = TC.face_aux_updrafts(state)
     aux_gm = TC.center_aux_grid_mean(state)
     @inbounds for i in 1:(up.n_updrafts)
         @inbounds for k in TC.real_face_indices(grid)
-            prog_up_f[i].w[k] = 0
+            aux_up_f[i].w[k] = 0
         end
 
         @inbounds for k in TC.real_center_indices(grid)
             aux_up[i].buoy[k] = 0
             # Simple treatment for now, revise when multiple updraft closures
             # become more well defined
             if up.prognostic
-                prog_up[i].area[k] = 0.0 #up.updraft_fraction/up.n_updrafts
+                aux_up[i].area[k] = 0.0 #up.updraft_fraction/up.n_updrafts
             else
-                prog_up[i].area[k] = up.updraft_fraction / up.n_updrafts
+                aux_up[i].area[k] = up.updraft_fraction / up.n_updrafts
             end
-            prog_up[i].q_tot[k] = prog_gm.q_tot[k]
-            prog_up[i].θ_liq_ice[k] = prog_gm.θ_liq_ice[k]
+            aux_up[i].q_tot[k] = prog_gm.q_tot[k]
+            aux_up[i].θ_liq_ice[k] = prog_gm.θ_liq_ice[k]
             aux_up[i].q_liq[k] = aux_gm.q_liq[k]
             aux_up[i].T[k] = aux_gm.T[k]
         end
 
-        prog_up[i].area[kc_surf] = up.updraft_fraction / up.n_updrafts
+        aux_up[i].area[kc_surf] = up.updraft_fraction / up.n_updrafts
     end
     return
 end
@@ -134,33 +133,32 @@ function initialize_updrafts_DryBubble(edmf, grid, state, up::TC.UpdraftVariable
         264.1574, 263.6518, 263.1461, 262.6451, 262.1476, 261.6524]
     #! format: on
 
-    prog_up = TC.center_prog_updrafts(state)
-    aux_up = TC.center_aux_updrafts(state)
-    prog_up_f = TC.face_prog_updrafts(state)
-    aux_gm = TC.center_aux_grid_mean(state)
-    prog_gm = TC.center_prog_grid_mean(state)
-    Area_in = TC.pyinterp(grid.zc, z_in, Area_in)
-    θ_liq_in = TC.pyinterp(grid.zc, z_in, θ_liq_in)
-    T_in = TC.pyinterp(grid.zc, z_in, T_in)
+    aux_up = center_aux_updrafts(state)
+    aux_up_f = face_aux_updrafts(state)
+    aux_gm = center_aux_grid_mean(state)
+    prog_gm = center_prog_grid_mean(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 TC.real_face_indices(grid)
+        @inbounds for k in real_face_indices(grid)
             if minimum(z_in) <= grid.zf[k] <= maximum(z_in)
-                prog_up_f[i].w[k] = 0.0
+                aux_up_f[i].w[k] = 0.0
             end
         end
 
-        @inbounds for k in TC.real_center_indices(grid)
+        @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
-                prog_up[i].θ_liq_ice[k] = θ_liq_in[k]
-                prog_up[i].q_tot[k] = 0.0
+                aux_up[i].area[k] = Area_in[k] #up.updraft_fraction/up.n_updrafts
+                aux_up[i].θ_liq_ice[k] = θ_liq_in[k]
+                aux_up[i].q_tot[k] = 0.0
                 aux_up[i].q_liq[k] = 0.0
 
                 # for now temperature is provided as diagnostics from LES
                 aux_up[i].T[k] = T_in[k]
             else
-                prog_up[i].area[k] = 0.0 #up.updraft_fraction/up.n_updrafts
-                prog_up[i].θ_liq_ice[k] = prog_gm.θ_liq_ice[k]
+                aux_up[i].area[k] = 0.0 #up.updraft_fraction/up.n_updrafts
+                aux_up[i].θ_liq_ice[k] = prog_gm.θ_liq_ice[k]
                 aux_up[i].T[k] = aux_gm.T[k]
             end
         end

integration_tests/utils/main.jl

@@ -62,7 +62,7 @@ cent_aux_vars_en_2m(FT) = (;
     rain_src = FT(0),
 )
 cent_aux_vars_up(FT) =
-    (; q_liq = FT(0), T = FT(0), RH = FT(0), buoy = FT(0), area_new = FT(0), q_tot_new = FT(0), θ_liq_ice_new = FT(0))
+    (; q_liq = FT(0), T = FT(0), RH = FT(0), buoy = FT(0), area = FT(0), q_tot = FT(0), θ_liq_ice = FT(0))
 cent_aux_vars_edmf(FT, n_up) = (;
     turbconv = (;
         bulk = (; area = FT(0), θ_liq_ice = FT(0), RH = FT(0), buoy = FT(0), q_tot = FT(0), q_liq = FT(0), T = FT(0)),
@@ -91,7 +91,7 @@ cent_aux_vars(FT, n_up) = (; aux_vars_ref_state(FT)..., cent_aux_vars_gm(FT)...,
 
 # Face only
 face_aux_vars_gm(FT) = ()
-face_aux_vars_up(FT) = (; w_new = FT(0))
+face_aux_vars_up(FT) = (; w = FT(0))
 
 face_aux_vars_edmf(FT, n_up) = (;
     turbconv = (;
@@ -125,7 +125,7 @@ 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), θ_liq_ice = FT(0), q_tot = FT(0))
-cent_prognostic_vars_up(FT) = (; area = FT(0), θ_liq_ice = FT(0), q_tot = FT(0))
+cent_prognostic_vars_up(FT) = (; ρarea = FT(0), ρaθ_liq_ice = FT(0), ρaq_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 = (;
@@ -136,7 +136,7 @@ cent_prognostic_vars_edmf(FT, n_up) = (;
 
 # Face only
 face_prognostic_vars(FT, n_up) = (; w = FT(0), face_prognostic_vars_edmf(FT, n_up)...)
-face_prognostic_vars_up(FT) = (; w = FT(0))
+face_prognostic_vars_up(FT) = (; ρaw = 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

@@ -1,4 +1,3 @@
-
 function initialize_io(up::UpdraftVariables, Stats::NetCDFIO_Stats)
     add_profile(Stats, "updraft_cloud_fraction")
 
@@ -26,7 +25,7 @@ end
 function upd_cloud_diagnostics(up::UpdraftVariables, grid, state)
     up.lwp = 0.0
 
-    prog_up = center_prog_updrafts(state)
+    aux_up = center_aux_updrafts(state)
     aux_up = center_aux_updrafts(state)
     ρ0_c = center_ref_state(state).ρ0
     @inbounds for i in 1:(up.n_updrafts)
@@ -36,14 +35,14 @@ function upd_cloud_diagnostics(up::UpdraftVariables, grid, state)
         up.cloud_cover[i] = 0.0
 
         @inbounds for k in real_center_indices(grid)
-            if prog_up[i].area[k] > 1e-3
+            if aux_up[i].area[k] > 1e-3
                 up.updraft_top[i] = max(up.updraft_top[i], grid.zc[k])
-                up.lwp += ρ0_c[k] * aux_up[i].q_liq[k] * prog_up[i].area[k] * grid.Δz
+                up.lwp += ρ0_c[k] * aux_up[i].q_liq[k] * aux_up[i].area[k] * grid.Δz
 
                 if aux_up[i].q_liq[k] > 1e-8
                     up.cloud_base[i] = min(up.cloud_base[i], grid.zc[k])
                     up.cloud_top[i] = max(up.cloud_top[i], grid.zc[k])
-                    up.cloud_cover[i] = max(up.cloud_cover[i], prog_up[i].area[k])
+                    up.cloud_cover[i] = max(up.cloud_cover[i], aux_up[i].area[k])
                 end
             end
         end
@@ -65,29 +64,28 @@ function compute_rain_formation_tendencies(
 )
     p0_c = center_ref_state(state).p0
     ρ0_c = center_ref_state(state).ρ0
-    prog_up = center_prog_updrafts(state)
     aux_up = center_aux_updrafts(state)
     prog_ra = center_prog_rain(state)
     tendencies_ra = center_tendencies_rain(state)
 
     @inbounds for i in 1:(up.n_updrafts)
         @inbounds for k in real_center_indices(grid)
             T_up = aux_up[i].T[k]
-            q_tot_up = prog_up[i].q_tot[k]
+            q_tot_up = aux_up[i].q_tot[k]
             ts_up = TD.PhaseEquil_pTq(param_set, p0_c[k], T_up, q_tot_up)
 
             # autoconversion and accretion
             mph = precipitation_formation(
                 param_set,
                 rain.rain_model,
                 prog_ra.qr[k],
-                prog_up[i].area[k],
+                aux_up[i].area[k],
                 ρ0_c[k],
                 dt,
                 ts_up,
             )
-            up_thermo.qt_tendency_rain_formation[i, k] = mph.qt_tendency * prog_up[i].area[k]
-            up_thermo.θ_liq_ice_tendency_rain_formation[i, k] = mph.θ_liq_ice_tendency * prog_up[i].area[k]
+            up_thermo.qt_tendency_rain_formation[i, k] = mph.qt_tendency * aux_up[i].area[k]
+            up_thermo.θ_liq_ice_tendency_rain_formation[i, k] = mph.θ_liq_ice_tendency * aux_up[i].area[k]
         end
     end
     # TODO - to be deleted once we sum all tendencies elsewhere

src/Operators.jl

@@ -576,9 +576,9 @@ function construct_tridiag_diffusion_en(
     ρ0_f = face_ref_state(state).ρ0
     aux_tc = center_aux_tc(state)
     w_en = face_aux_environment(state).w
-    prog_up = center_prog_updrafts(state)
-    prog_up_f = face_prog_updrafts(state)
+    aux_up_f = face_aux_updrafts(state)
     prog_en = center_prog_environment(state)
+    aux_up = center_aux_updrafts(state)
 
     ae = 1 .- aux_tc.bulk.area
     rho_ae_K_m = face_field(grid)
@@ -602,11 +602,11 @@ function construct_tridiag_diffusion_en(
         D_env = 0.0
 
         @inbounds for i in 1:n_updrafts
-            if prog_up[i].area[k] > minimum_area
+            if aux_up[i].area[k] > minimum_area
                 turb_entr = frac_turb_entr[i, k]
                 R_up = pressure_plume_spacing[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)
+                w_up_c = interpf2c(aux_up_f[i].w, grid, k)
+                D_env += ρ0_c[k] * aux_up[i].area[k] * w_up_c * (entr_sc[i, k] + turb_entr)
             else
                 D_env = 0.0
             end