Merge pull request #2937 from CliMA/as/deepatmos-init

Deep atmosphere init for baroclinic wave problem

Author: akshaysridhar

config/default_configs/default_config.yml

@@ -278,7 +278,7 @@ log_progress:
   help: "Log simulation progress, including wall-time estimates"
   value: true
 deep_atmosphere:
-  help: "Do not assume that the atmosphere is shallow, where the vertical columns are cylindrical"
+  help: "If true, use deep atmosphere equations and metric terms, otherwise assume columns are cylindrical (shallow atmosphere) [`true`, `false` (default)]"
   value: false
 restart_file:
   help: "Path to HDF5 file to use as simulation starting point"

config/model_configs/sphere_baroclinic_wave_rhoe.yml

@@ -1,11 +1,12 @@
 dt_save_state_to_disk: "2days"
 initial_condition: "DryBaroclinicWave"
+deep_atmosphere: true
 dt: "400secs"
 t_end: "10days"
 job_id: "sphere_baroclinic_wave_rhoe"
 diagnostics:
-  - short_name: [pfull, wa, va, rv, ke]
+  - short_name: [pfull, ua, wa, va, rv, ta, ke]
     period: 1days
-  - short_name: [pfull, wa, va, rv]
+  - short_name: [pfull, ua, wa, va, rv, ta, ke]
     period: 1days
     writer: h5

post_processing/ci_plots.jl

@@ -610,7 +610,10 @@ function make_plots(
     )
 end
 
-DryBaroWavePlots = Union{Val{:sphere_baroclinic_wave_rhoe}}
+DryBaroWavePlots = Union{
+    Val{:sphere_baroclinic_wave_rhoe},
+    Val{:sphere_baroclinic_wave_rhoe_deepatmos},
+}
 
 function make_plots(::DryBaroWavePlots, output_paths::Vector{<:AbstractString})
     simdirs = SimDir.(output_paths)
@@ -665,6 +668,7 @@ end
 
 function make_plots(
     ::Val{:sphere_baroclinic_wave_rhoe_equilmoist},
+    ::Val{:sphere_baroclinic_wave_rhoe_equilmoist_deepatmos},
     output_paths::Vector{<:AbstractString},
 )
     simdirs = SimDir.(output_paths)

regression_tests/ref_counter.jl

@@ -1,5 +1,10 @@
-159
+160
 
+# 160: 
+# - Introduces initial conditions for the baroclinic-wave
+#   test case in a deep-atmosphere configuration. Modifies
+#   existing config to use `deep_atmosphere` mode. 
+#
 # 159:
 # - Changed the boundary condition of edmf updraft properties
 #   to be dependent on the surface area

src/initial_conditions/initial_conditions.jl

@@ -352,7 +352,7 @@ end
 ## Baroclinic Wave
 ##
 
-function baroclinic_wave_values(z, ϕ, λ, params, perturb)
+function shallow_atmos_baroclinic_wave_values(z, ϕ, λ, params, perturb)
     FT = eltype(params)
     R_d = CAP.R_d(params)
     MSLP = CAP.MSLP(params)
@@ -418,7 +418,80 @@ function baroclinic_wave_values(z, ϕ, λ, params, perturb)
     return (; T_v, p, u, v)
 end
 
-function moist_baroclinic_wave_values(z, ϕ, λ, params, perturb)
+function deep_atmos_baroclinic_wave_values(z, ϕ, λ, params, perturb)
+    FT = eltype(params)
+    R_d = CAP.R_d(params)
+    MSLP = CAP.MSLP(params)
+    grav = CAP.grav(params)
+    Ω = CAP.Omega(params)
+    R = CAP.planet_radius(params)
+
+    # Constants from paper (See Table 1. in Ullrich et al (2014)) 
+    k = 3         # Power for temperature field
+    T_e = FT(310) # Surface temperature at the equator
+    T_p = FT(240) # Surface temperature at the pole
+    T_0 = FT(0.5) * (T_e + T_p)
+    Γ = FT(0.005) # Lapse rate
+    A = 1 / Γ  # (Eq 16)
+    B = (T_0 - T_p) / T_0 / T_p # (Eq 17)
+    C = FT(0.5) * (k + 2) * (T_e - T_p) / T_e / T_p # (Eq 17)
+    b = 2 # half-width parameter
+    H = R_d * T_0 / grav
+    z_t = FT(15e3) # Top of perturbation domain
+    λ_c = FT(20) # Geographical location (λ dim) of perturbation center
+    ϕ_c = FT(40) # Geographical location (ϕ dim) of perturbation center
+    d_0 = R / 6
+    V_p = FT(1)
+
+    # Virtual temperature and pressure
+    τ̃₁ =
+        A * Γ / T_0 * exp(Γ * z / T_0) +
+        B * (1 - 2 * (z / b / H)^2) * exp(-(z / b / H)^2)# (Eq 14)
+    τ̃₂ = C * (1 - 2 * (z / b / H)^2) * exp(-(z / b / H)^2) # (Eq 15)
+    ∫τ̃₁ = (A * (exp(Γ * z / T_0) - 1)) + B * z * exp(-(z / b / H)^2) # (Eq A1)
+    ∫τ̃₂ = C * z * exp(-(z / b / H)^2) # (Eq A2)
+    I_T =
+        ((z + R) / R * cosd(ϕ))^k -
+        (k / (k + 2)) * ((z + R) / R * cosd(ϕ))^(k + 2)
+    T_v = FT((R / (z + R))^2 * (τ̃₁ - τ̃₂ * I_T)^(-1)) # (Eq A3)
+    p = FT(MSLP * exp(-grav / R_d * (∫τ̃₁ - ∫τ̃₂ * I_T))) # (Eq A6)
+    # Horizontal velocity
+    U =
+        grav / R *
+        k *
+        T_v *
+        ∫τ̃₂ *
+        (((z + R) * cosd(ϕ) / R)^(k - 1) - ((R + z) * cosd(ϕ) / R)^(k + 1)) # wind-proxy (Eq A4)
+    u = FT(
+        -Ω * (R + z) * cosd(ϕ) +
+        sqrt((Ω * (R + z) * cosd(ϕ))^2 + (R + z) * cosd(ϕ) * U),
+    )
+    v = FT(0)
+    if perturb
+        F_z = (1 - 3 * (z / z_t)^2 + 2 * (z / z_t)^3) * (z ≤ z_t)
+        r = R * acos(sind(ϕ_c) * sind(ϕ) + cosd(ϕ_c) * cosd(ϕ) * cosd(λ - λ_c))
+        c3 = cos(π * r / 2 / d_0)^3
+        s1 = sin(π * r / 2 / d_0)
+        cond = (0 < r < d_0) * (r != R * pi)
+        u +=
+            -16 * V_p / 3 / sqrt(FT(3)) *
+            F_z *
+            c3 *
+            s1 *
+            (-sind(ϕ_c) * cosd(ϕ) + cosd(ϕ_c) * sind(ϕ) * cosd(λ - λ_c)) /
+            sin(r / R) * cond
+        v +=
+            16 * V_p / 3 / sqrt(FT(3)) *
+            F_z *
+            c3 *
+            s1 *
+            cosd(ϕ_c) *
+            sind(λ - λ_c) / sin(r / R) * cond
+    end
+    return (; T_v, p, u, v)
+end
+
+function moist_baroclinic_wave_values(z, ϕ, λ, params, perturb, deep_atmosphere)
     FT = eltype(params)
     MSLP = CAP.MSLP(params)
 
@@ -430,7 +503,14 @@ function moist_baroclinic_wave_values(z, ϕ, λ, params, perturb)
     ϕ_w = FT(40)
     ε = FT(0.608)
 
-    (; T_v, p, u, v) = baroclinic_wave_values(z, ϕ, λ, params, perturb)
+    if deep_atmosphere
+        (; p, T_v, u, v) =
+            deep_atmos_baroclinic_wave_values(z, ϕ, λ, params, perturb)
+    else
+        (; p, T_v, u, v) =
+            shallow_atmos_baroclinic_wave_values(z, ϕ, λ, params, perturb)
+    end
+
     q_tot =
         (p <= p_t) ? q_t : q_0 * exp(-(ϕ / ϕ_w)^4) * exp(-((p - MSLP) / p_w)^2)
     T = T_v / (1 + ε * q_tot) # This is the formula used in the paper.
@@ -442,21 +522,34 @@ function moist_baroclinic_wave_values(z, ϕ, λ, params, perturb)
 end
 
 """
-    DryBaroclinicWave(; perturb = true)
+    DryBaroclinicWave(; perturb = true, deep_atmosphere = false)
 
 An `InitialCondition` with a dry baroclinic wave, and with an optional
 perturbation to the horizontal velocity.
 """
 Base.@kwdef struct DryBaroclinicWave <: InitialCondition
     perturb::Bool = true
+    deep_atmosphere::Bool = false
 end
 
 function (initial_condition::DryBaroclinicWave)(params)
-    (; perturb) = initial_condition
+    (; perturb, deep_atmosphere) = initial_condition
     function local_state(local_geometry)
         thermo_params = CAP.thermodynamics_params(params)
         (; z, lat, long) = local_geometry.coordinates
-        (; p, T_v, u, v) = baroclinic_wave_values(z, lat, long, params, perturb)
+        if deep_atmosphere
+            (; p, T_v, u, v) =
+                deep_atmos_baroclinic_wave_values(z, lat, long, params, perturb)
+        else
+            (; p, T_v, u, v) = shallow_atmos_baroclinic_wave_values(
+                z,
+                lat,
+                long,
+                params,
+                perturb,
+            )
+        end
+
         return LocalState(;
             params,
             geometry = local_geometry,
@@ -468,22 +561,29 @@ function (initial_condition::DryBaroclinicWave)(params)
 end
 
 """
-    MoistBaroclinicWave(; perturb = true)
+    MoistBaroclinicWave(; perturb = true, deep_atmosphere = false)
 
 An `InitialCondition` with a moist baroclinic wave, and with an optional
 perturbation to the horizontal velocity.
 """
 Base.@kwdef struct MoistBaroclinicWave <: InitialCondition
     perturb::Bool = true
+    deep_atmosphere::Bool = false
 end
 
 function (initial_condition::MoistBaroclinicWave)(params)
-    (; perturb) = initial_condition
+    (; perturb, deep_atmosphere) = initial_condition
     function local_state(local_geometry)
         thermo_params = CAP.thermodynamics_params(params)
         (; z, lat, long) = local_geometry.coordinates
-        (; p, T, q_tot, u, v) =
-            moist_baroclinic_wave_values(z, lat, long, params, perturb)
+        (; p, T, q_tot, u, v) = moist_baroclinic_wave_values(
+            z,
+            lat,
+            long,
+            params,
+            perturb,
+            deep_atmosphere,
+        )
         return LocalState(;
             params,
             geometry = local_geometry,
@@ -495,23 +595,30 @@ function (initial_condition::MoistBaroclinicWave)(params)
 end
 
 """
-    MoistBaroclinicWaveWithEDMF(; perturb = true)
+    MoistBaroclinicWaveWithEDMF(; perturb = true, deep_atmosphere = false)
 
 The same `InitialCondition` as `MoistBaroclinicWave`, except with an initial TKE
 of 0 and an initial draft area fraction of 0.2.
 """
 Base.@kwdef struct MoistBaroclinicWaveWithEDMF <: InitialCondition
     perturb::Bool = true
+    deep_atmosphere::Bool = false
 end
 
 function (initial_condition::MoistBaroclinicWaveWithEDMF)(params)
-    (; perturb) = initial_condition
+    (; perturb, deep_atmosphere) = initial_condition
     function local_state(local_geometry)
         FT = eltype(params)
         thermo_params = CAP.thermodynamics_params(params)
         (; z, lat, long) = local_geometry.coordinates
-        (; p, T, q_tot, u, v) =
-            moist_baroclinic_wave_values(z, lat, long, params, perturb)
+        (; p, T, q_tot, u, v) = moist_baroclinic_wave_values(
+            z,
+            lat,
+            long,
+            params,
+            perturb,
+            deep_atmosphere,
+        )
         return LocalState(;
             params,
             geometry = local_geometry,

src/solver/type_getters.jl

@@ -318,10 +318,14 @@ function get_initial_condition(parsed_args)
     if parsed_args["initial_condition"] in [
         "DryBaroclinicWave",
         "MoistBaroclinicWave",
-        "DecayingProfile",
         "MoistBaroclinicWaveWithEDMF",
-        "MoistAdiabaticProfileEDMFX",
     ]
+        return getproperty(ICs, Symbol(parsed_args["initial_condition"]))(
+            parsed_args["perturb_initstate"],
+            parsed_args["deep_atmosphere"],
+        )
+    elseif parsed_args["initial_condition"] in
+           ["DecayingProfile", "MoistAdiabaticProfileEDMFX"]
         return getproperty(ICs, Symbol(parsed_args["initial_condition"]))(
             parsed_args["perturb_initstate"],
         )