Add RCEMIPII analytic sounding initial condition

Author: oalcabes

.buildkite/pipeline.yml

@@ -168,6 +168,13 @@ steps:
         agents:
           slurm_mem: 20GB
 
+      - label: ":computer: CRM rcemipii in a box with 1M"
+        command: >
+          julia --color=yes --project=.buildkite .buildkite/ci_driver.jl
+          --config_file $CONFIG_PATH/rcemipii_box_CRM_1M.yml
+          --job_id rcemipii_box_CRM_1M
+        artifact_paths: "rcemipii_box_CRM_1M/output_active/*"
+
       - label: ":genie: LES ISDAC in a box"
         command: >
           julia --color=yes --project=.buildkite .buildkite/ci_driver.jl

config/model_configs/rcemipii_box_CRM_1M.yml

@@ -0,0 +1,87 @@
+reference_job_id: les_box  # for plotting
+surface_setup: DefaultMoninObukhov
+surface_temperature: RCEMIPII
+insolation: rcemipii
+initial_condition: RCEMIPIIProfile_300
+config: box
+rad: allskywithclear
+co2_model: fixed
+approximate_linear_solve_iters: 2  # only valid when implicit_diffusion=true
+rayleigh_sponge: true  # set sponge level in toml file
+smagorinsky_lilly: "UV"
+implicit_diffusion: true
+
+# Reference diffusion setup
+vert_diff: "DecayWithHeightDiffusion"
+hyperdiff: "false"
+
+# microphysics
+### !! 1M and 2M !!
+cloud_model: grid_scale   # Ultimately wanted when `moist: nonequil`
+moist: nonequil
+precip_model: 1M
+### spatial discretization ###
+x_max: 96000.0  # 96km
+y_max: 96000.0  # 96km
+# x_max: 6000000.0  # 6000km
+# y_max: 400000.0  # 400km
+z_max: 30000.0  # 30km
+nh_poly: 3  # hor. poly deg -> # quad pts in 1D in a h. elem is Nq = nh + 1
+# z_elem: 43  # Note: Set sponge height (`zd_viscous`) in toml file to cover top ~5 points
+z_elem: 86
+dz_bottom: 30.0
+
+#~8km grid
+x_elem: 4
+y_elem: 4
+dt: 10secs
+# ~4km grid
+# x_elem: 8
+# y_elem: 8
+# dt: 5secs
+# ~2km grid
+# x_elem: 16
+# y_elem: 16
+# dt: 2secs
+# ~1km grid
+# x_elem: 32
+# y_elem: 32
+# dt: 1secs
+
+t_end: 10hours
+# t_end: 2days
+# t_end: 1days
+# t_end: 21days
+
+## >> During debugging, reduce these as appropriate
+dt_save_state_to_disk: 1days
+dt_rad: 1hours  # same as `dt_save_state_to_disk` for reproducibility
+check_nan_every: 1024  # for debugging: set to check NaNs every `n` iterations
+log_to_file: true
+## <<
+
+toml: [toml/rcemipii_box.toml]
+
+# enable_diagnostics: false
+output_default_diagnostics: false  # adds various 1h diagnostics. Only use if `dt_save_state_to_disk` is a multiple of 1hr
+diagnostics:
+  - short_name: [
+      wa, ua, va, ta, thetaa, ha, # dynamics & thermodynamics
+      hus, hur, cl, clw, cli,  # liquid
+      pr,  # precipitation
+      ke,  # kinetic energy for spectrum
+      # Smagorinsky diagnostics
+      Dh_smag, strainh_smag,  # horizontal
+      # Dv_smag, strainv_smag,  # vertical
+      # DecayWithHeight diffusion coefficient
+      edt,
+    ]
+    # period: 10mins
+    period: 1hours
+  ## 1M microphsics
+  - short_name: [husra, hussn]
+    # period: 10mins
+    period: 1hours
+  # 2M microphysics
+  # - short_name: [cdnc, ncra]
+  #   period: 10mins

docs/src/api.md

@@ -41,6 +41,9 @@ ClimaAtmos.InitialConditions.TRMM_LBA
 ClimaAtmos.InitialConditions.LifeCycleTan2018
 ClimaAtmos.InitialConditions.Bomex
 ClimaAtmos.InitialConditions.Soares
+ClimaAtmos.InitialConditions.RCEMIPIIProfile_295
+ClimaAtmos.InitialConditions.RCEMIPIIProfile_300
+ClimaAtmos.InitialConditions.RCEMIPIIProfile_305
 ```
 
 ### Helper

src/initial_conditions/initial_conditions.jl

@@ -358,6 +358,171 @@ function (initial_condition::RisingThermalBubbleProfile)(params)
     return local_state
 end
 
+"""
+    RCEMIPIIProfile_295()
+
+An [`InitialCondition`](@ref) following the sounding to initialize simulations for
+RCEMIPII as described by Wing et. al. (2018)
+(https://doi.org/10.5194/gmd-11-793-2018). This profile is for an initial
+temperature of 295 K. Note: this should be used for RCE_small and NOT
+RCE_large - RCE_large must be initialized with the final state of RCE_small.
+"""
+
+struct RCEMIPIIProfile_295 <: InitialCondition end
+function (::RCEMIPIIProfile_295)(params)
+    function local_state(local_geometry)
+        FT = eltype(params)
+        R_d = CAP.R_d(params)
+        grav = CAP.grav(params)
+        thermo_params = CAP.thermodynamics_params(params)
+
+        T_0 = FT(295) # K
+        q_0 = FT(12e-3) # kg kg-1
+
+        q_t = FT(10^(-14)) # kg kg -1
+        z_q1 = FT(4000) # m
+        z_q2 = FT(7500) # m
+        z_t = FT(15000) # m
+        Γ = FT(0.0067) # K m-1
+        p_0 = FT(101480) # Pa
+
+        T_v0 = T_0 * (FT(1) + FT(0.608) * q_0)
+        T_vt = T_v0 - Γ * z_t
+
+        p_t = p_0 * (T_vt / T_v0)^(grav / (R_d * Γ))
+
+        RCEMIPII_T_v(z) = z ≤ z_t ? T_v0 - Γ * z : T_vt
+        RCEMIPII_p(z) =
+            z ≤ z_t ? p_0 * ((T_v0 - Γ * z) / T_v0)^(grav / (R_d * Γ)) :
+            p_t * exp(-grav * (z - z_t) / (R_d * T_vt))
+        RCEMIPII_q(z) = z ≤ z_t ? q_0 * exp(-z / z_q1) * exp(-(z / z_q2)^2) : q_t
+
+        (; z) = local_geometry.coordinates
+
+        q = RCEMIPII_q(z)
+        T_v = RCEMIPII_T_v(z)
+        T = T_v / (FT(1) + FT(0.608) * q)
+        p = RCEMIPII_p(z)
+
+        return LocalState(;
+            params,
+            geometry = local_geometry,
+            thermo_state = TD.PhaseEquil_pTq(thermo_params, p, T, q),
+        )
+    end
+    return local_state
+end
+
+"""
+    RCEMIPIIProfile_300()
+
+An [`InitialCondition`](@ref) following the sounding to initialize simulations for
+RCEMIPII as described by Wing et. al. (2018)
+(https://doi.org/10.5194/gmd-11-793-2018). This profile is for an initial
+temperature of 300 K. Note: this should be used for RCE_small and NOT
+RCE_large - RCE_large must be initialized with the final state of RCE_small.
+"""
+
+struct RCEMIPIIProfile_300 <: InitialCondition end
+function (::RCEMIPIIProfile_300)(params)
+    function local_state(local_geometry)
+        FT = eltype(params)
+        R_d = CAP.R_d(params)
+        grav = CAP.grav(params)
+        thermo_params = CAP.thermodynamics_params(params)
+
+        T_0 = FT(300) # K
+        q_0 = FT(18.65e-3) # kg kg-1
+
+        q_t = FT(10^(-14)) # kg kg -1
+        z_q1 = FT(4000) # m
+        z_q2 = FT(7500) # m
+        z_t = FT(15000) # m
+        Γ = FT(0.0067) # K m-1
+        p_0 = FT(101480) # Pa
+
+        T_v0 = T_0 * (FT(1) + FT(0.608) * q_0)
+        T_vt = T_v0 - Γ * z_t
+
+        p_t = p_0 * (T_vt / T_v0)^(grav / (R_d * Γ))
+
+        RCEMIPII_T_v(z) = z ≤ z_t ? T_v0 - Γ * z : T_vt
+        RCEMIPII_p(z) =
+            z ≤ z_t ? p_0 * ((T_v0 - Γ * z) / T_v0)^(grav / (R_d * Γ)) :
+            p_t * exp(-grav * (z - z_t) / (R_d * T_vt))
+        RCEMIPII_q(z) = z ≤ z_t ? q_0 * exp(-z / z_q1) * exp(-(z / z_q2)^2) : q_t
+
+        (; z) = local_geometry.coordinates
+
+        q = RCEMIPII_q(z)
+        T_v = RCEMIPII_T_v(z)
+        T = T_v / (FT(1) + FT(0.608) * q)
+        p = RCEMIPII_p(z)
+
+        return LocalState(;
+            params,
+            geometry = local_geometry,
+            thermo_state = TD.PhaseEquil_pTq(thermo_params, p, T, q),
+        )
+    end
+    return local_state
+end
+
+"""
+    RCEMIPIIProfile_305()
+
+An [`InitialCondition`](@ref) following the sounding to initialize simulations for
+RCEMIPII as described by Wing et. al. (2018)
+(https://doi.org/10.5194/gmd-11-793-2018). This profile is for an initial
+temperature of 305 K. Note: this should be used for RCE_small and NOT
+RCE_large - RCE_large must be initialized with the final state of RCE_small.
+"""
+
+struct RCEMIPIIProfile_305 <: InitialCondition end
+function (::RCEMIPIIProfile_305)(params)
+    function local_state(local_geometry)
+        FT = eltype(params)
+        R_d = CAP.R_d(params)
+        grav = CAP.grav(params)
+        thermo_params = CAP.thermodynamics_params(params)
+
+        T_0 = FT(305) # K
+        q_0 = FT(24e-3) # kg kg-1
+
+        q_t = FT(10^(-14)) # kg kg -1
+        z_q1 = FT(4000) # m
+        z_q2 = FT(7500) # m
+        z_t = FT(15000) # m
+        Γ = FT(0.0067) # K m-1
+        p_0 = FT(101480) # Pa
+
+        T_v0 = T_0 * (FT(1) + FT(0.608) * q_0)
+        T_vt = T_v0 - Γ * z_t
+
+        p_t = p_0 * (T_vt / T_v0)^(grav / (R_d * Γ))
+
+        RCEMIPII_T_v(z) = z ≤ z_t ? T_v0 - Γ * z : T_vt
+        RCEMIPII_p(z) =
+            z ≤ z_t ? p_0 * ((T_v0 - Γ * z) / T_v0)^(grav / (R_d * Γ)) :
+            p_t * exp(-grav * (z - z_t) / (R_d * T_vt))
+        RCEMIPII_q(z) = z ≤ z_t ? q_0 * exp(-z / z_q1) * exp(-(z / z_q2)^2) : q_t
+
+        (; z) = local_geometry.coordinates
+
+        q = RCEMIPII_q(z)
+        T_v = RCEMIPII_T_v(z)
+        T = T_v / (FT(1) + FT(0.608) * q)
+        p = RCEMIPII_p(z)
+
+        return LocalState(;
+            params,
+            geometry = local_geometry,
+            thermo_state = TD.PhaseEquil_pTq(thermo_params, p, T, q),
+        )
+    end
+    return local_state
+end
+
 """
     overwrite_initial_conditions!(initial_condition, args...)
 

src/solver/type_getters.jl

@@ -414,6 +414,9 @@ function get_initial_condition(parsed_args, atmos)
         "DryDensityCurrentProfile",
         "RisingThermalBubbleProfile",
         "PrecipitatingColumn",
+        "RCEMIPIIProfile_295",
+        "RCEMIPIIProfile_300",
+        "RCEMIPIIProfile_305",
     ]
         return getproperty(ICs, Symbol(parsed_args["initial_condition"]))()
     elseif isfile(parsed_args["initial_condition"])

toml/rcemipii_box.toml

@@ -0,0 +1,34 @@
+# toml/rcemipii_box.toml
+
+[c_smag]
+value = 1.0
+
+[condensation_evaporation_timescale]
+value = 50
+
+[sublimation_deposition_timescale]
+value = 50
+
+[angular_velocity_planet_rotation]
+value = 0
+
+[idealized_ocean_albedo]
+value = 0.07
+
+[mean_sea_level_pressure]
+value = 101480
+
+[D_0_diffusion]
+value = 5
+
+[H_diffusion]
+value = 800
+
+[SST_mean]
+value = 300
+
+[SST_delta]
+value = 1.25
+
+[SST_wavelength]
+value = 6e6