Merge #1826

1826: fix turb_source bug and set CBL_model similar to SBL_model r=charleskawczynski a=yairchn

### Description
This PR take the changes to convective_bl_model.jl from #1817 that do not effect the output of a test and adds them to master separately. 

- [x] Code follows the [style guidelines](https://clima.github.io/ClimateMachine.jl/latest/DevDocs/CodeStyle/) OR N/A.
- [x] Unit tests are included OR N/A.
- [x] Code is exercised in an integration test OR N/A.
- [x] Documentation has been added/updated OR N/A.


Co-authored-by: yairchn <[email protected]>

experiments/AtmosLES/convective_bl_les.jl

@@ -63,28 +63,16 @@ function main()
         driver_config,
         init_on_cpu = true,
         Courant_number = CFLmax,
-        CFL_direction = HorizontalDirection(),
     )
     dgn_config = config_diagnostics(driver_config)
 
-    cbtmarfilter = GenericCallbacks.EveryXSimulationSteps(1) do
-        Filters.apply!(
-            solver_config.Q,
-            ("moisture.ρq_tot",),
-            solver_config.dg.grid,
-            TMARFilter(),
-        )
-        nothing
-    end
-
     check_cons = (
         ClimateMachine.ConservationCheck("ρ", "1mins", FT(0.0001)),
         ClimateMachine.ConservationCheck("ρe", "1mins", FT(0.0025)),
     )
 
     result = ClimateMachine.invoke!(
         solver_config;
-        user_callbacks = (cbtmarfilter,),
         diagnostics_config = dgn_config,
         check_cons = check_cons,
         check_euclidean_distance = true,

experiments/AtmosLES/convective_bl_model.jl

@@ -126,7 +126,7 @@ atmos_source!(::ConvectiveBLSponge, args...) = nothing
 """
   Initial Condition for ConvectiveBoundaryLayer LES
 """
-function init_convective_bl!(problem, bl, state, aux, localgeo, t)
+function init_problem!(problem, bl, state, aux, localgeo, t)
     (x, y, z) = localgeo.coord
 
     # Problem floating point precision
@@ -167,7 +167,9 @@ function init_convective_bl!(problem, bl, state, aux, localgeo, t)
     state.ρ = ρ
     state.ρu = SVector(ρu, ρv, ρw)
     state.ρe = ρe_tot
-    state.moisture.ρq_tot = ρ * q_tot
+    if !(bl.moisture isa DryModel)
+        state.moisture.ρq_tot = ρ * q_tot
+    end
 
     if z <= FT(400) # Add random perturbations to bottom 400m of model
         state.ρe += rand() * ρe_tot / 100
@@ -194,24 +196,28 @@ function convective_bl_model(
     zmax,
     surface_flux;
     turbconv = NoTurbConv(),
+    moisture_model = "dry",
 ) where {FT}
 
-    ics = init_convective_bl!     # Initial conditions
+    ics = init_problem!     # Initial conditions
 
     C_smag = FT(0.23)     # Smagorinsky coefficient
     C_drag = FT(0.001)    # Momentum exchange coefficient
     z_sponge = FT(2560)     # Start of sponge layer
+
     α_max = FT(0.75)       # Strength of sponge layer (timescale)
+
     γ = 2                  # Strength of sponge layer (exponent)
     u_geostrophic = FT(4)        # Eastward relaxation speed
     u_slope = FT(0)              # Slope of altitude-dependent relaxation speed
     v_geostrophic = FT(0)        # Northward relaxation speed
     f_coriolis = FT(1.031e-4) # Coriolis parameter
     u_star = FT(0.3)
     q_sfc = FT(0)
+    moisture_flux = FT(0)
 
     # Assemble source components
-    source = (
+    source_default = (
         Gravity(),
         ConvectiveBLSponge(
             FT,
@@ -230,8 +236,26 @@ function convective_bl_model(
             u_slope,
             v_geostrophic,
         ),
+        turbconv_sources(turbconv)...,
     )
 
+    if moisture_model == "dry"
+        source = source_default
+        moisture = DryModel()
+    elseif moisture_model == "equilibrium"
+        source = source_default
+        moisture = EquilMoist{FT}(; maxiter = 5, tolerance = FT(0.1))
+    elseif moisture_model == "nonequilibrium"
+        source = (source_default..., CreateClouds()...)
+        moisture = NonEquilMoist()
+    else
+        @warn @sprintf(
+            """
+%s: unrecognized moisture_model in source terms, using the defaults""",
+            moisture_model,
+        )
+        source = source_default
+    end
     # Set up problem initial and boundary conditions
     if surface_flux == "prescribed"
         energy_bc = PrescribedEnergyFlux((state, aux, t) -> LHF + SHF)
@@ -254,10 +278,21 @@ function convective_bl_model(
         )
     end
 
-    # Set up problem initial and boundary conditions
-    moisture_flux = FT(0)
-    problem = AtmosProblem(
-        boundaryconditions = (
+    if moisture_model == "dry"
+        boundary_conditions = (
+            AtmosBC(
+                momentum = Impenetrable(DragLaw(
+                    # normPu_int is the internal horizontal speed
+                    # P represents the projection onto the horizontal
+                    (state, aux, t, normPu_int) -> (u_star / normPu_int)^2,
+                )),
+                energy = energy_bc,
+                turbconv = turbconv_bcs(turbconv)[1],
+            ),
+            AtmosBC(turbconv = turbconv_bcs(turbconv)[2]),
+        )
+    else
+        boundary_conditions = (
             AtmosBC(
                 momentum = Impenetrable(DragLaw(
                     # normPu_int is the internal horizontal speed
@@ -269,20 +304,25 @@ function convective_bl_model(
                 turbconv = turbconv_bcs(turbconv)[1],
             ),
             AtmosBC(turbconv = turbconv_bcs(turbconv)[2]),
-        ),
+        )
+    end
+
+    problem = AtmosProblem(
         init_state_prognostic = ics,
+        boundaryconditions = boundary_conditions,
     )
 
     # Assemble model components
     model = AtmosModel{FT}(
-        AtmosLESConfigType,
+        config_type,
         param_set;
         problem = problem,
         turbulence = SmagorinskyLilly{FT}(C_smag),
-        moisture = EquilMoist{FT}(; maxiter = 5, tolerance = FT(0.1)),
+        moisture = moisture,
         source = source,
         turbconv = turbconv,
     )
+
     return model
 end