diff --git a/.buildkite/longruns/pipeline.yml b/.buildkite/longruns/pipeline.yml index 967ad2fb4c9..2b8e0fc1351 100644 --- a/.buildkite/longruns/pipeline.yml +++ b/.buildkite/longruns/pipeline.yml @@ -210,19 +210,6 @@ steps: env: JOB_NAME: "longrun_aquaplanet_rhoe_equil_highres_allsky_ft64" - - - group: "TurbulenceConvection" - steps: - - - label: ":balloon: Compressible single column EDMF TRMM_LBA" - command: "julia --color=yes --project=examples examples/hybrid/driver.jl --config_file $CONFIG_PATH/$$JOB_NAME.yml" - artifact_paths: "$$JOB_NAME/*" - agents: - slurm_mem: 20GB - slurm_time: 24:00:00 - env: - JOB_NAME: "longrun_compressible_edmf_trmm" - - group: "DYAMOND" steps: diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index 39c70da0c39..d88a300d459 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -359,22 +359,6 @@ steps: agents: slurm_mem: 20GB - - label: ":computer: aquaplanet (ρe_tot) equilmoist with EDMF" - command: > - julia --color=yes --project=examples examples/hybrid/driver.jl - --config_file $CONFIG_PATH/sphere_baroclinic_wave_rhoe_equilmoist_edmf.yml - - julia --color=yes --project=examples post_processing/remap/remap_pipeline.jl - --data_dir sphere_baroclinic_wave_rhoe_equilmoist_edmf - --out_dir sphere_baroclinic_wave_rhoe_equilmoist_edmf - - julia --color=yes --project=examples post_processing/plot/plot_pipeline.jl - --nc_dir sphere_baroclinic_wave_rhoe_equilmoist_edmf - --fig_dir sphere_baroclinic_wave_rhoe_equilmoist_edmf --case_name dry_held_suarez - artifact_paths: "sphere_baroclinic_wave_rhoe_equilmoist_edmf/*" - agents: - slurm_mem: 20GB - - label: ":computer: aquaplanet (ρe_tot) equilmoist allsky radiation monin_obukhov varying insolation gravity wave (raw_topo) high top zonally asymmetric" command: > julia --color=yes --project=examples examples/hybrid/driver.jl @@ -726,147 +710,6 @@ steps: agents: slurm_mem: 20GB - - group: "TurbulenceConvection" - steps: - - - label: ":balloon: SCM EDMF GABLS (Dry, JFNK IMEX EDMF)" - command: > - julia --color=yes --project=examples examples/hybrid/driver.jl - --config_file $CONFIG_PATH/edmf_gabls_jfnk_imex.yml - - julia --color=yes --project=examples regression_tests/test_mse.jl - --job_id edmf_gabls_jfnk_imex --out_dir edmf_gabls_jfnk_imex - artifact_paths: "edmf_gabls_jfnk_imex/*" - agents: - slurm_mem: 20GB - - - label: ":balloon: SCM EDMF: Nieuwstadt (Dry)" - command: > - julia --color=yes --project=examples examples/hybrid/driver.jl - --config_file $CONFIG_PATH/edmf_nieuwstadt.yml - - julia --color=yes --project=examples regression_tests/test_mse.jl - --job_id edmf_nieuwstadt --out_dir edmf_nieuwstadt - artifact_paths: "edmf_nieuwstadt/*" - agents: - slurm_mem: 20GB - - - label: ":balloon: SCM EDMF: Nieuwstadt (Dry, JFNK IMEX EDMF)" - command: > - julia --color=yes --project=examples examples/hybrid/driver.jl - --config_file $CONFIG_PATH/edmf_nieuwstadt_jfnk_imex.yml - - julia --color=yes --project=examples regression_tests/test_mse.jl - --job_id edmf_nieuwstadt_jfnk_imex --out_dir edmf_nieuwstadt_jfnk_imex - artifact_paths: "edmf_nieuwstadt_jfnk_imex/*" - soft_fail: true - agents: - slurm_mem: 20GB - - - label: ":balloon: SCM EDMF Soares" - command: > - julia --color=yes --project=examples examples/hybrid/driver.jl - --config_file $CONFIG_PATH/edmf_soares.yml - - julia --color=yes --project=examples regression_tests/test_mse.jl - --job_id edmf_soares --out_dir edmf_soares - artifact_paths: "edmf_soares/*" - agents: - slurm_mem: 20GB - - - label: ":balloon: SCM EDMF Soares (JFNK IMEX EDMF)" - command: > - julia --color=yes --project=examples examples/hybrid/driver.jl - --config_file $CONFIG_PATH/edmf_soares_jfnk_imex.yml - - julia --color=yes --project=examples regression_tests/test_mse.jl - --job_id edmf_soares_jfnk_imex --out_dir edmf_soares_jfnk_imex - artifact_paths: "edmf_soares_jfnk_imex/*" - soft_fail: true - agents: - slurm_mem: 20GB - - - label: ":balloon: SCM EDMF Bomex (Default: Newton's method)" - command: > - julia --color=yes --project=examples examples/hybrid/driver.jl - --config_file $CONFIG_PATH/edmf_bomex.yml - - julia --color=yes --project=examples regression_tests/test_mse.jl - --job_id edmf_bomex --out_dir edmf_bomex - artifact_paths: "edmf_bomex/*" - agents: - slurm_mem: 20GB - - - label: ":balloon: SCM EDMF Bomex (JFNK)" - command: > - julia --color=yes --project=examples examples/hybrid/driver.jl - --config_file $CONFIG_PATH/edmf_bomex_jfnk.yml - - julia --color=yes --project=examples regression_tests/test_mse.jl - --job_id edmf_bomex_jfnk --out_dir edmf_bomex_jfnk - artifact_paths: "edmf_bomex_jfnk/*" - agents: - slurm_mem: 20GB - - - label: ":balloon: SCM EDMF Bomex (JFNK IMEX EDMF)" - command: > - julia --color=yes --project=examples examples/hybrid/driver.jl - --config_file $CONFIG_PATH/edmf_bomex_jfnk_imex.yml - artifact_paths: "edmf_bomex_jfnk_imex/*" - agents: - slurm_mem: 20GB - - label: ":balloon: SCM EDMF DYCOMS_RF01" - command: > - julia --color=yes --project=examples examples/hybrid/driver.jl - --config_file $CONFIG_PATH/edmf_dycoms_rf01.yml - - julia --color=yes --project=examples regression_tests/test_mse.jl - --job_id edmf_dycoms_rf01 --out_dir edmf_dycoms_rf01 - artifact_paths: "edmf_dycoms_rf01/*" - agents: - slurm_mem: 20GB - - - label: ":balloon: SCM EDMF LifeCycleTan2018" - command: > - julia --color=yes --project=examples examples/hybrid/driver.jl - --config_file $CONFIG_PATH/edmf_life_cycle_tan2018.yml - - julia --color=yes --project=examples regression_tests/test_mse.jl - --job_id edmf_life_cycle_tan2018 --out_dir edmf_life_cycle_tan2018 - artifact_paths: "edmf_life_cycle_tan2018/*" - agents: - slurm_mem: 20GB - - - label: ":balloon: SCM EDMF Rico" - command: > - julia --color=yes --project=examples examples/hybrid/driver.jl - --config_file $CONFIG_PATH/edmf_rico.yml - - julia --color=yes --project=examples regression_tests/test_mse.jl - --job_id edmf_rico --out_dir edmf_rico - artifact_paths: "edmf_rico/*" - agents: - slurm_mem: 20GB - - - label: ":balloon: SCM EDMF TRMM_LBA (Short, 1-moment micro, sponge)" - command: > - julia --color=yes --project=examples examples/hybrid/driver.jl - --config_file $CONFIG_PATH/edmf_trmm.yml - - julia --color=yes --project=examples regression_tests/test_mse.jl - --job_id edmf_trmm --out_dir edmf_trmm - artifact_paths: "edmf_trmm/*" - agents: - slurm_mem: 20GB - - - label: ":dart: SCM EDMF Consistency test Bomex" - command: > - julia --color=yes --project=examples examples/hybrid/driver.jl - --config_file $CONFIG_PATH/edmf_bomex_consistency.yml - artifact_paths: "edmf_bomex_consistency/*" - agents: - slurm_mem: 20GB - - group: "GPU" steps: @@ -1004,14 +847,6 @@ steps: agents: slurm_mem: 20GB - - label: ":fire: Flame graph: perf target (edmf)" - command: > - julia --color=yes --project=perf perf/flame.jl - --config_file $PERF_CONFIG_PATH/flame/perf_target_edmf.yml - artifact_paths: "flame_perf_target_edmf/*" - agents: - slurm_mem: 20GB - - label: ":fire: Flame graph: perf target (Threaded)" command: > julia --threads 8 --color=yes --project=perf perf/flame.jl diff --git a/config/longrun_configs/longrun_compressible_edmf_trmm.yml b/config/longrun_configs/longrun_compressible_edmf_trmm.yml deleted file mode 100644 index 7fcaac0c345..00000000000 --- a/config/longrun_configs/longrun_compressible_edmf_trmm.yml +++ /dev/null @@ -1,17 +0,0 @@ -config: "column" -FLOAT_TYPE: "Float64" -hyperdiff: "false" -moist: "equil" -rad: "TRMM_LBA" -turbconv: "edmf" -turbconv_case: "TRMM_LBA" -precip_model: "1M" -dt_save_to_sol: "5mins" -z_elem: 90 -z_stretch: false -z_max: 17000 -rayleigh_sponge: true -job_id: "longrun_compressible_edmf_trmm" -dt: "0.2secs" -t_end: "3.8hours" -toml: [longrun_compressible_edmf_trmm.toml] diff --git a/config/model_configs/edmf_bomex.yml b/config/model_configs/edmf_bomex.yml deleted file mode 100644 index f0f29890aa7..00000000000 --- a/config/model_configs/edmf_bomex.yml +++ /dev/null @@ -1,20 +0,0 @@ -edmf_coriolis: Bomex -dt_save_to_disk: 5mins -hyperdiff: "false" -z_elem: 60 -dt: 20secs -debugging_tc: true -surface_setup: Bomex -turbconv_case: Bomex -t_end: 6hours -turbconv: "edmf" -z_stretch: false -config: column -subsidence: Bomex -FLOAT_TYPE: Float64 -z_max: 3000.0 -regression_test: true -ls_adv: Bomex -dt_save_to_sol: 5mins -job_id: edmf_bomex -moist: equil diff --git a/config/model_configs/edmf_bomex_consistency.yml b/config/model_configs/edmf_bomex_consistency.yml deleted file mode 100644 index 1b97ce4700b..00000000000 --- a/config/model_configs/edmf_bomex_consistency.yml +++ /dev/null @@ -1,20 +0,0 @@ -test_edmf_consistency: true -edmf_coriolis: "Bomex" -dt_save_to_disk: "5mins" -hyperdiff: "false" -z_elem: 60 -dt: "6secs" -debugging_tc: true -surface_setup: "Bomex" -turbconv_case: "Bomex" -t_end: "6hours" -turbconv: "edmf" -z_stretch: false -config: "column" -subsidence: "Bomex" -FLOAT_TYPE: "Float64" -z_max: 3000.0 -ls_adv: "Bomex" -dt_save_to_sol: "5mins" -job_id: "edmf_bomex_consistency" -moist: "equil" diff --git a/config/model_configs/edmf_bomex_jfnk.yml b/config/model_configs/edmf_bomex_jfnk.yml deleted file mode 100644 index 29a5af3a12d..00000000000 --- a/config/model_configs/edmf_bomex_jfnk.yml +++ /dev/null @@ -1,22 +0,0 @@ -edmf_coriolis: "Bomex" -dt_save_to_disk: "5mins" -hyperdiff: "false" -z_elem: 60 -dt: "15secs" -debugging_tc: true -surface_setup: "Bomex" -turbconv_case: "Bomex" -t_end: "6hours" -turbconv: "edmf" -z_stretch: false -config: "column" -subsidence: "Bomex" -use_krylov_method: true -FLOAT_TYPE: "Float64" -z_max: 3000.0 -reference_job_id: "edmf_bomex" -regression_test: true -ls_adv: "Bomex" -dt_save_to_sol: "5mins" -job_id: "edmf_bomex_jfnk" -moist: "equil" diff --git a/config/model_configs/edmf_bomex_jfnk_imex.yml b/config/model_configs/edmf_bomex_jfnk_imex.yml deleted file mode 100644 index e41c0ce32ad..00000000000 --- a/config/model_configs/edmf_bomex_jfnk_imex.yml +++ /dev/null @@ -1,23 +0,0 @@ -edmf_coriolis: "Bomex" -dt_save_to_disk: "5mins" -hyperdiff: "false" -z_elem: 60 -dt: "30secs" -imex_edmf_gm: true -surface_setup: "Bomex" -debugging_tc: true -turbconv_case: "Bomex" -t_end: "6hours" -turbconv: "edmf" -z_stretch: false -imex_edmf_turbconv: true -config: "column" -subsidence: "Bomex" -use_krylov_method: true -FLOAT_TYPE: "Float64" -z_max: 3000.0 -reference_job_id: "edmf_bomex" -ls_adv: "Bomex" -dt_save_to_sol: "5mins" -job_id: "edmf_bomex_jfnk_imex" -moist: "equil" diff --git a/config/model_configs/edmf_dycoms_rf01.yml b/config/model_configs/edmf_dycoms_rf01.yml deleted file mode 100644 index 6483a104dda..00000000000 --- a/config/model_configs/edmf_dycoms_rf01.yml +++ /dev/null @@ -1,20 +0,0 @@ -edmf_coriolis: "DYCOMS_RF01" -moist: "equil" -dt_save_to_disk: "2mins" -hyperdiff: "false" -z_elem: 30 -dt: "6secs" -debugging_tc: true -surface_setup: "DYCOMS_RF01" -turbconv_case: "DYCOMS_RF01" -t_end: "4hours" -turbconv: "edmf" -z_stretch: false -config: "column" -subsidence: "DYCOMS" -FLOAT_TYPE: "Float64" -z_max: 1500.0 -regression_test: true -dt_save_to_sol: "5mins" -job_id: "edmf_dycoms_rf01" -rad: "DYCOMS_RF01" diff --git a/config/model_configs/edmf_gabls_jfnk_imex.yml b/config/model_configs/edmf_gabls_jfnk_imex.yml deleted file mode 100644 index b59ef514bc5..00000000000 --- a/config/model_configs/edmf_gabls_jfnk_imex.yml +++ /dev/null @@ -1,21 +0,0 @@ -config: column -turbconv: edmf -debugging_tc: true -regression_test: true -FLOAT_TYPE: Float64 -hyperdiff: "false" -turbconv_case: GABLS -edmf_coriolis: GABLS -surface_setup: GABLS -z_elem: 8 -z_stretch: false -z_max: 400 -t_end: 9hours -dt_save_to_sol: 5mins -dt_save_to_disk: 5mins -moist: equil -imex_edmf_turbconv: true -imex_edmf_gm: true -use_krylov_method: true -job_id: edmf_gabls_jfnk_imex -dt: 600secs diff --git a/config/model_configs/edmf_life_cycle_tan2018.yml b/config/model_configs/edmf_life_cycle_tan2018.yml deleted file mode 100644 index 27d85950d5f..00000000000 --- a/config/model_configs/edmf_life_cycle_tan2018.yml +++ /dev/null @@ -1,20 +0,0 @@ -edmf_coriolis: "LifeCycleTan2018" -dt_save_to_disk: "5mins" -hyperdiff: "false" -z_elem: 75 -dt: "6secs" -debugging_tc: true -surface_setup: "LifeCycleTan2018" -turbconv_case: "LifeCycleTan2018" -t_end: "6hours" -turbconv: "edmf" -z_stretch: false -config: "column" -subsidence: "LifeCycleTan2018" -FLOAT_TYPE: "Float64" -z_max: 3000.0 -regression_test: true -ls_adv: "LifeCycleTan2018" -dt_save_to_sol: "5mins" -job_id: "edmf_life_cycle_tan2018" -moist: "equil" diff --git a/config/model_configs/edmf_nieuwstadt.yml b/config/model_configs/edmf_nieuwstadt.yml deleted file mode 100644 index 227007b0c4d..00000000000 --- a/config/model_configs/edmf_nieuwstadt.yml +++ /dev/null @@ -1,17 +0,0 @@ -config: column -turbconv: edmf -debugging_tc: true -regression_test: true -FLOAT_TYPE: Float64 -hyperdiff: "false" -turbconv_case: Nieuwstadt -surface_setup: Nieuwstadt -z_elem: 200 -z_stretch: false -z_max: 3750 -t_end: 8hours -dt_save_to_sol: 5mins -dt_save_to_disk: 5mins -moist: equil -job_id: edmf_nieuwstadt -dt: 0.5secs diff --git a/config/model_configs/edmf_nieuwstadt_jfnk_imex.yml b/config/model_configs/edmf_nieuwstadt_jfnk_imex.yml deleted file mode 100644 index 9cce2c9f70b..00000000000 --- a/config/model_configs/edmf_nieuwstadt_jfnk_imex.yml +++ /dev/null @@ -1,20 +0,0 @@ -config: column -turbconv: edmf -debugging_tc: true -FLOAT_TYPE: Float64 -hyperdiff: "false" -turbconv_case: Nieuwstadt -surface_setup: Nieuwstadt -z_elem: 200 -z_stretch: false -z_max: 3750 -t_end: 8hours -dt_save_to_sol: 5mins -dt_save_to_disk: 5mins -moist: equil -imex_edmf_turbconv: true -imex_edmf_gm: true -use_krylov_method: true -reference_job_id: edmf_nieuwstadt -job_id: edmf_nieuwstadt_jfnk_imex -dt: 50secs diff --git a/config/model_configs/edmf_rico.yml b/config/model_configs/edmf_rico.yml deleted file mode 100644 index 5feaca6d993..00000000000 --- a/config/model_configs/edmf_rico.yml +++ /dev/null @@ -1,20 +0,0 @@ -dt_save_to_disk: "8mins" -hyperdiff: "false" -z_elem: 100 -dt: "0.75secs" -debugging_tc: true -surface_setup: "Rico" -turbconv_case: "Rico" -t_end: "24hours" -turbconv: "edmf" -z_stretch: false -config: "column" -subsidence: "Rico" -FLOAT_TYPE: "Float64" -z_max: 4000.0 -precip_model: "1M" -regression_test: true -ls_adv: "Rico" -dt_save_to_sol: "5mins" -job_id: "edmf_rico" -moist: "equil" diff --git a/config/model_configs/edmf_soares.yml b/config/model_configs/edmf_soares.yml deleted file mode 100644 index 0697164f827..00000000000 --- a/config/model_configs/edmf_soares.yml +++ /dev/null @@ -1,17 +0,0 @@ -dt_save_to_disk: "5mins" -hyperdiff: "false" -z_elem: 150 -dt: "0.25secs" -debugging_tc: true -surface_setup: "Soares" -turbconv_case: "Soares" -t_end: "8hours" -turbconv: "edmf" -z_stretch: false -config: "column" -FLOAT_TYPE: "Float64" -z_max: 3750.0 -regression_test: true -dt_save_to_sol: "5mins" -job_id: "edmf_soares" -moist: "equil" diff --git a/config/model_configs/edmf_soares_jfnk_imex.yml b/config/model_configs/edmf_soares_jfnk_imex.yml deleted file mode 100644 index 1b67953b24b..00000000000 --- a/config/model_configs/edmf_soares_jfnk_imex.yml +++ /dev/null @@ -1,20 +0,0 @@ -dt_save_to_disk: "5mins" -hyperdiff: "false" -z_elem: 150 -dt: "50secs" -imex_edmf_gm: true -surface_setup: "Soares" -debugging_tc: true -turbconv_case: "Soares" -t_end: "8hours" -turbconv: "edmf" -z_stretch: false -imex_edmf_turbconv: true -config: "column" -use_krylov_method: true -FLOAT_TYPE: "Float64" -z_max: 3750.0 -reference_job_id: "edmf_soares" -dt_save_to_sol: "5mins" -job_id: "edmf_soares_jfnk_imex" -moist: "equil" diff --git a/config/model_configs/edmf_trmm.yml b/config/model_configs/edmf_trmm.yml deleted file mode 100644 index 8e0d1eaa028..00000000000 --- a/config/model_configs/edmf_trmm.yml +++ /dev/null @@ -1,21 +0,0 @@ -moist: "equil" -dt_save_to_disk: "1mins" -rayleigh_sponge: true -hyperdiff: "false" -z_elem: 82 -dt: "0.1secs" -debugging_tc: true -surface_setup: "TRMM_LBA" -turbconv_case: "TRMM_LBA" -t_end: "1hours" -turbconv: "edmf" -z_stretch: false -config: "column" -FLOAT_TYPE: "Float64" -z_max: 16400.0 -precip_model: "1M" -regression_test: true -dt_save_to_sol: "5mins" -job_id: "edmf_trmm" -rad: "TRMM_LBA" -toml: [toml/edmf_trmm.toml] diff --git a/config/model_configs/sphere_baroclinic_wave_rhoe_equilmoist_edmf.yml b/config/model_configs/sphere_baroclinic_wave_rhoe_equilmoist_edmf.yml deleted file mode 100644 index 22902af9f2f..00000000000 --- a/config/model_configs/sphere_baroclinic_wave_rhoe_equilmoist_edmf.yml +++ /dev/null @@ -1,13 +0,0 @@ -job_id: "sphere_baroclinic_wave_rhoe_equilmoist_edmf" -vert_diff: "true" -surface_setup: "DefaultExchangeCoefficients" -moist: "dry" -dt: "1secs" -t_end: "20secs" -FLOAT_TYPE: "Float64" -turbconv: "edmf" -config: "sphere" -z_stretch: false -post_process: false -dt_save_to_sol: "10secs" -dt_save_to_disk: "10secs" \ No newline at end of file diff --git a/config/perf_configs/flame/perf_target_edmf.yml b/config/perf_configs/flame/perf_target_edmf.yml deleted file mode 100644 index 8dd6088c683..00000000000 --- a/config/perf_configs/flame/perf_target_edmf.yml +++ /dev/null @@ -1,5 +0,0 @@ -job_id: "flame_perf_target_edmf" -precip_model: "nothing" -target_job: "sphere_baroclinic_wave_rhoe_equilmoist_edmf" -rad: "nothing" -apply_limiter: false \ No newline at end of file diff --git a/config/perf_configs/jet_n_failures.yml b/config/perf_configs/jet_n_failures.yml index 6c8f6142a93..36b0100ca6d 100644 --- a/config/perf_configs/jet_n_failures.yml +++ b/config/perf_configs/jet_n_failures.yml @@ -1,4 +1,4 @@ -target_job: "sphere_baroclinic_wave_rhoe_equilmoist_edmf" +target_job: "sphere_baroclinic_wave_rhoe_equilmoist_edmfx" moist: "dry" precip_model: ~ rad: ~ diff --git a/examples/hybrid/driver.jl b/examples/hybrid/driver.jl index 7e9db267861..a8d769ae0de 100644 --- a/examples/hybrid/driver.jl +++ b/examples/hybrid/driver.jl @@ -119,14 +119,8 @@ if !CA.is_distributed(config.comms_ctx) && !is_edmfx && !(atmos.model_config isa CA.SphericalModel) ENV["GKSwstype"] = "nul" # avoid displaying plots - if CA.is_column_without_edmf(config.parsed_args) + if CA.is_column_without_edmfx(config.parsed_args) custom_postprocessing(sol, simulation.output_dir, p) - elseif CA.is_column_edmf(config.parsed_args) - postprocessing_edmf( - sol, - simulation.output_dir, - config.parsed_args["fps"], - ) elseif CA.is_solid_body(config.parsed_args) postprocessing(sol, simulation.output_dir, config.parsed_args["fps"]) elseif atmos.model_config isa CA.BoxModel diff --git a/perf/flame.jl b/perf/flame.jl index 80ab8d826be..5d69e2f0d0e 100644 --- a/perf/flame.jl +++ b/perf/flame.jl @@ -59,7 +59,6 @@ allocs_limit["flame_perf_target_tracers"] = 212496 allocs_limit["flame_perf_target_edmfx"] = 304064 allocs_limit["flame_perf_diagnostics"] = 3024344 allocs_limit["flame_perf_target_diagnostic_edmfx"] = 762784 -allocs_limit["flame_perf_target_edmf"] = 12459299664 allocs_limit["flame_perf_target_threaded"] = 6175664 allocs_limit["flame_perf_target_callbacks"] = 46413904 allocs_limit["flame_perf_gw"] = 4911463328 diff --git a/post_processing/post_processing_funcs.jl b/post_processing/post_processing_funcs.jl index 75ac542c624..2fcad3b2ee2 100644 --- a/post_processing/post_processing_funcs.jl +++ b/post_processing/post_processing_funcs.jl @@ -205,10 +205,6 @@ function safe_index(ius, t) end end -function postprocessing_edmf(sol, output_dir, fps) - profile_animation(sol, output_dir, fps) -end - function custom_postprocessing(sol, output_dir, p) thermo_params = CAP.thermodynamics_params(params) get_var(i, var) = Fields.single_field(sol.u[i], var) diff --git a/src/ClimaAtmos.jl b/src/ClimaAtmos.jl index 1f0c7fcac38..9f7a5795fa4 100644 --- a/src/ClimaAtmos.jl +++ b/src/ClimaAtmos.jl @@ -22,7 +22,7 @@ include(joinpath("parameterized_tendencies", "radiation", "RRTMGPInterface.jl")) import .RRTMGPInterface as RRTMGPI include(joinpath("parameterized_tendencies", "radiation", "radiation.jl")) -include(joinpath("TurbulenceConvection_deprecated", "TurbulenceConvection.jl")) +include(joinpath("TurbulenceConvection", "TurbulenceConvection.jl")) import .TurbulenceConvection as TC include(joinpath("cache", "edmf_precomputed_quantities.jl")) @@ -101,7 +101,6 @@ include( include(joinpath("parameterized_tendencies", "sponge", "rayleigh_sponge.jl")) include(joinpath("parameterized_tendencies", "sponge", "viscous_sponge.jl")) include(joinpath("prognostic_equations", "advection.jl")) -include(joinpath("dycore_equations_deprecated", "sgs_flux_tendencies.jl")) include(joinpath("cache", "temporary_quantities.jl")) include(joinpath("cache", "cache.jl")) diff --git a/src/TurbulenceConvection_deprecated/Parameters.jl b/src/TurbulenceConvection/Parameters.jl similarity index 100% rename from src/TurbulenceConvection_deprecated/Parameters.jl rename to src/TurbulenceConvection/Parameters.jl diff --git a/src/TurbulenceConvection/TurbulenceConvection.jl b/src/TurbulenceConvection/TurbulenceConvection.jl new file mode 100644 index 00000000000..66ae8d3726c --- /dev/null +++ b/src/TurbulenceConvection/TurbulenceConvection.jl @@ -0,0 +1,9 @@ +module TurbulenceConvection + +include("Parameters.jl") +import .Parameters as TCP +const APS = TCP.AbstractTurbulenceConvectionParameters + +Base.broadcastable(param_set::APS) = tuple(param_set) + +end diff --git a/src/TurbulenceConvection_deprecated/EDMF_functions.jl b/src/TurbulenceConvection_deprecated/EDMF_functions.jl deleted file mode 100644 index 9692e843cb1..00000000000 --- a/src/TurbulenceConvection_deprecated/EDMF_functions.jl +++ /dev/null @@ -1,633 +0,0 @@ -function update_cloud_frac(edmf::EDMFModel, state::State) - # update grid-mean cloud fraction and cloud cover - aux_bulk = center_aux_bulk(state) - aux_gm = center_aux_grid_mean(state) - aux_en = center_aux_environment(state) - a_up_bulk = aux_bulk.area - # update grid-mean cloud fraction and cloud cover - @. aux_gm.cloud_fraction = - aux_en.area * aux_en.cloud_fraction + - a_up_bulk * aux_bulk.cloud_fraction -end - -function compute_implicit_turbconv_tendencies!( - edmf::EDMFModel, - grid::Grid, - state::State, -) - compute_implicit_up_tendencies!(edmf, grid, state) - return nothing -end - -function compute_explicit_turbconv_tendencies!( - edmf::EDMFModel, - grid::Grid, - state::State, -) - compute_explicit_up_tendencies!(edmf, grid, state) - return nothing -end - -function compute_sgs_flux!( - edmf::EDMFModel, - grid::Grid, - state::State, - param_set::APS, -) - (; surface_conditions) = state - thermo_params = TCP.thermodynamics_params(param_set) - N_up = n_updrafts(edmf) - FT = float_type(state) - prog_gm = center_prog_grid_mean(state) - aux_gm = center_aux_grid_mean(state) - prog_gm_f = face_prog_grid_mean(state) - aux_gm_f = face_aux_grid_mean(state) - aux_en = center_aux_environment(state) - aux_en_f = face_aux_environment(state) - aux_up = center_aux_updrafts(state) - aux_tc_f = face_aux_turbconv(state) - aux_up_f = face_aux_updrafts(state) - prog_up_f = face_prog_updrafts(state) - ρ_f = aux_gm_f.ρ - ρ_c = prog_gm.ρ - p_c = center_aux_grid_mean_p(state) - kf_surf = kf_surface(grid) - kc_surf = kc_surface(grid) - massflux = aux_tc_f.massflux - massflux_h = aux_tc_f.massflux_h - massflux_qt = aux_tc_f.massflux_qt - aux_tc = center_aux_turbconv(state) - - wvec = CC.Geometry.WVector - ∇c = CCO.DivergenceF2C() - - # TODO: we shouldn't need to call parent here - a_en = aux_en.area - w_en = aux_en_f.w - w_gm = prog_gm_f.u₃ - h_tot_gm = aux_gm.h_tot - q_tot_gm = aux_gm.q_tot - a_en_bcs = (; - bottom = CCO.SetValue( - 1 - sum( - i -> area_surface_bc(surface_conditions, edmf), - 1:n_updrafts(edmf), - ), - ), - top = CCO.Extrapolate(), - ) - Ifae = CCO.InterpolateC2F(; a_en_bcs...) - If = CCO.InterpolateC2F(; - bottom = CCO.SetValue(FT(0)), - top = CCO.SetValue(FT(0)), - ) - Ic = CCO.InterpolateF2C() - - # compute total enthalpies - ts_en = center_aux_environment(state).ts - ts_gm = center_aux_grid_mean_ts(state) - @. h_tot_gm = - TD.total_specific_enthalpy(thermo_params, ts_gm, prog_gm.ρe_tot / ρ_c) - # Compute the mass flux and associated scalar fluxes - @. massflux = ρ_f * Ifae(a_en) * (w_en - w_gm) - @. massflux_h = - ρ_f * Ifae(a_en) * (w_en - w_gm) * (If(aux_en.h_tot) - If(h_tot_gm)) - @. massflux_qt = - ρ_f * Ifae(a_en) * (w_en - w_gm) * (If(aux_en.q_tot) - If(q_tot_gm)) - @inbounds for i in 1:N_up - a_up_bcs = a_up_boundary_conditions(surface_conditions, edmf) - Ifau = CCO.InterpolateC2F(; a_up_bcs...) - a_up = aux_up[i].area - w_up_i = prog_up_f[i].w - @. aux_up_f[i].massflux = ρ_f * Ifau(a_up) * (w_up_i - w_gm) - @. massflux_h += - ρ_f * ( - Ifau(a_up) * - (w_up_i - w_gm) * - (If(aux_up[i].h_tot) - If(h_tot_gm)) - ) - @. massflux_qt += - ρ_f * ( - Ifau(a_up) * - (w_up_i - w_gm) * - (If(aux_up[i].q_tot) - If(q_tot_gm)) - ) - end - - massflux_h[kf_surf] = zero(eltype(massflux_h)) - massflux_qt[kf_surf] = zero(eltype(massflux_qt)) - - diffusive_flux_h = aux_tc_f.diffusive_flux_h - diffusive_flux_qt = aux_tc_f.diffusive_flux_qt - diffusive_flux_uₕ = aux_tc_f.diffusive_flux_uₕ - - sgs_flux_h_tot = aux_gm_f.sgs_flux_h_tot - sgs_flux_q_tot = aux_gm_f.sgs_flux_q_tot - sgs_flux_uₕ = aux_gm_f.sgs_flux_uₕ - - @. sgs_flux_h_tot = diffusive_flux_h + massflux_h - @. sgs_flux_q_tot = diffusive_flux_qt + massflux_qt - @. sgs_flux_uₕ = diffusive_flux_uₕ # + massflux_u - - # apply surface BC as SGS flux at lowest level - sgs_flux_h_tot[kf_surf] = surface_conditions.ρ_flux_h_tot - sgs_flux_q_tot[kf_surf] = surface_conditions.ρ_flux_q_tot - sgs_flux_uₕ[kf_surf] = - edmf.zero_uv_fluxes ? zero(surface_conditions.ρ_flux_uₕ) : - surface_conditions.ρ_flux_uₕ - return nothing -end - -function compute_diffusive_fluxes( - edmf::EDMFModel, - grid::Grid, - state::State, - param_set::APS, -) - (; surface_conditions) = state - thermo_params = TCP.thermodynamics_params(param_set) - FT = float_type(state) - aux_bulk = center_aux_bulk(state) - aux_tc_f = face_aux_turbconv(state) - aux_en_f = face_aux_environment(state) - aux_en = center_aux_environment(state) - aux_gm = center_aux_grid_mean(state) - aux_gm_f = face_aux_grid_mean(state) - KM = center_aux_turbconv(state).KM - KH = center_aux_turbconv(state).KH - aeKM = center_aux_turbconv(state).ϕ_temporary - aeKH = center_aux_turbconv(state).ψ_temporary - prog_gm_uₕ = grid_mean_uₕ(state) - - ρ_f = aux_gm_f.ρ - a_en = aux_en.area - @. aeKM = a_en * KM - @. aeKH = a_en * KH - kc_surf = kc_surface(grid) - kc_toa = kc_top_of_atmos(grid) - kf_surf = kf_surface(grid) - prog_gm = center_prog_grid_mean(state) - ts_gm = center_aux_grid_mean_ts(state) - IfKH = CCO.InterpolateC2F(; - bottom = CCO.SetValue(aeKH[kc_surf]), - top = CCO.SetValue(aeKH[kc_toa]), - ) - IfKM = CCO.InterpolateC2F(; - bottom = CCO.SetValue(aeKM[kc_surf]), - top = CCO.SetValue(aeKM[kc_toa]), - ) - Ic = CCO.InterpolateF2C() - - @. aux_tc_f.ρ_ae_KH = IfKH(aeKH) * ρ_f - @. aux_tc_f.ρ_ae_KM = IfKM(aeKM) * ρ_f - - aeKHq_tot_bc = - -surface_conditions.ρ_flux_q_tot / a_en[kc_surf] / - aux_tc_f.ρ_ae_KH[kf_surf] - aeKHh_tot_bc = - -surface_conditions.ρ_flux_h_tot / a_en[kc_surf] / - aux_tc_f.ρ_ae_KH[kf_surf] - ∇q_tot_en = CCO.GradientC2F(; - bottom = CCO.SetGradient(aeKHq_tot_bc), - top = CCO.SetGradient(zero(aeKHq_tot_bc)), - ) - ∇h_tot_en = CCO.GradientC2F(; - bottom = CCO.SetGradient(aeKHh_tot_bc), - top = CCO.SetGradient(zero(aeKHh_tot_bc)), - ) - ρ_flux_uₕ = - edmf.zero_uv_fluxes ? zero(surface_conditions.ρ_flux_uₕ) : - surface_conditions.ρ_flux_uₕ - ∇uₕ_gm = CCO.GradientC2F(; - bottom = CCO.SetGradient( - -ρ_flux_uₕ / a_en[kc_surf] / aux_tc_f.ρ_ae_KM[kf_surf], - ), - top = CCO.SetGradient(zero(ρ_flux_uₕ)), - ) - - @. aux_tc_f.diffusive_flux_qt = -aux_tc_f.ρ_ae_KH * ∇q_tot_en(aux_en.q_tot) - @. aux_tc_f.diffusive_flux_h = -aux_tc_f.ρ_ae_KH * ∇h_tot_en(aux_en.h_tot) - @. aux_tc_f.diffusive_flux_uₕ = -aux_tc_f.ρ_ae_KM * ∇uₕ_gm(prog_gm_uₕ) - - return nothing -end - -function affect_filter!( - edmf::EDMFModel, - grid::Grid, - state::State, - param_set::APS, - t::Real, -) - ### - ### Filters - ### - set_edmf_surface_bc(edmf, grid, state, param_set) - filter_updraft_vars(edmf, grid, state, param_set) - return nothing -end - -function set_edmf_surface_bc( - edmf::EDMFModel, - grid::Grid, - state::State, - param_set::APS, -) - (; surface_conditions) = state - thermo_params = TCP.thermodynamics_params(param_set) - FT = float_type(state) - N_up = n_updrafts(edmf) - kc_surf = kc_surface(grid) - kf_surf = kf_surface(grid) - prog_gm = center_prog_grid_mean(state) - prog_gm_uₕ = grid_mean_uₕ(state) - prog_up = center_prog_updrafts(state) - prog_up_f = face_prog_updrafts(state) - ts_gm = center_aux_grid_mean_ts(state) - ρ_c = prog_gm.ρ - p_c = TD.air_pressure.(thermo_params, ts_gm) - C123 = CCG.Covariant123Vector - Ic = CCO.InterpolateF2C() - @inbounds for i in 1:N_up - θ_surf = θ_surface_bc(grid, state, edmf, i, param_set) - q_surf = q_surface_bc(grid, state, edmf, i, param_set) - a_surf = area_surface_bc(surface_conditions, edmf) - e_kin = @. LA.norm_sqr( - C123(prog_gm_uₕ) + C123(Ic(CCG.WVector(prog_up_f[i].w))), - ) / 2 - e_pot_surf = geopotential(thermo_params, grid.zc.z[kc_surf]) - ts_up_i_surf = - TD.PhaseEquil_pθq(thermo_params, p_c[kc_surf], θ_surf, q_surf) - e_tot_surf = TD.total_energy( - thermo_params, - ts_up_i_surf, - e_kin[kc_surf], - e_pot_surf, - ) - prog_up[i].ρarea[kc_surf] = ρ_c[kc_surf] * a_surf - prog_up[i].ρae_tot[kc_surf] = prog_up[i].ρarea[kc_surf] * e_tot_surf - prog_up[i].ρaq_tot[kc_surf] = prog_up[i].ρarea[kc_surf] * q_surf - prog_up_f[i].w[kf_surf] = CCG.Covariant3Vector( - CCG.WVector(FT(0)), - CC.Fields.local_geometry_field(axes(prog_up_f))[kf_surf], - ) - end - return nothing -end - -function surface_helper(grid::Grid, state::State) - (; surface_conditions) = state - FT = float_type(state) - zLL::FT = grid.zc[Cent(1)].z - ρLL = center_prog_grid_mean(state).ρ[Cent(1)] - return (; - ustar = surface_conditions.ustar, - zLL, - oblength = surface_conditions.obukhov_length, - ρLL, - ) -end - -function a_up_boundary_conditions(surface_conditions, edmf::EDMFModel) - a_surf = area_surface_bc(surface_conditions, edmf) - return (; bottom = CCO.SetValue(a_surf), top = CCO.Extrapolate()) -end - -function area_surface_bc(surface_conditions, edmf::EDMFModel) - FT = typeof(edmf.minimum_area) - return surface_conditions.buoyancy_flux > 0 ? - edmf.surface_area / n_updrafts(edmf) : FT(0) -end - -function θ_surface_bc( - grid::Grid, - state::State, - edmf::EDMFModel, - i::Int, - param_set::APS, -) - FT = eltype(edmf) - (; surface_conditions) = state - thermo_params = TCP.thermodynamics_params(param_set) - aux_gm = center_aux_grid_mean(state) - kc_surf = kc_surface(grid) - ts_gm = center_aux_grid_mean_ts(state) - c_p = TD.cp_m(thermo_params, ts_gm[kc_surf]) - (; ustar, zLL, oblength, ρLL) = surface_helper(grid, state) - surface_conditions.buoyancy_flux > 0 || return aux_gm.θ_liq_ice[kc_surf] - a_total = edmf.surface_area - a_ = area_surface_bc(surface_conditions, edmf) - ρθ_liq_ice_flux = surface_conditions.ρ_flux_θ - h_var = get_surface_variance( - ρθ_liq_ice_flux / ρLL, - ustar, - zLL, - oblength, - Fields.local_geometry_field(Fields.level(grid.fs, Fields.half))[], - ) - surface_scalar_coeff = percentile_bounds_mean_norm( - 1 - a_total + (i - 1) * a_, - 1 - a_total + i * a_, - ) - return aux_gm.θ_liq_ice[kc_surf] + surface_scalar_coeff * sqrt(h_var) -end - -function q_surface_bc( - grid::Grid, - state::State, - edmf::EDMFModel, - i::Int, - param_set, -) - (; surface_conditions) = state - aux_gm = center_aux_grid_mean(state) - kc_surf = kc_surface(grid) - surface_conditions.buoyancy_flux > 0 || return aux_gm.q_tot[kc_surf] - a_total = edmf.surface_area - a_ = area_surface_bc(surface_conditions, edmf) - (; ustar, zLL, oblength, ρLL) = surface_helper(grid, state) - ρq_tot_flux = surface_conditions.ρ_flux_q_tot - qt_var = get_surface_variance( - ρq_tot_flux / ρLL, - ustar, - zLL, - oblength, - Fields.local_geometry_field(Fields.level(grid.fs, Fields.half))[], - ) - surface_scalar_coeff = percentile_bounds_mean_norm( - 1 - a_total + (i - 1) * a_, - 1 - a_total + i * a_, - ) - return aux_gm.q_tot[kc_surf] + surface_scalar_coeff * sqrt(qt_var) -end - -function compute_implicit_up_tendencies!( - edmf::EDMFModel, - grid::Grid, - state::State, -) - N_up = n_updrafts(edmf) - kc_surf = kc_surface(grid) - kf_surf = kf_surface(grid) - FT = float_type(state) - - prog_up = center_prog_updrafts(state) - prog_up_f = face_prog_updrafts(state) - tendencies_up = center_tendencies_updrafts(state) - tendencies_up_f = face_tendencies_updrafts(state) - aux_up = center_aux_updrafts(state) - ts = center_aux_grid_mean_ts(state) - p_c = center_aux_grid_mean_p(state) - prog_gm = center_prog_grid_mean(state) - ρ_c = prog_gm.ρ - - # Solve for updraft area fraction - - Ic = CCO.InterpolateF2C() - ∇c = CCO.DivergenceF2C() - LBF = CCO.LeftBiasedC2F(; bottom = CCO.SetValue(CCG.WVector(FT(0)))) - - @inbounds for i in 1:N_up - w_up = prog_up_f[i].w - - ρarea = prog_up[i].ρarea - ρae_tot = prog_up[i].ρae_tot - ρaq_tot = prog_up[i].ρaq_tot - - tends_ρarea = tendencies_up[i].ρarea - tends_ρae_tot = tendencies_up[i].ρae_tot - tends_ρaq_tot = tendencies_up[i].ρaq_tot - - volume_term = - @. -p_c / ρ_c * (-(∇c(LBF(Ic(CCG.WVector(w_up)) * ρarea)))) - @. tends_ρarea += -∇c(LBF(Ic(CCG.WVector(w_up)) * ρarea)) - @. tends_ρae_tot += - -∇c(LBF(Ic(CCG.WVector(w_up)) * ρarea * aux_up[i].h_tot)) + - volume_term - @. tends_ρaq_tot += -∇c(LBF(Ic(CCG.WVector(w_up)) * ρaq_tot)) - - tends_ρarea[kc_surf] = 0 - tends_ρae_tot[kc_surf] = 0 - tends_ρaq_tot[kc_surf] = 0 - end - - # Solve for updraft velocity - - LBC = CCO.LeftBiasedF2C(; bottom = CCO.SetValue(FT(0))) - prog_bcs = (; - bottom = CCO.SetGradient(CCG.WVector(FT(0))), - top = CCO.SetGradient(CCG.WVector(FT(0))), - ) - grad_f = CCO.GradientC2F(; prog_bcs...) - - @inbounds for i in 1:N_up - w_up = prog_up_f[i].w - tends_w = tendencies_up_f[i].w - @. tends_w += -grad_f(LBC(LA.norm_sqr(CCG.WVector(w_up)) / 2)) - tends_w[kf_surf] = zero(tends_w[kf_surf]) - end - - return nothing -end - -function compute_explicit_up_tendencies!( - edmf::EDMFModel, - grid::Grid, - state::State, -) - N_up = n_updrafts(edmf) - kc_surf = kc_surface(grid) - kf_surf = kf_surface(grid) - FT = float_type(state) - - aux_up = center_aux_updrafts(state) - aux_en = center_aux_environment(state) - aux_en_f = face_aux_environment(state) - aux_up_f = face_aux_updrafts(state) - prog_up = center_prog_updrafts(state) - prog_up_f = face_prog_updrafts(state) - tendencies_up = center_tendencies_updrafts(state) - tendencies_up_f = face_tendencies_updrafts(state) - prog_gm = center_prog_grid_mean(state) - p_c = center_aux_grid_mean_p(state) - ρ_c = prog_gm.ρ - - Ic = CCO.InterpolateF2C() - # We know that, since W = 0 at z = 0, BCs for entr, detr, - # and buoyancy should not matter in the end - zero_bcs = (; bottom = CCO.SetValue(FT(0)), top = CCO.SetValue(FT(0))) - I0f = CCO.InterpolateC2F(; zero_bcs...) - - @inbounds for i in 1:N_up - - w_up = prog_up_f[i].w - w_en = aux_en_f.w - # Augment the tendencies of updraft area, tracers and vertical velocity - - # entrainment and detrainment - could be moved to implicit - volume_term_entr = @. -p_c / ρ_c * - prog_up[i].ρarea * - Ic(wcomponent(CCG.WVector(w_up))) * - (aux_up[i].entr - aux_up[i].detr) - @. tendencies_up[i].ρarea += - prog_up[i].ρarea * - Ic(wcomponent(CCG.WVector(w_up))) * - (aux_up[i].entr - aux_up[i].detr) - @. tendencies_up[i].ρae_tot += - prog_up[i].ρarea * - aux_en.h_tot * - Ic(wcomponent(CCG.WVector(w_up))) * - aux_up[i].entr - - prog_up[i].ρarea * - aux_up[i].h_tot * - Ic(wcomponent(CCG.WVector(w_up))) * - aux_up[i].detr + volume_term_entr - @. tendencies_up[i].ρaq_tot += - prog_up[i].ρarea * - aux_en.q_tot * - Ic(wcomponent(CCG.WVector(w_up))) * - aux_up[i].entr - - prog_up[i].ρaq_tot * - Ic(wcomponent(CCG.WVector(w_up))) * - aux_up[i].detr - @. tendencies_up_f[i].w += - w_up * I0f(aux_up[i].entr) * (wcomponent(CCG.WVector(w_en - w_up))) - - # precipitation formation - @. tendencies_up[i].ρae_tot += - prog_gm.ρ * aux_up[i].e_tot_tendency_precip_formation - @. tendencies_up[i].ρaq_tot += - prog_gm.ρ * aux_up[i].qt_tendency_precip_formation - - # buoyancy and pressure - @. tendencies_up_f[i].w += CCG.Covariant3Vector( - CCG.WVector(aux_up_f[i].buoy + aux_up_f[i].nh_pressure), - ) - - # TODO - to be removed? - tendencies_up[i].ρarea[kc_surf] = 0 - tendencies_up[i].ρae_tot[kc_surf] = 0 - tendencies_up[i].ρaq_tot[kc_surf] = 0 - tendencies_up_f[i].w[kf_surf] = zero(tendencies_up_f[i].w[kf_surf]) - end - return nothing -end - -function filter_updraft_vars( - edmf::EDMFModel, - grid::Grid, - state::State, - param_set::APS, -) - (; surface_conditions) = state - N_up = n_updrafts(edmf) - kc_surf = kc_surface(grid) - FT = float_type(state) - N_up = n_updrafts(edmf) - - thermo_params = TCP.thermodynamics_params(param_set) - prog_up = center_prog_updrafts(state) - prog_gm = center_prog_grid_mean(state) - prog_gm_uₕ = grid_mean_uₕ(state) - aux_gm_f = face_aux_grid_mean(state) - aux_bulk = center_aux_bulk(state) - prog_up_f = face_prog_updrafts(state) - ρ_c = prog_gm.ρ - ρ_f = aux_gm_f.ρ - ts_gm = center_aux_grid_mean_ts(state) - p_c = TD.air_pressure.(thermo_params, ts_gm) - If = CCO.InterpolateC2F(; - bottom = CCO.Extrapolate(), - top = CCO.Extrapolate(), - ) - @. ρ_f = If(ρ_c) - a_min = edmf.minimum_area - a_max = edmf.max_area - - @inbounds for i in 1:N_up - @. aux_bulk.filter_flag_1 = ifelse(prog_up[i].ρarea < FT(0), 1, 0) - @. aux_bulk.filter_flag_3 = ifelse(prog_up[i].ρaq_tot < FT(0), 1, 0) - @. aux_bulk.filter_flag_4 = ifelse(prog_up[i].ρarea > ρ_c * a_max, 1, 0) - - @. prog_up[i].ρarea = max(prog_up[i].ρarea, 0) #flag_1 - @. prog_up[i].ρaq_tot = max(prog_up[i].ρaq_tot, 0) #flag_3 - @. prog_up[i].ρarea = min(prog_up[i].ρarea, ρ_c * a_max) #flag_4 - end - @inbounds for i in 1:N_up - @. prog_up_f[i].w = CCG.Covariant3Vector( - CCG.WVector(max(wcomponent(CCG.WVector(prog_up_f[i].w)), 0)), - ) - a_up_bcs = a_up_boundary_conditions(surface_conditions, edmf) - If = CCO.InterpolateC2F(; a_up_bcs...) - @. prog_up_f[i].w = - Int(If(prog_up[i].ρarea) >= ρ_f * a_min) * prog_up_f[i].w - end - - Δz = Fields.Δz_field(axes(ρ_c)) - z = Fields.coordinate_field(axes(ρ_c)).z - Δz1 = Spaces.level(Δz, 1) - @inbounds for i in 1:N_up - # this is needed to make sure Rico is unchanged. - # TODO : look into it further to see why - # a similar filtering of ρaθ_liq_ice breaks the simulation - @. prog_up[i].ρaq_tot = ifelse( - z > Δz1 / 2, - ifelse( - prog_up[i].ρarea / ρ_c < a_min, - 0, - max(prog_up[i].ρaq_tot, 0), - ), - prog_up[i].ρaq_tot, - ) - @. prog_up[i].ρarea = ifelse( - z > Δz1 / 2, - ifelse(prog_up[i].ρarea / ρ_c < a_min, 0, max(prog_up[i].ρarea, 0)), - prog_up[i].ρarea, - ) - @. prog_up[i].ρae_tot = ifelse( - z > Δz1 / 2, - ifelse(prog_up[i].ρarea / ρ_c < a_min, 0, prog_up[i].ρae_tot), - prog_up[i].ρae_tot, - ) - end - - Ic = CCO.InterpolateF2C() - C123 = CCG.Covariant123Vector - @inbounds for i in 1:N_up - @. prog_up[i].ρarea = ifelse( - Ic(wcomponent(CCG.WVector(prog_up_f[i].w))) <= 0, - FT(0), - prog_up[i].ρarea, - ) - @. prog_up[i].ρae_tot = ifelse( - Ic(wcomponent(CCG.WVector(prog_up_f[i].w))) <= 0, - FT(0), - prog_up[i].ρae_tot, - ) - @. prog_up[i].ρaq_tot = ifelse( - Ic(wcomponent(CCG.WVector(prog_up_f[i].w))) <= 0, - FT(0), - prog_up[i].ρaq_tot, - ) - - θ_surf = θ_surface_bc(grid, state, edmf, i, param_set) - q_surf = q_surface_bc(grid, state, edmf, i, param_set) - a_surf = area_surface_bc(surface_conditions, edmf) - e_kin = @. LA.norm_sqr( - C123(prog_gm_uₕ) + C123(Ic(CCG.WVector(prog_up_f[i].w))), - ) / 2 - e_pot_surf = geopotential(thermo_params, grid.zc.z[kc_surf]) - ts_up_i_surf = - TD.PhaseEquil_pθq(thermo_params, p_c[kc_surf], θ_surf, q_surf) - e_tot_surf = TD.total_energy( - thermo_params, - ts_up_i_surf, - e_kin[kc_surf], - e_pot_surf, - ) - prog_up[i].ρarea[kc_surf] = ρ_c[kc_surf] * a_surf - prog_up[i].ρae_tot[kc_surf] = prog_up[i].ρarea[kc_surf] * e_tot_surf - prog_up[i].ρaq_tot[kc_surf] = prog_up[i].ρarea[kc_surf] * q_surf - end - return nothing -end diff --git a/src/TurbulenceConvection_deprecated/Fields.jl b/src/TurbulenceConvection_deprecated/Fields.jl deleted file mode 100644 index 9319e5c523e..00000000000 --- a/src/TurbulenceConvection_deprecated/Fields.jl +++ /dev/null @@ -1,146 +0,0 @@ -# A complementary struct to ClimaCore's `PlusHalf` type. -struct Cent{I <: Integer} - i::I -end - -Base.:+(h::Cent) = h -Base.:-(h::Cent) = Cent(-h.i - one(h.i)) -Base.:+(i::Integer, h::Cent) = Cent(i + h.i) -Base.:+(h::Cent, i::Integer) = Cent(h.i + i) -Base.:+(h1::Cent, h2::Cent) = h1.i + h2.i + one(h1.i) -Base.:-(i::Integer, h::Cent) = Cent(i - h.i - one(h.i)) -Base.:-(h::Cent, i::Integer) = Cent(h.i - i) -Base.:-(h1::Cent, h2::Cent) = h1.i - h2.i - -Base.:<=(h1::Cent, h2::Cent) = h1.i <= h2.i -Base.:<(h1::Cent, h2::Cent) = h1.i < h2.i -Base.max(h1::Cent, h2::Cent) = Cent(max(h1.i, h2.i)) -Base.min(h1::Cent, h2::Cent) = Cent(min(h1.i, h2.i)) - -wcomponent(x::CCG.WVector) = x.w - -const FDFields = Union{ - CC.Fields.ExtrudedFiniteDifferenceField, - CC.Fields.FiniteDifferenceField, -} - -const FaceFields = Union{ - CC.Fields.FaceExtrudedFiniteDifferenceField, - CC.Fields.FaceFiniteDifferenceField, -} - -const CenterFields = Union{ - CC.Fields.CenterExtrudedFiniteDifferenceField, - CC.Fields.CenterFiniteDifferenceField, -} - -Base.@propagate_inbounds Base.getindex(field::FDFields, i::Integer) = - Base.getproperty(field, i) - -Base.@propagate_inbounds Base.getindex(field::CenterFields, i::Cent) = - Base.getindex(CC.Fields.field_values(field), i.i) -Base.@propagate_inbounds Base.setindex!(field::CenterFields, v, i::Cent) = - Base.setindex!(CC.Fields.field_values(field), v, i.i) - -Base.@propagate_inbounds Base.getindex(field::FaceFields, i::CCO.PlusHalf) = - Base.getindex(CC.Fields.field_values(field), i.i) -Base.@propagate_inbounds Base.setindex!(field::FaceFields, v, i::CCO.PlusHalf) = - Base.setindex!(CC.Fields.field_values(field), v, i.i) - -Base.@propagate_inbounds Base.getindex(field::FaceFields, ::Cent) = - error("Attempting to getindex with a center index (Cent) into a Face field") -Base.@propagate_inbounds Base.getindex(field::CenterFields, ::CCO.PlusHalf) = - error( - "Attempting to getindex with a face index (PlusHalf) into a Center field", - ) - -Base.@propagate_inbounds Base.setindex!(field::FaceFields, v, ::Cent) = - error("Attempting to setindex with a center index (Cent) into a Face field") -Base.@propagate_inbounds Base.setindex!( - field::CenterFields, - v, - ::CCO.PlusHalf, -) = error( - "Attempting to setindex with a face index (PlusHalf) into a Center field", -) - -function Base.cumsum(field::CC.Fields.FiniteDifferenceField) - Base.cumsum( - parent(CC.Fields.weighted_jacobian(field)) .* vec(field), - dims = 1, - ) -end -function Base.cumsum!( - fieldout::CC.Fields.FiniteDifferenceField, - fieldin::CC.Fields.FiniteDifferenceField, -) - Base.cumsum!( - fieldout, - parent(CC.Fields.weighted_jacobian(fieldin)) .* vec(fieldin), - dims = 1, - ) -end - -# TODO: move these things into ClimaCore - -const CallableZType = Union{Function, Dierckx.Spline1D} - -function set_z!( - field::CC.Fields.Field, - u::CallableZType = x -> x, - v::CallableZType = y -> y, -) - z = CC.Fields.coordinate_field(axes(field)).z - @. field = CCG.Covariant12Vector(CCG.UVVector(u(z), v(z))) -end - -function set_z!(field::CC.Fields.Field, u::Real, v::Real) - lg = CC.Fields.local_geometry_field(axes(field)) - uconst(coord) = u - vconst(coord) = v - @. field = CCG.Covariant12Vector(CCG.UVVector(uconst(lg), vconst(lg))) -end - -# TODO: move to ClimaCore. - -""" - column_findfirstvalue(fn, field::Fields.ColumnField) - -A (;value, z)::NamedTuple containing the value and first z-location -at which the condition, returned by `fn`, is met. If no value -is found, the _last_ value is returned. -""" -function column_findfirstvalue(fn, field::Fields.ColumnField) - space = axes(field) - zfield = Fields.coordinate_field(space).z - li = Operators.left_idx(space) - ri = Operators.right_idx(space) - @inbounds for j in li:ri - level_field = Spaces.level(field, j) - if fn(level_field) - return (; value = level_field, z = Spaces.level(zfield, j)) - end - end - return (; value = Spaces.level(field, ri), z = Spaces.level(zfield, ri)) -end - -""" - column_findlastvalue(fn, field::Fields.ColumnField) - -A (;value, z)::NamedTuple containing the value and last z-location -at which the condition, returned by `fn`, is met. If no value -is found, the _first_ value is returned. -""" -function column_findlastvalue(fn, field::Fields.ColumnField) - space = axes(field) - zfield = Fields.coordinate_field(space).z - li = Operators.left_idx(space) - ri = Operators.right_idx(space) - @inbounds for j in ri:-1:li - level_field = Spaces.level(field, j) - if fn(level_field) - return (; value = level_field, z = Spaces.level(zfield, j)) - end - end - return (; value = Spaces.level(field, li), z = Spaces.level(zfield, li)) -end diff --git a/src/TurbulenceConvection_deprecated/Grid.jl b/src/TurbulenceConvection_deprecated/Grid.jl deleted file mode 100644 index 1fe56658fc3..00000000000 --- a/src/TurbulenceConvection_deprecated/Grid.jl +++ /dev/null @@ -1,41 +0,0 @@ -struct Grid{NZ, CS, FS, SC, SF} - cs::CS - fs::FS - zc::SC - zf::SF - function Grid(space::CC.Spaces.CenterFiniteDifferenceSpace) - nz = length(space) - cs = space - fs = CC.Spaces.FaceFiniteDifferenceSpace(cs) - zc = CC.Fields.coordinate_field(cs) - zf = CC.Fields.coordinate_field(fs) - CS = typeof(cs) - FS = typeof(fs) - SC = typeof(zc) - SF = typeof(zf) - return new{nz, CS, FS, SC, SF}(cs, fs, zc, zf) - end -end - -Grid(mesh::CC.Meshes.IntervalMesh) = - Grid(CC.Spaces.CenterFiniteDifferenceSpace(mesh)) - -function Grid(Δz::FT, nz::Int) where {FT <: AbstractFloat} - z₀, z₁ = FT(0), FT(nz * Δz) - - domain = CC.Domains.IntervalDomain( - CC.Geometry.ZPoint{FT}(z₀), - CC.Geometry.ZPoint{FT}(z₁), - boundary_tags = (:bottom, :top), - ) - - mesh = CC.Meshes.IntervalMesh(domain, nelems = nz) - return Grid(mesh) -end - -n_cells(::Grid{NZ}) where {NZ} = NZ - -# Index of the first interior cell above the surface -kc_surface(grid::Grid) = Cent(1) -kf_surface(grid::Grid) = CCO.PlusHalf(1) -kc_top_of_atmos(grid::Grid) = Cent(n_cells(grid)) diff --git a/src/TurbulenceConvection_deprecated/Operators.jl b/src/TurbulenceConvection_deprecated/Operators.jl deleted file mode 100644 index 4412999379a..00000000000 --- a/src/TurbulenceConvection_deprecated/Operators.jl +++ /dev/null @@ -1,31 +0,0 @@ -##### -##### Masked operators -##### - -""" - shrink_mask!(f, mask) - -Shrinks the subdomains of `f` where `mask == 1`. - -Example: - -```julia -using Test -f = Int[0, 0, 0, 1, 1, 1, 0, 0, 1, 1] -mask = Bool[0, 0, 0, 1, 1, 1, 0, 0, 1, 1] -@test shrink_mask!(f, mask) == Bool[0, 0, 0, 0, 1, 0, 0, 0, 0, 1] -``` -""" -function shrink_mask!(f::AbstractArray, mask::AbstractArray) - @inbounds for (i, m) in enumerate(mask) - if i == 1 || i == length(mask) - f[i] = m - elseif m == 1 && mask[i - 1] == 0 - f[i] = 0 - elseif m == 1 && mask[i + 1] == 0 - f[i] = 0 - else - f[i] = m - end - end -end diff --git a/src/TurbulenceConvection_deprecated/TurbulenceConvection.jl b/src/TurbulenceConvection_deprecated/TurbulenceConvection.jl deleted file mode 100644 index fd177b26ac0..00000000000 --- a/src/TurbulenceConvection_deprecated/TurbulenceConvection.jl +++ /dev/null @@ -1,129 +0,0 @@ -module TurbulenceConvection - -# Import ClimaAtmos types -import ..DryModel -import ..EquilMoistModel -import ..NonEquilMoistModel -import ..RadiationDYCOMS_RF01 -import ..RadiationTRMM_LBA -import ..AbstractPrecipitationModel -import ..Microphysics0Moment -import ..Microphysics1Moment - -import ..compute_kinetic! - -import ClimaCore as CC -import ClimaCore.Geometry as CCG -import ClimaCore.Operators as Operators -import ClimaCore.Fields as Fields -import ClimaCore.Spaces as Spaces -import ClimaCore.Geometry: ⊗, _norm_sqr -import ClimaCore.Operators as CCO -import LinearAlgebra as LA -import DocStringExtensions -import StaticArrays as SA -import StatsBase -import Dierckx -import LambertW -import Thermodynamics as TD -import Distributions -import FastGaussQuadrature -import CloudMicrophysics as CM -import CloudMicrophysics.Microphysics0M as CM0 -import CloudMicrophysics.Microphysics1M as CM1 -import Random - -const liq_type = CM.CommonTypes.LiquidType() -const ice_type = CM.CommonTypes.IceType() -const rain_type = CM.CommonTypes.RainType() -const snow_type = CM.CommonTypes.SnowType() - -include("Parameters.jl") -import .Parameters as TCP -const APS = TCP.AbstractTurbulenceConvectionParameters - -Base.broadcastable(param_set::APS) = tuple(param_set) - -#= - debug_state(state, code_location::String) - -A simple function for debugging the entire state, -specifically for when quantities that should remain -positive-definite become negative. - -=# -function debug_state(state, code_location::String) - prog_gm = center_prog_grid_mean(state) - aux_gm = center_aux_grid_mean(state) - prog_gm_f = face_prog_grid_mean(state) - aux_gm_f = face_aux_grid_mean(state) - - prog_up = center_prog_updrafts(state) - aux_up = center_aux_updrafts(state) - prog_up_f = face_prog_updrafts(state) - aux_up_f = face_aux_updrafts(state) - - prog_en = center_prog_environment(state) - aux_en = center_aux_environment(state) - aux_en_f = face_aux_environment(state) - - ###### - ###### Positive-definite variables - ###### - - vars_positive = [ - vec(prog_gm_f.u₃), - vec(prog_up[1].ρarea), - vec(prog_up_f[1].w), - vec(aux_en.area), - vec(aux_en.θ_liq_ice), - ] - vars = vars_positive - vars_conds = map(v -> any(v .< 0), vars) - - if any(vars_conds) - @show code_location - for (i, vc, v) in zip(1:length(vars), vars_conds, vars) - vc || continue - @show i, v - end - @show vars_conds - error("Negative state for positive-definite field(s)") - end - - ###### - ###### All listed variables - ###### - vars = vars_positive - vars_conds = map(v -> any(isnan.(v)) || any(isinf.(v)), vars) - - if any(vars_conds) - @show code_location - for (i, vc, v) in zip(1:length(vars), vars_conds, vars) - vc || continue - @show i, v - end - @show vars_conds - error("Nan/Inf state for field(s)") - end -end - -include("Grid.jl") -include("dycore_api.jl") -include("Fields.jl") -include("types.jl") -include("Operators.jl") - -include("turbulence_functions.jl") -include("utility_functions.jl") -include("variables.jl") -include("update_aux.jl") -include("EDMF_functions.jl") -include("thermodynamics.jl") -include("closures/microphysics.jl") -include("closures/perturbation_pressure.jl") -include("closures/entr_detr.jl") -include("closures/mixing_length.jl") -include("closures/buoyancy_gradients.jl") - -end diff --git a/src/TurbulenceConvection_deprecated/closures/buoyancy_gradients.jl b/src/TurbulenceConvection_deprecated/closures/buoyancy_gradients.jl deleted file mode 100644 index 42922e39484..00000000000 --- a/src/TurbulenceConvection_deprecated/closures/buoyancy_gradients.jl +++ /dev/null @@ -1,105 +0,0 @@ -""" - buoyancy_gradients( - param_set, - moisture_model, - bg_model::EnvBuoyGrad{FT, EBG} - ) where {FT <: Real, EBG <: AbstractEnvBuoyGradClosure} - -Returns the vertical buoyancy gradients in the environment, as well as in its dry and cloudy volume fractions. -The dispatch on EnvBuoyGrad type is performed at the EnvBuoyGrad construction time, and the analytical solutions -used here are consistent for both mean fields and conditional fields obtained from assumed distributions -over the conserved thermodynamic variables. -""" -function buoyancy_gradients( - param_set::APS, - moisture_model, - bg_model::EnvBuoyGrad{FT, EBG}, -) where {FT <: Real, EBG <: AbstractEnvBuoyGradClosure} - - thermo_params = TCP.thermodynamics_params(param_set) - g = TCP.grav(param_set) - molmass_ratio = TCP.molmass_ratio(param_set) - R_d = TCP.R_d(param_set) - R_v = TCP.R_v(param_set) - - phase_part = TD.PhasePartition(FT(0), FT(0), FT(0)) # assuming R_d = R_m - Π = TD.exner_given_pressure(thermo_params, bg_model.p, phase_part) - - ∂b∂θv = g * (R_d * bg_model.ρ / bg_model.p) * Π - - if bg_model.en_cld_frac > 0.0 - ts_sat = if moisture_model isa DryModel - TD.PhaseDry_pθ(thermo_params, bg_model.p, bg_model.θ_liq_ice_sat) - elseif moisture_model isa EquilMoistModel - TD.PhaseEquil_pθq( - thermo_params, - bg_model.p, - bg_model.θ_liq_ice_sat, - bg_model.qt_sat, - ) - elseif moisture_model isa NonEquilMoistModel - error("Unsupported moisture model") - end - - phase_part = TD.PhasePartition(thermo_params, ts_sat) - lh = TD.latent_heat_liq_ice(thermo_params, phase_part) - cp_m = TD.cp_m(thermo_params, ts_sat) - ∂b∂θl_sat = ( - ∂b∂θv * ( - 1 + - molmass_ratio * - (1 + lh / R_v / bg_model.t_sat) * - bg_model.qv_sat - bg_model.qt_sat - ) / ( - 1 + - lh * lh / cp_m / R_v / bg_model.t_sat / bg_model.t_sat * - bg_model.qv_sat - ) - ) - ∂b∂qt_sat = - (lh / cp_m / bg_model.t_sat * ∂b∂θl_sat - ∂b∂θv) * bg_model.θ_sat - else - ∂b∂θl_sat = FT(0) - ∂b∂qt_sat = FT(0) - end - - ∂b∂z, ∂b∂z_unsat, ∂b∂z_sat = - buoyancy_gradient_chain_rule(bg_model, ∂b∂θv, ∂b∂θl_sat, ∂b∂qt_sat) - return GradBuoy{FT}(∂b∂z, ∂b∂z_unsat, ∂b∂z_sat) -end - -""" - buoyancy_gradient_chain_rule( - bg_model::EnvBuoyGrad{FT, EBG}, - ∂b∂θv::FT, - ∂b∂θl_sat::FT, - ∂b∂qt_sat::FT, - ) where {FT <: Real, EBG <: AbstractEnvBuoyGradClosure} - -Returns the vertical buoyancy gradients in the environment, as well as in its dry and cloudy volume fractions, -from the partial derivatives with respect to thermodynamic variables in dry and cloudy volumes. -""" -function buoyancy_gradient_chain_rule( - bg_model::EnvBuoyGrad{FT, EBG}, - ∂b∂θv::FT, - ∂b∂θl_sat::FT, - ∂b∂qt_sat::FT, -) where {FT <: Real, EBG <: AbstractEnvBuoyGradClosure} - if bg_model.en_cld_frac > FT(0) - ∂b∂z_θl_sat = ∂b∂θl_sat * bg_model.∂θl∂z_sat - ∂b∂z_qt_sat = ∂b∂qt_sat * bg_model.∂qt∂z_sat - else - ∂b∂z_θl_sat = FT(0) - ∂b∂z_qt_sat = FT(0) - end - - ∂b∂z_unsat = - bg_model.en_cld_frac < FT(1) ? ∂b∂θv * bg_model.∂θv∂z_unsat : FT(0) - - ∂b∂z_sat = ∂b∂z_θl_sat + ∂b∂z_qt_sat - ∂b∂z = - (1 - bg_model.en_cld_frac) * ∂b∂z_unsat + - bg_model.en_cld_frac * ∂b∂z_sat - - return ∂b∂z, ∂b∂z_unsat, ∂b∂z_sat -end diff --git a/src/TurbulenceConvection_deprecated/closures/entr_detr.jl b/src/TurbulenceConvection_deprecated/closures/entr_detr.jl deleted file mode 100644 index c0368b1e762..00000000000 --- a/src/TurbulenceConvection_deprecated/closures/entr_detr.jl +++ /dev/null @@ -1,26 +0,0 @@ -""" - function pi_groups_detrainment!(gm_tke, up_area, up_RH, en_area, en_tke, en_RH) - - gm_tke - grid mean turbulent kinetic energy - - up_area - updraft area - - up_RH - updraft relative humidity - - en_area - environment area - - en_tke - environment turbulent kinetic energy - - en_RH - environment relative humidity - - Computes detrainment based on Π-groups -""" -function pi_groups_detrainment!( - gm_tke::FT, - up_area::FT, - up_RH::FT, - en_area::FT, - en_tke::FT, - en_RH::FT, -) where {FT} - - Π2 = (gm_tke - en_area * en_tke) / (gm_tke + eps(FT)) - Π4 = up_RH - en_RH - - return up_area > FT(0) ? max(-0.0102 + 0.0612 * Π2 + 0.0827 * Π4, FT(0)) : - FT(0) -end diff --git a/src/TurbulenceConvection_deprecated/closures/microphysics.jl b/src/TurbulenceConvection_deprecated/closures/microphysics.jl deleted file mode 100644 index 702aed6537f..00000000000 --- a/src/TurbulenceConvection_deprecated/closures/microphysics.jl +++ /dev/null @@ -1,236 +0,0 @@ -""" -Computes the tendencies to qt and θ_liq_ice due to precipitation formation -(autoconversion + accretion) -""" -function precipitation_formation( - param_set::APS, - precip_model::AbstractPrecipitationModel, - prog_gm, - area::FT, - z::FT, - Δt::Real, - ts, -) where {FT} - thermo_params = TCP.thermodynamics_params(param_set) - - microphys_params = TCP.microphysics_params(param_set) - # TODO - when using adaptive timestepping we are limiting the source terms - # with the previous timestep Δt - qt_tendency = FT(0) - ql_tendency = FT(0) - qi_tendency = FT(0) - qr_tendency = FT(0) - qs_tendency = FT(0) - θ_liq_ice_tendency = FT(0) - e_tot_tendency = FT(0) - - if area > 0 - - q = TD.PhasePartition(thermo_params, ts) - - Π_m = TD.exner(thermo_params, ts) - c_pm = TD.cp_m(thermo_params, ts) - L_v0 = TCP.LH_v0(param_set) - L_s0 = TCP.LH_s0(param_set) - I_l = TD.internal_energy_liquid(thermo_params, ts) - I_i = TD.internal_energy_ice(thermo_params, ts) - I = TD.internal_energy(thermo_params, ts) - Φ = geopotential(thermo_params, z) - - if precip_model isa Microphysics0Moment - qsat = TD.q_vap_saturation(thermo_params, ts) - λ = TD.liquid_fraction(thermo_params, ts) - - S_qt = - -min( - (q.liq + q.ice) / Δt, - -CM0.remove_precipitation(microphys_params, q, qsat), - ) - - qr_tendency -= S_qt * λ - qs_tendency -= S_qt * (1 - λ) - qt_tendency += S_qt - ql_tendency += S_qt * λ - qi_tendency += S_qt * (1 - λ) - θ_liq_ice_tendency -= - S_qt / Π_m / c_pm * (L_v0 * λ + L_s0 * (1 - λ)) - e_tot_tendency += - S_qt * - e_tot_0M_precipitation_sources_helper(thermo_params, ts, Φ) - end - - if precip_model isa Microphysics1Moment - T = TD.air_temperature(thermo_params, ts) - T_fr = TCP.T_freeze(param_set) - c_vl = TCP.cv_l(param_set) - c_vm = TD.cv_m(thermo_params, ts) - Rm = TD.gas_constant_air(thermo_params, ts) - Lf = TD.latent_heat_fusion(thermo_params, ts) - - ρ = prog_gm.ρ - qr = max(FT(0), prog_gm.ρq_rai) / ρ - qs = max(FT(0), prog_gm.ρq_sno) / ρ - - # Autoconversion of cloud ice to snow is done with a simplified rate. - # The saturation adjustment scheme prevents using the - # 1-moment snow autoconversion rate that assumes - # that the supersaturation is present in the domain. - S_qt_rain = - -min( - q.liq / Δt, - CM1.conv_q_liq_to_q_rai(microphys_params, q.liq), - ) - S_qt_snow = - -min( - q.ice / Δt, - CM1.conv_q_ice_to_q_sno_no_supersat( - microphys_params, - q.ice, - ), - ) - qr_tendency -= S_qt_rain - qs_tendency -= S_qt_snow - qt_tendency += S_qt_rain + S_qt_snow - ql_tendency += S_qt_rain - qi_tendency += S_qt_snow - θ_liq_ice_tendency -= - 1 / Π_m / c_pm * (L_v0 * S_qt_rain + L_s0 * S_qt_snow) - e_tot_tendency += S_qt_rain * (I_l + Φ) + S_qt_snow * (I_i + Φ) - - # accretion cloud water + rain - S_qr = min( - q.liq / Δt, - CM1.accretion( - microphys_params, - liq_type, - rain_type, - q.liq, - qr, - ρ, - ), - ) - qr_tendency += S_qr - qt_tendency -= S_qr - ql_tendency -= S_qr - θ_liq_ice_tendency += S_qr / Π_m / c_pm * L_v0 - e_tot_tendency -= S_qr * (I_l + Φ) - - # accretion cloud ice + snow - S_qs = min( - q.ice / Δt, - CM1.accretion( - microphys_params, - ice_type, - snow_type, - q.ice, - qs, - ρ, - ), - ) - qs_tendency += S_qs - qt_tendency -= S_qs - qi_tendency -= S_qs - θ_liq_ice_tendency += S_qs / Π_m / c_pm * L_s0 - e_tot_tendency -= S_qs * (I_i + Φ) - - # sink of cloud water via accretion cloud water + snow - S_qt = - -min( - q.liq / Δt, - CM1.accretion( - microphys_params, - liq_type, - snow_type, - q.liq, - qs, - ρ, - ), - ) - if T < T_fr # cloud droplets freeze to become snow) - qs_tendency -= S_qt - qt_tendency += S_qt - ql_tendency += S_qt - θ_liq_ice_tendency -= S_qt / Π_m / c_pm * Lf * (1 + Rm / c_vm) - e_tot_tendency += S_qt * (I_i + Φ) - else # snow melts, both cloud water and snow become rain - α::FT = c_vl / Lf * (T - T_fr) - qt_tendency += S_qt - ql_tendency += S_qt - qs_tendency += S_qt * α - qr_tendency -= S_qt * (1 + α) - θ_liq_ice_tendency += - S_qt / Π_m / c_pm * (Lf * (1 + Rm / c_vm) * α - L_v0) - e_tot_tendency += S_qt * ((1 + α) * I_l - α * I_i + Φ) - end - - # sink of cloud ice via accretion cloud ice - rain - S_qt = - -min( - q.ice / Δt, - CM1.accretion( - microphys_params, - ice_type, - rain_type, - q.ice, - qr, - ρ, - ), - ) - # sink of rain via accretion cloud ice - rain - S_qr = - -min( - qr / Δt, - CM1.accretion_rain_sink(microphys_params, q.ice, qr, ρ), - ) - qt_tendency += S_qt - qi_tendency += S_qt - qr_tendency += S_qr - qs_tendency += -(S_qt + S_qr) - θ_liq_ice_tendency -= - 1 / Π_m / c_pm * (S_qr * Lf * (1 + Rm / c_vm) + S_qt * L_s0) - e_tot_tendency += S_qt * (I_i + Φ) - e_tot_tendency -= S_qr * Lf - - # accretion rain - snow - if T < T_fr - S_qs = min( - qr / Δt, - CM1.accretion_snow_rain( - microphys_params, - snow_type, - rain_type, - qs, - qr, - ρ, - ), - ) - else - S_qs = - -min( - qs / Δt, - CM1.accretion_snow_rain( - microphys_params, - rain_type, - snow_type, - qr, - qs, - ρ, - ), - ) - end - qs_tendency += S_qs - qr_tendency -= S_qs - θ_liq_ice_tendency += S_qs * Lf / Π_m / c_vm - e_tot_tendency += S_qs * Lf - end - end - return PrecipFormation{FT}( - θ_liq_ice_tendency, - e_tot_tendency, - qt_tendency, - ql_tendency, - qi_tendency, - qr_tendency, - qs_tendency, - ) -end diff --git a/src/TurbulenceConvection_deprecated/closures/mixing_length.jl b/src/TurbulenceConvection_deprecated/closures/mixing_length.jl deleted file mode 100644 index 46075ff7f14..00000000000 --- a/src/TurbulenceConvection_deprecated/closures/mixing_length.jl +++ /dev/null @@ -1,85 +0,0 @@ -""" -function mixing_length(param_set, ustar_surf, zc, tke_surf, ∂b∂z, tke, - obukhov_length, Shear², Pr, b_exch) -where: -- `param_set`: set with model parameters -- `ustar_surf`: friction velocity -- `zc`: height -- `tke_surf`: env kinetic energy at first cell center -- `∂b∂z`: buoyancy gradient -- `tke`: env turbulent kinetic energy -- `obukhov_length`: surface Monin Obukhov length -- `Shear²`: shear term -- `Pr`: Prandtl number -- `b_exch`: subdomain echange term - -Returns mixing length as a smooth minimum between -wall-constrained length scale, -production-dissipation balanced length scale and -effective static stability length scale. -""" -function mixing_length( - param_set, - ustar_surf::FT, - zc::FT, - tke_surf::FT, - ∂b∂z::FT, - tke::FT, - obukhov_length::FT, - Shear²::FT, - Pr::FT, - b_exch::FT, -) where {FT} - - c_m = TCP.tke_ed_coeff(param_set) - c_d = TCP.tke_diss_coeff(param_set) - smin_ub = TCP.smin_ub(param_set) - smin_rm = TCP.smin_rm(param_set) - l_max = TCP.l_max(param_set) - c_b = TCP.static_stab_coeff(param_set) - vkc = TCP.von_karman_const(param_set) - - # compute the l_W - the wall constraint mixing length - # which imposes an upper limit on the size of eddies near the surface - # kz scale (surface layer) - if obukhov_length < 0.0 #unstable - l_W = - vkc * zc / (sqrt(tke_surf / ustar_surf / ustar_surf) * c_m) * - min((1 - 100 * zc / obukhov_length)^FT(0.2), 1 / vkc) - else # neutral or stable - l_W = vkc * zc / (sqrt(tke_surf / ustar_surf / ustar_surf) * c_m) - end - - # compute l_TKE - the production-dissipation balanced length scale - a_pd = c_m * (Shear² - ∂b∂z / Pr) * sqrt(tke) - # Dissipation term - c_neg = c_d * tke * sqrt(tke) - if abs(a_pd) > eps(FT) && 4 * a_pd * c_neg > -b_exch * b_exch - l_TKE = max( - -b_exch / 2 / a_pd + - sqrt(b_exch * b_exch + 4 * a_pd * c_neg) / 2 / a_pd, - 0, - ) - elseif abs(a_pd) < eps(FT) && abs(b_exch) > eps(FT) - l_TKE = c_neg / b_exch - else - l_TKE = FT(0) - end - - # compute l_N - the effective static stability length scale. - N_eff = sqrt(max(∂b∂z, 0)) - if N_eff > 0.0 - l_N = min(sqrt(max(c_b * tke, 0)) / N_eff, l_max) - else - l_N = l_max - end - - # add limiters - l = SA.SVector( - (l_N < eps(FT) || l_N > l_max) ? l_max : l_N, - (l_TKE < eps(FT) || l_TKE > l_max) ? l_max : l_TKE, - (l_W < eps(FT) || l_W > l_max) ? l_max : l_W, - ) - # get soft minimum - return lamb_smooth_minimum(l, smin_ub, smin_rm) -end diff --git a/src/TurbulenceConvection_deprecated/closures/perturbation_pressure.jl b/src/TurbulenceConvection_deprecated/closures/perturbation_pressure.jl deleted file mode 100644 index 6aacf0553b0..00000000000 --- a/src/TurbulenceConvection_deprecated/closures/perturbation_pressure.jl +++ /dev/null @@ -1,48 +0,0 @@ -function compute_updraft_top(updraft_area) - (; value, z) = column_findlastvalue(a -> a[] > 1e-3, updraft_area) - return z[] -end - -""" - compute_nh_pressure!( - param_set, - buoy_up, - w_up, - div_w_up, - w_up_en_diff, - updraft_top, - ) -where - - `param_set`: contains model parameters - - `buoy_up`: updraft buoyancy - - `w_up`: updraft vertical velocity - - `div_w_up`: updraft vertical velocity divergence - - `w_up_en_diff`: difference between updraft and environment vertical velocity - - `updraft_top`: height at which updraft terminates - -Returns the updraft perturbation pressure gradient advection and drag -following [He2020](@cite) -""" -function compute_nh_pressure!( - param_set, - buoy_up, - w_up, - div_w_up, - w_up_en_diff, - updraft_top, -) - # factor multiplier for pressure buoyancy terms (effective buoyancy is (1-α_b)) - α_b = TCP.pressure_normalmode_buoy_coeff1(param_set) - # factor multiplier for pressure advection - α_a = TCP.pressure_normalmode_adv_coeff(param_set) - # factor multiplier for pressure drag - α_d = TCP.pressure_normalmode_drag_coeff(param_set) - - # Independence of aspect ratio hardcoded: α₂_asp_ratio² = FT(0) - - H_up_min = TCP.min_updraft_top(param_set) - plume_scale_height = max(updraft_top, H_up_min) - - return -α_b * buoy_up + α_a * w_up * div_w_up - - α_d * w_up_en_diff * abs(w_up_en_diff) / plume_scale_height -end diff --git a/src/TurbulenceConvection_deprecated/dycore_api.jl b/src/TurbulenceConvection_deprecated/dycore_api.jl deleted file mode 100644 index 23d7ad1ff16..00000000000 --- a/src/TurbulenceConvection_deprecated/dycore_api.jl +++ /dev/null @@ -1,82 +0,0 @@ -##### -##### Dycore API -##### - -abstract type FieldLocation end -struct CentField <: FieldLocation end -struct FaceField <: FieldLocation end -struct SingleValuePerColumn <: FieldLocation end - -field_loc(::CentField) = :c -field_loc(::FaceField) = :f -field_loc(::SingleValuePerColumn) = :svpc - -#= -This file provides a list of methods that TurbulenceConvection.jl -expects that the host dycore supports. This is experimental, as -we're not sure how the data structures / flow control will shake out. -=# - -##### -##### Grid mean fields -##### - -""" Prognostic fields for the host model """ -prognostic(state, fl) = getproperty(state.prog, field_loc(fl)) - -center_prog_grid_mean(state) = prognostic(state, CentField()) -face_prog_grid_mean(state) = prognostic(state, FaceField()) - -""" Auxiliary fields for the host model """ -aux(state, fl) = getproperty(state.aux, field_loc(fl)) - -center_aux_grid_mean_ts(state) = state.p.ᶜts[state.colidx] -center_aux_grid_mean_p(state) = state.p.ᶜp[state.colidx] -center_aux_grid_mean_e_kin(state) = state.p.ᶜK[state.colidx] -center_aux_grid_mean(state) = aux(state, CentField()) -face_aux_grid_mean(state) = aux(state, FaceField()) - -""" Tendency fields for the host model """ -tendencies(state, fl) = getproperty(state.tendencies, field_loc(fl)) - -center_tendencies_grid_mean(state) = tendencies(state, CentField()) - -##### -##### TurbulenceConvection fields -##### - -#= Prognostic fields for TurbulenceConvection =# -prognostic_turbconv(state, fl) = prognostic(state, fl).turbconv - -center_prog_updrafts(state) = prognostic_turbconv(state, CentField()).up -face_prog_updrafts(state) = prognostic_turbconv(state, FaceField()).up -center_prog_environment(state) = prognostic_turbconv(state, CentField()).en - -#= Auxiliary fields for TurbulenceConvection =# -aux_turbconv(state, fl) = aux(state, fl).turbconv - -center_aux_turbconv(state) = aux_turbconv(state, CentField()) -face_aux_turbconv(state) = aux_turbconv(state, FaceField()) -center_aux_updrafts(state) = aux_turbconv(state, CentField()).up -face_aux_updrafts(state) = aux_turbconv(state, FaceField()).up -center_aux_environment(state) = aux_turbconv(state, CentField()).en -center_aux_bulk(state) = aux_turbconv(state, CentField()).bulk -face_aux_bulk(state) = aux_turbconv(state, FaceField()).bulk -face_aux_environment(state) = aux_turbconv(state, FaceField()).en - -#= Tendency fields for TurbulenceConvection =# -tendencies_turbconv(state, fl) = tendencies(state, fl).turbconv - -center_tendencies_turbconv(state) = tendencies_turbconv(state, CentField()) -face_tendencies_turbconv(state) = tendencies_turbconv(state, FaceField()) -center_tendencies_updrafts(state) = tendencies_turbconv(state, CentField()).up -center_tendencies_environment(state) = - tendencies_turbconv(state, CentField()).en -face_tendencies_updrafts(state) = tendencies_turbconv(state, FaceField()).up - -physical_grid_mean_uₕ(state) = CCG.UVVector.(grid_mean_uₕ(state)) -physical_grid_mean_u(state) = map(x -> x.u, physical_grid_mean_uₕ(state)) -physical_grid_mean_v(state) = map(x -> x.v, physical_grid_mean_uₕ(state)) -grid_mean_uₕ(state) = center_prog_grid_mean(state).uₕ - -tendencies_grid_mean_uₕ(state) = center_tendencies_grid_mean(state).uₕ diff --git a/src/TurbulenceConvection_deprecated/thermodynamics.jl b/src/TurbulenceConvection_deprecated/thermodynamics.jl deleted file mode 100644 index 70fc394994d..00000000000 --- a/src/TurbulenceConvection_deprecated/thermodynamics.jl +++ /dev/null @@ -1,10 +0,0 @@ -geopotential(thermo_params::TD.Parameters.ThermodynamicsParameters, z::Real) = - z * TD.Parameters.grav(thermo_params) - -function enthalpy(h_tot::FT, e_kin::FT, e_pot::FT) where {FT} - return h_tot - e_kin - e_pot -end - -function enthalpy(mse::FT, e_pot::FT) where {FT} - return mse - e_pot -end diff --git a/src/TurbulenceConvection_deprecated/turbulence_functions.jl b/src/TurbulenceConvection_deprecated/turbulence_functions.jl deleted file mode 100644 index 80a4efd902b..00000000000 --- a/src/TurbulenceConvection_deprecated/turbulence_functions.jl +++ /dev/null @@ -1,61 +0,0 @@ -function buoyancy_c(thermo_params, ρ::FT, ρ_i::FT) where {FT} - g::FT = TD.Parameters.grav(thermo_params) - return g * (ρ - ρ_i) / ρ -end - -# MO scaling of near surface tke and scalar variance -function get_surface_tke(param_set, ustar::FT, zLL::FT, oblength::FT) where {FT} - κ_star² = TCP.κ_star²(param_set) - if oblength < 0 - return ( - (κ_star² + cbrt(zLL / oblength * zLL / oblength)) * ustar * ustar - ) - else - return κ_star² * ustar * ustar - end -end -function get_surface_variance( - flux, - ustar::FT, - zLL::FT, - oblength::FT, - local_geometry_data, -) where {FT} - c_star = -flux / ustar - if oblength < 0 - return 4 * - _norm_sqr(c_star, local_geometry_data) * - (1 - FT(8.3) * zLL / oblength)^(-FT(2 / 3)) - else - return 4 * _norm_sqr(c_star, local_geometry_data) - end -end - -function gradient_Richardson_number( - param_set, - ∂b∂z::FT, - Shear²::FT, - ϵ::FT, -) where {FT} - Ri_c = TCP.Ri_crit(param_set) - return min(∂b∂z / max(Shear², ϵ), Ri_c) -end - -# Turbulent Prandtl number given the obukhov length sign and the gradient Richardson number -function turbulent_Prandtl_number( - param_set, - obukhov_length::FT, - ∇Ri::FT, -) where {FT} - ω_pr = TCP.Prandtl_number_scale(param_set) - Pr_n = TCP.Prandtl_number_0(param_set) - if obukhov_length > 0 && ∇Ri > 0 #stable - # CSB (Dan Li, 2019, eq. 75), where ω_pr = ω_1 + 1 = 53.0/13.0 - prandtl_nvec = - Pr_n * - (2 * ∇Ri / (1 + ω_pr * ∇Ri - sqrt((1 + ω_pr * ∇Ri)^2 - 4 * ∇Ri))) - else - prandtl_nvec = Pr_n - end - return prandtl_nvec -end diff --git a/src/TurbulenceConvection_deprecated/types.jl b/src/TurbulenceConvection_deprecated/types.jl deleted file mode 100644 index 3534cdb9142..00000000000 --- a/src/TurbulenceConvection_deprecated/types.jl +++ /dev/null @@ -1,255 +0,0 @@ -""" - PrecipFormation - -Storage for tendencies due to precipitation formation - -$(DocStringExtensions.FIELDS) -""" -Base.@kwdef struct PrecipFormation{FT} - θ_liq_ice_tendency::FT = FT(0) - e_tot_tendency::FT = FT(0) - qt_tendency::FT = FT(0) - ql_tendency::FT = FT(0) - qi_tendency::FT = FT(0) - qr_tendency::FT = FT(0) - qs_tendency::FT = FT(0) -end - -""" - GradBuoy - -Environmental buoyancy gradients. - -$(DocStringExtensions.FIELDS) -""" -Base.@kwdef struct GradBuoy{FT} - "environmental vertical buoyancy gradient" - ∂b∂z::FT = FT(0) - "vertical buoyancy gradient in the unsaturated part of the environment" - ∂b∂z_unsat::FT = FT(0) - "vertical buoyancy gradient in the saturated part of the environment" - ∂b∂z_sat::FT = FT(0) -end - -abstract type AbstractEnvBuoyGradClosure end -struct BuoyGradMean <: AbstractEnvBuoyGradClosure end -struct BuoyGradQuadratures <: AbstractEnvBuoyGradClosure end - -Base.broadcastable(x::BuoyGradMean) = tuple(x) -Base.broadcastable(x::BuoyGradQuadratures) = tuple(x) - -""" - EnvBuoyGrad - -Variables used in the environmental buoyancy gradient computation. - -$(DocStringExtensions.FIELDS) -""" -Base.@kwdef struct EnvBuoyGrad{FT, EBC <: AbstractEnvBuoyGradClosure} - "temperature in the saturated part" - t_sat::FT - "vapor specific humidity in the saturated part" - qv_sat::FT - "total specific humidity in the saturated part" - qt_sat::FT - "potential temperature in the saturated part" - θ_sat::FT - "liquid ice potential temperature in the saturated part" - θ_liq_ice_sat::FT - "virtual potential temperature gradient in the non saturated part" - ∂θv∂z_unsat::FT - "total specific humidity gradient in the saturated part" - ∂qt∂z_sat::FT - "liquid ice potential temperature gradient in the saturated part" - ∂θl∂z_sat::FT - "reference pressure" - p::FT - "cloud fraction" - en_cld_frac::FT - "density" - ρ::FT -end -function EnvBuoyGrad( - ::EBG, - t_sat::FT, - args..., -) where {FT <: Real, EBG <: AbstractEnvBuoyGradClosure} - return EnvBuoyGrad{FT, EBG}(t_sat, args...) -end - -""" - MinDisspLen - -Minimum dissipation model - -$(DocStringExtensions.FIELDS) -""" -Base.@kwdef struct MinDisspLen{FT} - "height" - z::FT - "obukhov length" - obukhov_length::FT - "surface TKE values" - tke_surf::FT - "u star - surface velocity scale" - ustar::FT - "turbulent Prandtl number" - Pr::FT - "reference pressure" - p::FT - "vertical buoyancy gradient struct" - ∇b::GradBuoy{FT} - "env shear" - Shear²::FT - "environment turbulent kinetic energy" - tke::FT - "Updraft tke source" - b_exch::FT -end - -const_or_func(s::Function, t::Real) = s(t) -const_or_func(s::Dierckx.Spline1D, t::Real) = s(t) -const_or_func(s, t::Real) = s - -struct EDMFModel{N_up, FT, MM, PM, EBGC} - surface_area::FT - max_area::FT - minimum_area::FT - moisture_model::MM - precip_model::PM - bg_closure::EBGC - zero_uv_fluxes::Bool -end -function EDMFModel( - ::Type{FT}, - moisture_model, - precip_model, - parsed_args, - turbconv_params, -) where {FT} - - tc_case = parsed_args["turbconv_case"] - zero_uv_fluxes = any(tcc -> tcc == tc_case, ["TRMM_LBA", "ARM_SGP"]) - # Set the number of updrafts (1) - n_updrafts = turbconv_params.updraft_number - - surface_area = turbconv_params.surface_area - max_area = turbconv_params.max_area - minimum_area = turbconv_params.min_area - - bg_closure = BuoyGradMean() - if !(moisture_model isa EquilMoistModel) - @warn string( - "SGS model only supports equilibrium moisture choice.", - "EDMF will carry zeros, but avoid saturation adjustment calls.", - ) - end - - MM = typeof(moisture_model) - PM = typeof(precip_model) - EBGC = typeof(bg_closure) - return EDMFModel{n_updrafts, FT, MM, PM, EBGC}( - surface_area, - max_area, - minimum_area, - moisture_model, - precip_model, - bg_closure, - zero_uv_fluxes, - ) -end - -parameter_set(obj) = obj.param_set -n_updrafts(::EDMFModel{N_up}) where {N_up} = N_up -Base.eltype(::EDMFModel{N_up, FT}) where {N_up, FT} = FT - -Base.broadcastable(edmf::EDMFModel) = tuple(edmf) - -function Base.summary(io::IO, edmf::EDMFModel) - pns = string.(propertynames(edmf)) - buf = maximum(length.(pns)) - keys = propertynames(edmf) - vals = repeat.(" ", map(s -> buf - length(s) + 2, pns)) - bufs = (; zip(keys, vals)...) - print(io, '\n') - for pn in propertynames(edmf) - prop = getproperty(edmf, pn) - # Skip some data: - prop isa Bool && continue - prop isa NTuple && continue - prop isa Int && continue - prop isa Float64 && continue - prop isa Float32 && continue - s = string( - " ", # needed for some reason - getproperty(bufs, pn), - '`', - string(pn), - '`', - "::", - '`', - typeof(prop), - '`', - '\n', - ) - print(io, s) - end -end - - -struct State{P, A, T, CACHE, C, SC} - prog::P - aux::A - tendencies::T - p::CACHE - colidx::C - surface_conditions::SC -end - -Grid(state::State) = Grid(axes(state.prog.c)) - -float_type(state::State) = eltype(state.prog) -# float_type(field::CC.Fields.Field) = CC.Spaces.undertype(axes(field)) -float_type(field::CC.Fields.Field) = eltype(parent(field)) - -import ClimaCore.Fields as Fields -import ClimaCore.Spaces as Spaces - - -Base.@propagate_inbounds function field_vector_column( - fv::Fields.FieldVector{T}, - colidx::Fields.ColumnIndex, -) where {T} - values = map(x -> x[colidx], Fields._values(fv)) - return Fields.FieldVector{T, typeof(values)}(values) -end - -function tc_column_state(prog, p, tendencies, colidx, t) - prog_column = field_vector_column(prog, colidx) - aux_column = field_vector_column(p.edmf_cache.aux, colidx) - tends_column = field_vector_column(tendencies, colidx) - surface_conditions = CC.column(p.sfc_conditions, colidx)[] - return State( - prog_column, - aux_column, - tends_column, - p, - colidx, - surface_conditions, - ) -end - -function tc_column_state(prog, p, tendencies::Nothing, colidx, t) - prog_column = field_vector_column(prog, colidx) - aux_column = field_vector_column(p.edmf_cache.aux, colidx) - tends_column = nothing - surface_conditions = CC.column(p.sfc_conditions, colidx)[] - return State( - prog_column, - aux_column, - tends_column, - p, - colidx, - surface_conditions, - ) -end diff --git a/src/TurbulenceConvection_deprecated/update_aux.jl b/src/TurbulenceConvection_deprecated/update_aux.jl deleted file mode 100644 index 91ecab92767..00000000000 --- a/src/TurbulenceConvection_deprecated/update_aux.jl +++ /dev/null @@ -1,557 +0,0 @@ -function update_aux!( - edmf::EDMFModel, - grid::Grid, - state::State, - param_set::APS, - t::Real, - Δt::Real, -) - ##### - ##### Unpack common variables - ##### - surface_conditions = state.surface_conditions - a_bulk_bcs = (; - bottom = CCO.SetValue( - sum( - i -> area_surface_bc(surface_conditions, edmf), - 1:n_updrafts(edmf), - ), - ), - top = CCO.Extrapolate(), - ) - Ifb = CCO.InterpolateC2F(; a_bulk_bcs...) - thermo_params = TCP.thermodynamics_params(param_set) - microphys_params = TCP.microphysics_params(param_set) - N_up = n_updrafts(edmf) - kc_surf = kc_surface(grid) - kf_surf = kf_surface(grid) - c_m = TCP.tke_ed_coeff(param_set) - KM = center_aux_turbconv(state).KM - KH = center_aux_turbconv(state).KH - oblength = surface_conditions.obukhov_length - FT = float_type(state) - prog_gm = center_prog_grid_mean(state) - prog_gm_f = face_prog_grid_mean(state) - aux_up = center_aux_updrafts(state) - aux_up_f = face_aux_updrafts(state) - aux_en = center_aux_environment(state) - aux_en_f = face_aux_environment(state) - aux_gm = center_aux_grid_mean(state) - aux_gm_f = face_aux_grid_mean(state) - aux_tc_f = face_aux_turbconv(state) - aux_tc = center_aux_turbconv(state) - aux_bulk = center_aux_bulk(state) - prog_en = center_prog_environment(state) - prog_up = center_prog_updrafts(state) - prog_up_f = face_prog_updrafts(state) - ρ_f = aux_gm_f.ρ - p_c = center_aux_grid_mean_p(state) - ρ_c = prog_gm.ρ - zc = Fields.coordinate_field(axes(ρ_c)).z - aux_en_unsat = aux_en.unsat - aux_en_sat = aux_en.sat - mph = center_aux_turbconv(state).mph - m_entr_detr = aux_tc.ϕ_temporary - ∇m_entr_detr = aux_tc.ψ_temporary - wvec = CC.Geometry.WVector - max_area = edmf.max_area - ts_gm = center_aux_grid_mean_ts(state) - ts_env = center_aux_environment(state).ts - e_kin = center_aux_grid_mean_e_kin(state) - - prog_gm_uₕ = grid_mean_uₕ(state) - Ic = CCO.InterpolateF2C() - ##### - ##### center variables - ##### - C123 = CCG.Covariant123Vector - - @inbounds for i in 1:N_up - @. aux_up[i].e_kin = - LA.norm_sqr( - C123(prog_gm_uₕ) + C123(Ic(CCG.WVector(prog_up_f[i].w))), - ) / 2 - end - - e_pot(z) = geopotential(thermo_params, z) - thresh_area(prog_up, ρ_c) = prog_up[i].ρarea / ρ_c[k] >= edmf.minimum_area - @inbounds for i in 1:N_up - @. aux_up[i].e_tot = ifelse( - prog_up[i].ρarea / ρ_c >= edmf.minimum_area, - prog_up[i].ρae_tot / prog_up[i].ρarea, - aux_gm.e_tot, - ) - @. aux_up[i].q_tot = ifelse( - prog_up[i].ρarea / ρ_c >= edmf.minimum_area, - prog_up[i].ρaq_tot / prog_up[i].ρarea, - aux_gm.q_tot, - ) - @. aux_up[i].area = ifelse( - prog_up[i].ρarea / ρ_c >= edmf.minimum_area, - prog_up[i].ρarea / ρ_c, - 0, - ) - @. aux_up[i].e_kin = ifelse( - prog_up[i].ρarea / ρ_c >= edmf.minimum_area, - aux_up[i].e_kin, - e_kin, - ) - ##### - ##### Set primitive variables - ##### - e_int = @. aux_up[i].e_tot - aux_up[i].e_kin - e_pot(zc) - if edmf.moisture_model isa DryModel - @. aux_up[i].ts = TD.PhaseDry_pe(thermo_params, p_c, e_int) - elseif edmf.moisture_model isa EquilMoistModel - @. aux_up[i].ts = - TD.PhaseEquil_peq(thermo_params, p_c, e_int, aux_up[i].q_tot) - elseif edmf.moisture_model isa NonEquilMoistModel - error("Unsupported moisture model") - end - - ts_up = aux_up[i].ts - @. aux_up[i].θ_liq_ice = TD.liquid_ice_pottemp(thermo_params, ts_up) - @. aux_up[i].h_tot = - TD.total_specific_enthalpy(thermo_params, ts_up, aux_up[i].e_tot) - @. aux_up[i].q_liq = TD.liquid_specific_humidity(thermo_params, ts_up) - @. aux_up[i].q_ice = TD.ice_specific_humidity(thermo_params, ts_up) - @. aux_up[i].T = TD.air_temperature(thermo_params, ts_up) - @. aux_up[i].RH = TD.relative_humidity(thermo_params, ts_up) - - end - ##### - ##### compute bulk - ##### - @. aux_bulk.area = 0 - @inbounds for i in 1:N_up - @. aux_bulk.area += aux_up[i].area - end - - @. aux_bulk.q_tot = 0 - @. aux_bulk.h_tot = 0 - @inbounds for i in 1:N_up - @. aux_bulk.q_tot += aux_up[i].area * aux_up[i].q_tot / aux_bulk.area - @. aux_bulk.h_tot += aux_up[i].area * aux_up[i].h_tot / aux_bulk.area - end - @inbounds for i in 1:N_up - @. aux_bulk.q_tot = - ifelse(aux_bulk.area > 0, aux_bulk.q_tot, aux_gm.q_tot) - @. aux_bulk.h_tot = - ifelse(aux_bulk.area > 0, aux_bulk.h_tot, aux_gm.h_tot) - end - - @. aux_en.area = 1 - aux_bulk.area - @. aux_en.tke = prog_en.ρatke / (ρ_c * aux_en.area) - - @. aux_en_f.w = prog_gm_f.u₃ / (1 - Ifb(aux_bulk.area)) - @inbounds for i in 1:N_up - @. aux_en_f.w -= - Ifb(aux_up[i].area) * prog_up_f[i].w / (1 - Ifb(aux_bulk.area)) - end - - @. aux_en.e_kin = - LA.norm_sqr(C123(prog_gm_uₕ) + C123(Ic(wvec(aux_en_f.w)))) / 2 - - ##### - ##### decompose_environment - ##### - val1(aux_bulk) = 1 / (1 - aux_bulk.area) - val2(aux_bulk) = aux_bulk.area * val1(aux_bulk) - #Yair - this is here to prevent negative QT - @. aux_en.q_tot = - max(val1(aux_bulk) * aux_gm.q_tot - val2(aux_bulk) * aux_bulk.q_tot, 0) - @. aux_en.h_tot = - val1(aux_bulk) * aux_gm.h_tot - val2(aux_bulk) * aux_bulk.h_tot - - if edmf.moisture_model isa DryModel - @. aux_en.ts = TD.PhaseDry_ph( - thermo_params, - p_c, - enthalpy(aux_en.h_tot, e_pot(zc), aux_en.e_kin), - ) - elseif edmf.moisture_model isa EquilMoistModel - @. aux_en.ts = TD.PhaseEquil_phq( - thermo_params, - p_c, - enthalpy(aux_en.h_tot, e_pot(zc), aux_en.e_kin), - aux_en.q_tot, - ) - elseif edmf.moisture_model isa NonEquilMoistModel - error("Add support got non-equilibrium thermo states") - end - ts_en = aux_en.ts - @. aux_en.θ_liq_ice = TD.liquid_ice_pottemp(thermo_params, ts_en) - @. aux_en.e_tot = - TD.total_energy(thermo_params, ts_en, aux_en.e_kin, e_pot(zc)) - @. aux_en.T = TD.air_temperature(thermo_params, ts_en) - @. aux_en.θ_virt = TD.virtual_pottemp(thermo_params, ts_en) - @. aux_en.θ_dry = TD.dry_pottemp(thermo_params, ts_en) - @. aux_en.q_liq = TD.liquid_specific_humidity(thermo_params, ts_en) - @. aux_en.q_ice = TD.ice_specific_humidity(thermo_params, ts_en) - @. aux_en.RH = TD.relative_humidity(thermo_params, ts_en) - - # compute rain formation in the environment (autoconversion and accretion) - @. mph = precipitation_formation( - param_set, - edmf.precip_model, - prog_gm, - aux_en.area, - zc, - Δt, - ts_en, - ) - # TODO: move ..._tendency_precip_formation to diagnostics - @. aux_en.e_tot_tendency_precip_formation = mph.e_tot_tendency * aux_en.area - @. aux_en.qt_tendency_precip_formation = mph.qt_tendency * aux_en.area - @. aux_en.qr_tendency_precip_formation = mph.qr_tendency * aux_en.area - @. aux_en.qs_tendency_precip_formation = mph.qs_tendency * aux_en.area - - @. aux_en.cloud_fraction = - ifelse(TD.has_condensate(thermo_params, ts_en), 1, 0) - - # TODO: shouldn't we always populate aux_en_sat and aux_en_unsat? - # Otherwise we may be using values from other timesteps / stages - # update_sat_unsat - conditional_assign!(aux_en_sat.θ_dry, ts_en, TD.dry_pottemp, thermo_params) - conditional_assign!( - aux_en_sat.θ_liq_ice, - ts_en, - TD.liquid_ice_pottemp, - thermo_params, - ) - conditional_assign!(aux_en_sat.T, ts_en, TD.air_temperature, thermo_params) - conditional_assign!( - aux_en_sat.q_tot, - ts_en, - TD.total_specific_humidity, - thermo_params, - ) - conditional_assign!( - aux_en_sat.q_vap, - ts_en, - TD.vapor_specific_humidity, - thermo_params, - ) - - conditional_assign!( - aux_en_unsat.θ_dry, - ts_en, - TD.dry_pottemp, - thermo_params, - x -> !x, - ) - conditional_assign!( - aux_en_unsat.θ_virt, - ts_en, - TD.virtual_pottemp, - thermo_params, - x -> !x, - ) - conditional_assign!( - aux_en_unsat.q_tot, - ts_en, - TD.total_specific_humidity, - thermo_params, - x -> !x, - ) - - # compute the buoyancy - LBF_ρ = CCO.LeftBiasedC2F(; bottom = CCO.SetValue(ρ_f[kf_surf])) - #LBF_ρ = CCO.LeftBiasedC2F(; bottom = CCO.SetValue(ρ_c[kc_surf])) - LBF_a = CCO.LeftBiasedC2F(; bottom = CCO.SetValue(edmf.surface_area)) - @. aux_en_f.buoy = buoyancy_c( - thermo_params, - ρ_f, - LBF_ρ(TD.air_density(thermo_params, aux_en.ts)), - ) - @inbounds for i in 1:N_up - @. aux_up_f[i].buoy = buoyancy_c( - thermo_params, - ρ_f, - LBF_ρ(TD.air_density(thermo_params, aux_up[i].ts)), - ) - end - @. aux_gm_f.buoy = (1.0 - LBF_a(aux_bulk.area)) * aux_en_f.buoy - @inbounds for i in 1:N_up - @. aux_gm_f.buoy += LBF_a(aux_up[i].area) * aux_up_f[i].buoy - end - @inbounds for i in 1:N_up - @. aux_up_f[i].buoy -= aux_gm_f.buoy - end - # Only needed for diagnostics/plotting - @. aux_en_f.buoy -= aux_gm_f.buoy - @. aux_tc_f.bulk.buoy_up1 = aux_up_f[1].buoy - - - ##### - ##### compute bulk thermodynamics - ##### - @. aux_bulk.q_liq = 0 - @. aux_bulk.q_ice = 0 - @. aux_bulk.T = 0 - @inbounds for i in 1:N_up - @. aux_bulk.q_liq += aux_up[i].area * aux_up[i].q_liq / aux_bulk.area - @. aux_bulk.q_ice += aux_up[i].area * aux_up[i].q_ice / aux_bulk.area - @. aux_bulk.T += aux_up[i].area * aux_up[i].T / aux_bulk.area - end - @. aux_bulk.q_liq = ifelse(aux_bulk.area > 0, aux_bulk.q_liq, 0) - @. aux_bulk.q_ice = ifelse(aux_bulk.area > 0, aux_bulk.q_ice, 0) - @. aux_bulk.T = ifelse(aux_bulk.area > 0, aux_bulk.T, aux_en.T) - - ##### - ##### update_GMV_diagnostics - ##### - @. aux_gm.q_liq = - aux_bulk.area * aux_bulk.q_liq + (1 - aux_bulk.area) * aux_en.q_liq - @. aux_gm.q_ice = - aux_bulk.area * aux_bulk.q_ice + (1 - aux_bulk.area) * aux_en.q_ice - @. aux_gm.T = aux_bulk.area * aux_bulk.T + (1 - aux_bulk.area) * aux_en.T - - @. aux_bulk.cloud_fraction = ifelse( - TD.has_condensate(aux_bulk.q_liq + aux_bulk.q_ice) && - aux_bulk.area > 1e-3, - 1, - 0, - ) - - @. aux_gm.tke = aux_en.area * aux_en.tke - @. aux_gm.tke += - 0.5 * - aux_en.area * - LA.norm_sqr(C123(Ic(wvec(aux_en_f.w - prog_gm_f.u₃)))) - @inbounds for i in 1:N_up - @. aux_gm.tke += - 0.5 * - aux_up[i].area * - LA.norm_sqr(C123(Ic(wvec(prog_up_f[i].w - prog_gm_f.u₃)))) - end - - ##### - ##### face variables: diagnose primitive, diagnose env and compute bulk - ##### - parent(aux_tc_f.bulk.w) .= 0 - @inbounds for i in 1:N_up - a_up = aux_up[i].area - a_up_bcs = a_up_boundary_conditions(surface_conditions, edmf) - Ifu = CCO.InterpolateC2F(; a_up_bcs...) - @. aux_tc_f.bulk.w += ifelse( - Ifb(aux_bulk.area) > 0, - Ifu(a_up) * prog_up_f[i].w / Ifb(aux_bulk.area), - CCG.Covariant3Vector(FT(0)), - ) - end - - ##### - ##### compute_updraft_closures - ##### - # entrainment and detrainment - @inbounds for i in 1:N_up - @. aux_up[i].entr = FT(5e-4) - @. aux_up[i].detr = pi_groups_detrainment!( - aux_gm.tke, - aux_up[i].area, - aux_up[i].RH, - aux_en.area, - aux_en.tke, - aux_en.RH, - ) - end - # updraft pressure - w_bcs = - (; bottom = CCO.SetValue(wvec(FT(0))), top = CCO.SetValue(wvec(FT(0)))) - ∇ = CCO.DivergenceC2F(; w_bcs...) - @inbounds for i in 1:N_up - updraft_top = compute_updraft_top(aux_up[i].area) - @. aux_up_f[i].nh_pressure = compute_nh_pressure!( - param_set, - aux_up_f[i].buoy, - wcomponent(CCG.WVector(prog_up_f[i].w)), - ∇(Ic(CCG.WVector(prog_up_f[i].w))), - wcomponent(CCG.WVector(prog_up_f[i].w - aux_en_f.w)), - updraft_top, - ) - end - - ##### - ##### compute_eddy_diffusivities_tke - ##### - - # Subdomain exchange term - Ic = CCO.InterpolateF2C() - b_exch = center_aux_turbconv(state).b_exch - parent(b_exch) .= 0 - w_en = aux_en_f.w - @inbounds for i in 1:N_up - w_up = prog_up_f[i].w - @. b_exch += - aux_up[i].area * - Ic(wcomponent(CCG.WVector(w_up))) * - aux_up[i].detr / aux_en.area * - (1 / 2 * (Ic(wcomponent(CCG.WVector(w_up - w_en))))^2 - aux_en.tke) - end - - Shear² = center_aux_turbconv(state).Shear² - ∂qt∂z = center_aux_turbconv(state).∂qt∂z - ∂θl∂z = center_aux_turbconv(state).∂θl∂z - ∂θv∂z = center_aux_turbconv(state).∂θv∂z - ∂qt∂z_sat = center_aux_turbconv(state).∂qt∂z_sat - ∂θl∂z_sat = center_aux_turbconv(state).∂θl∂z_sat - ∂θv∂z_unsat = center_aux_turbconv(state).∂θv∂z_unsat - - ∇0_bcs = (; bottom = CCO.Extrapolate(), top = CCO.Extrapolate()) - If0 = CCO.InterpolateC2F(; ∇0_bcs...) - uvw = face_aux_turbconv(state).uvw - - uₕ_gm = grid_mean_uₕ(state) - w_en = aux_en_f.w - # compute shear - k̂ = center_aux_turbconv(state).k̂ - - # TODO: Will need to be changed with topography - local_geometry = CC.Fields.local_geometry_field(axes(ρ_c)) - @. k̂ = CCG.Contravariant3Vector(CCG.WVector(FT(1)), local_geometry) - Ifuₕ = - CCO.InterpolateC2F(bottom = CCO.Extrapolate(), top = CCO.Extrapolate()) - ∇uvw = CCO.GradientF2C() - @. uvw = C123(Ifuₕ(uₕ_gm)) + C123(CCG.WVector(w_en)) - @. Shear² = LA.norm_sqr(adjoint(∇uvw(uvw)) * k̂) - - ∇c = CCO.DivergenceF2C() - q_tot_en = aux_en.q_tot - θ_liq_ice_en = aux_en.θ_liq_ice - θ_virt_en = aux_en.θ_virt - @. ∂qt∂z = ∇c(wvec(If0(q_tot_en))) - @. ∂θl∂z = ∇c(wvec(If0(θ_liq_ice_en))) - @. ∂θv∂z = ∇c(wvec(If0(θ_virt_en))) - - # Second order approximation: Use dry and cloudy environmental fields. - cf = aux_en.cloud_fraction - shm = copy(cf) - pshm = parent(shm) - shrink_mask!(pshm, vec(cf)) - - # Since NaN*0 ≠ 0, we need to conditionally replace - # our gradients by their default values. - @. ∂qt∂z_sat = ∇c(wvec(If0(aux_en_sat.q_tot))) - @. ∂θl∂z_sat = ∇c(wvec(If0(aux_en_sat.θ_liq_ice))) - @. ∂θv∂z_unsat = ∇c(wvec(If0(aux_en_unsat.θ_virt))) - @. ∂qt∂z_sat = ifelse(shm == 0, ∂qt∂z, ∂qt∂z_sat) - @. ∂θl∂z_sat = ifelse(shm == 0, ∂θl∂z, ∂θl∂z_sat) - @. ∂θv∂z_unsat = ifelse(shm == 0, ∂θv∂z, ∂θv∂z_unsat) - - bg = center_aux_turbconv(state).buoy_grad - - # buoyancy_gradients - if edmf.bg_closure isa BuoyGradMean - # First order approximation: Use environmental mean fields. - @. bg = buoyancy_gradients( - param_set, - edmf.moisture_model, - EnvBuoyGrad( - edmf.bg_closure, - aux_en.T, # t_sat - TD.vapor_specific_humidity(thermo_params, ts_env), # qv_sat - aux_en.q_tot, # qt_sat - aux_en.θ_dry, # θ_sat - aux_en.θ_liq_ice, # θ_liq_ice_sat - ∂θv∂z, # ∂θv∂z_unsat - ∂qt∂z, # ∂qt∂z_sat - ∂θl∂z, # ∂θl∂z_sat - p_c, # p - aux_en.cloud_fraction, # en_cld_frac - ρ_c, # ρ - ), - ) - elseif edmf.bg_closure isa BuoyGradQuadratures - @. bg = buoyancy_gradients( - param_set, - edmf.moisture_model, - EnvBuoyGrad( - edmf.bg_closure, - aux_en_sat.T, # t_sat - aux_en_sat.q_vap, # qv_sat - aux_en_sat.q_tot, # qt_sat - aux_en_sat.θ_dry, # θ_sat - aux_en_sat.θ_liq_ice, # θ_liq_ice_sat - ∂θv∂z_unsat, # ∂θv∂z_unsat - ∂qt∂z_sat, # ∂qt∂z_sat - ∂θl∂z_sat, # ∂θl∂z_sat - p_c, # p - aux_en.cloud_fraction, # en_cld_frac - ρ_c, # ρ - ), - ) - else - error( - "Something went wrong. The buoyancy gradient model is not specified", - ) - end - - # Limiting stratification scale (Deardorff, 1976) - # compute ∇Ri and Pr - @. aux_tc.prandtl_nvec = turbulent_Prandtl_number( - param_set, - oblength, - gradient_Richardson_number(param_set, bg.∂b∂z, Shear², FT(eps(FT))), - ) - - tke_surf = aux_en.tke[kc_surf] - - @. aux_tc.mixing_length = mixing_length( - param_set, - surface_conditions.ustar, - zc, - tke_surf, - bg.∂b∂z, - aux_en.tke, - oblength, - Shear², - aux_tc.prandtl_nvec, - b_exch, - ) - - @. KM = c_m * aux_tc.mixing_length * sqrt(max(aux_en.tke, 0)) - @. KH = KM / aux_tc.prandtl_nvec - - compute_diffusive_fluxes(edmf, grid, state, param_set) - - # compute precipitation formation tendencies from updrafts - @. aux_bulk.e_tot_tendency_precip_formation = 0 - @. aux_bulk.qt_tendency_precip_formation = 0 - @. aux_bulk.qr_tendency_precip_formation = 0 - @. aux_bulk.qs_tendency_precip_formation = 0 - @inbounds for i in 1:N_up - # autoconversion and accretion - @. mph = precipitation_formation( - param_set, - edmf.precip_model, - prog_gm, - aux_up[i].area, - zc, - Δt, - aux_up[i].ts, - ) - @. aux_up[i].θ_liq_ice_tendency_precip_formation = - mph.θ_liq_ice_tendency * aux_up[i].area - @. aux_up[i].e_tot_tendency_precip_formation = - mph.e_tot_tendency * aux_up[i].area - @. aux_up[i].qt_tendency_precip_formation = - mph.qt_tendency * aux_up[i].area - @. aux_up[i].qr_tendency_precip_formation = - mph.qr_tendency * aux_up[i].area - @. aux_up[i].qs_tendency_precip_formation = - mph.qs_tendency * aux_up[i].area - - # TODO - to be deleted once we sum all tendencies elsewhere - @. aux_bulk.e_tot_tendency_precip_formation += - aux_up[i].e_tot_tendency_precip_formation - @. aux_bulk.qt_tendency_precip_formation += - aux_up[i].qt_tendency_precip_formation - @. aux_bulk.qr_tendency_precip_formation += - aux_up[i].qr_tendency_precip_formation - @. aux_bulk.qs_tendency_precip_formation += - aux_up[i].qs_tendency_precip_formation - end - - return nothing -end diff --git a/src/TurbulenceConvection_deprecated/utility_functions.jl b/src/TurbulenceConvection_deprecated/utility_functions.jl deleted file mode 100644 index 279e45bbbfe..00000000000 --- a/src/TurbulenceConvection_deprecated/utility_functions.jl +++ /dev/null @@ -1,88 +0,0 @@ -"""Compute average between two percentiles of a standard Gaussian""" -function percentile_bounds_mean_norm( - low_percentile::FT, - high_percentile::FT, -) where {FT <: Real} - D = Distributions - gauss_int(x) = -exp(-x * x / 2) / sqrt(2 * pi) - xp_low = D.quantile(D.Normal(), low_percentile) - xp_high = D.quantile(D.Normal(), high_percentile) - return (gauss_int(xp_high) - gauss_int(xp_low)) / - (D.cdf(D.Normal(), xp_high) - D.cdf(D.Normal(), xp_low)) -end - -function logistic(x, slope, mid) - return 1 / (1 + exp(-slope * (x - mid))) -end - -# lambert_2_over_e(::Type{FT}) where {FT} = FT(LambertW.lambertw(FT(2) / FT(MathConstants.e))) -lambert_2_over_e(::Type{FT}) where {FT} = FT(0.46305551336554884) # since we can evaluate - -function lamb_smooth_minimum( - l::SA.SVector, - lower_bound::FT, - upper_bound::FT, -) where {FT} - x_min = minimum(l) - λ_0 = max(x_min * lower_bound / lambert_2_over_e(FT), upper_bound) - - num = sum(l_i -> l_i * exp(-(l_i - x_min) / λ_0), l) - den = sum(l_i -> exp(-(l_i - x_min) / λ_0), l) - smin = num / den - return smin -end - -mean_nc_data(data, group, var, imin, imax) = - StatsBase.mean(data.group[group][var][:][:, imin:imax], dims = 2)[:] - -""" - compare(a::Field, a::Field) - compare(a::FieldVector, b::FieldVector) - -Recursively compare two identically structured -`Field`s, or `FieldVector`s, with potentially different -data. If `!(maximum(abs.(parent(a) .- parent(b))) == 0.0)` -for any single field, then `compare` will print out the fields -and `err`. This can be helpful for debugging where and why -two `Field`/`FieldVector`s are different. -""" -function compare(a::FV, b::FV, pn0 = "") where {FV <: CC.Fields.FieldVector} - for pn in propertynames(a) - pa = getproperty(a, pn) - pb = getproperty(b, pn) - compare(pa, pb, "$pn0.$pn") - end -end - -function compare(a::F, b::F, pn0 = "") where {F <: CC.Fields.Field} - if isempty(propertynames(a)) - err = abs.(parent(a) .- parent(b)) - if !(maximum(err) == 0.0) - println("--- Comparing field $pn0") - @show a - @show b - @show err - end - else - for pn in propertynames(a) - pa = getproperty(a, pn) - pb = getproperty(b, pn) - compare(pa, pb, "$pn0.$pn") - end - end -end - -function conditional_assign!( - v, - ts, - f::F, - thermo_params, - cond::CF = identity, -) where {F, CF} - @. v = ifelse( - cond(TD.has_condensate(thermo_params, ts)), - f(thermo_params, ts), - v, - ) - nothing -end diff --git a/src/TurbulenceConvection_deprecated/variables.jl b/src/TurbulenceConvection_deprecated/variables.jl deleted file mode 100644 index 9a7057afa67..00000000000 --- a/src/TurbulenceConvection_deprecated/variables.jl +++ /dev/null @@ -1,137 +0,0 @@ -##### -##### Fields -##### -import ClimaCore.Geometry: ⊗ - -# Helpers for adding empty thermodynamic state fields: -thermo_state_zeros(::DryModel, FT) = zero(TD.PhaseDry{FT}) -thermo_state_zeros(::EquilMoistModel, FT) = zero(TD.PhaseEquil{FT}) -thermo_state_zeros(::NonEquilMoistModel, FT) = zero(TD.PhaseNonEquil{FT}) - -##### Auxiliary fields - -# Center only -cent_aux_vars_up(FT, local_geometry, atmos) = (; - ts = thermo_state_zeros(atmos.moisture_model, FT), - q_liq = FT(0), - q_ice = FT(0), - T = FT(0), - RH = FT(0), - area = FT(0), - q_tot = FT(0), - θ_liq_ice = FT(0), - e_tot = FT(0), - e_kin = FT(0), - h_tot = FT(0), - θ_liq_ice_tendency_precip_formation = FT(0), - e_tot_tendency_precip_formation = FT(0), - qt_tendency_precip_formation = FT(0), - qr_tendency_precip_formation = FT(0), - qs_tendency_precip_formation = FT(0), - detr = FT(0), - entr = FT(0), -) -cent_aux_vars_edmf(::Type{FT}, local_geometry, atmos) where {FT} = (; - turbconv = (; - ϕ_temporary = FT(0), - ψ_temporary = FT(0), - k̂ = CCG.Contravariant3Vector(CCG.WVector(FT(1)), local_geometry), - bulk = (; - area = FT(0), - h_tot = FT(0), - q_tot = FT(0), - q_liq = FT(0), - q_ice = FT(0), - T = FT(0), - cloud_fraction = FT(0), - e_tot_tendency_precip_formation = FT(0), - qt_tendency_precip_formation = FT(0), - qr_tendency_precip_formation = FT(0), - qs_tendency_precip_formation = FT(0), - filter_flag_1 = FT(0), # tmp flag for testing filters - filter_flag_2 = FT(0), # tmp flag for testing filters - filter_flag_3 = FT(0), # tmp flag for testing filters - filter_flag_4 = FT(0), # tmp flag for testing filters - ), - up = ntuple( - i -> cent_aux_vars_up(FT, local_geometry, atmos), - Val(n_updrafts(atmos.turbconv_model)), - ), - en = (; - ts = thermo_state_zeros(atmos.moisture_model, FT), - area = FT(0), - q_tot = FT(0), - q_liq = FT(0), - q_ice = FT(0), - θ_liq_ice = FT(0), - e_tot = FT(0), - e_kin = FT(0), - h_tot = FT(0), - θ_virt = FT(0), - θ_dry = FT(0), - RH = FT(0), - T = FT(0), - cloud_fraction = FT(0), - tke = FT(0), - #θ_liq_ice_tendency_precip_formation = FT(0), - e_tot_tendency_precip_formation = FT(0), - qt_tendency_precip_formation = FT(0), - qr_tendency_precip_formation = FT(0), - qs_tendency_precip_formation = FT(0), - unsat = (; q_tot = FT(0), θ_dry = FT(0), θ_virt = FT(0)), - sat = (; - T = FT(0), - q_vap = FT(0), - q_tot = FT(0), - θ_dry = FT(0), - θ_liq_ice = FT(0), - ), - ), - mph = PrecipFormation{FT}(), - buoy_grad = GradBuoy{FT}(), - KM = FT(0), - KH = FT(0), - mixing_length = FT(0), - prandtl_nvec = FT(0), - # Variable Prandtl number initialized as neutral value. - b_exch = FT(0), - w_up_c = FT(0), - w_en_c = FT(0), - Shear² = FT(0), - ∂θv∂z = FT(0), - ∂qt∂z = FT(0), - ∂θl∂z = FT(0), - ∂θv∂z_unsat = FT(0), - ∂qt∂z_sat = FT(0), - ∂θl∂z_sat = FT(0), - ) -) - -# Face only -face_aux_vars_up(FT, local_geometry) = (; - nh_pressure = FT(0), - massflux = CCG.Covariant3Vector(FT(0)), - buoy = FT(0), -) -face_aux_vars_edmf(::Type{FT}, local_geometry, edmf) where {FT} = (; - turbconv = (; - bulk = (; w = CCG.Covariant3Vector(FT(0)), buoy_up1 = FT(0)), - ρ_ae_KM = FT(0), - ρ_ae_KH = FT(0), - ρ_ae_K = FT(0), - en = (; w = CCG.Covariant3Vector(FT(0)), buoy = FT(0)), - up = ntuple( - i -> face_aux_vars_up(FT, local_geometry), - Val(n_updrafts(edmf)), - ), - massflux = CCG.Covariant3Vector(FT(0)), - massflux_h = CCG.Covariant3Vector(FT(0)), - massflux_qt = CCG.Covariant3Vector(FT(0)), - ϕ_temporary = CCG.Covariant3Vector(FT(0)), - diffusive_flux_h = CCG.Covariant3Vector(FT(0)), - diffusive_flux_qt = CCG.Covariant3Vector(FT(0)), - diffusive_flux_uₕ = CCG.Covariant3Vector(FT(0)) ⊗ - CCG.Covariant12Vector(FT(0), FT(0)), - uvw = CCG.Covariant123Vector(CCG.WVector(FT(0)), local_geometry), - ) -) diff --git a/src/cache/cache.jl b/src/cache/cache.jl index 94607c0a2b8..aea5faa3820 100644 --- a/src/cache/cache.jl +++ b/src/cache/cache.jl @@ -163,13 +163,6 @@ function additional_cache( ), edmfx_nh_pressure_cache(Y, atmos.turbconv_model), (; Δt = dt), - turbconv_cache( - Y, - turbconv_model, - atmos, - params, - parsed_args, - initial_condition, - ), + (; turbconv_model), ) end diff --git a/src/cache/precomputed_quantities.jl b/src/cache/precomputed_quantities.jl index 88e8e95d0cd..2be28412bcd 100644 --- a/src/cache/precomputed_quantities.jl +++ b/src/cache/precomputed_quantities.jl @@ -358,9 +358,9 @@ values of the first updraft. function output_diagnostic_sgs_quantities(Y, p, t) thermo_params = CAP.thermodynamics_params(p.params) (; ᶜρaʲs, ᶜtsʲs) = p - ᶠu³⁺ = p.ᶠu³ʲs[1] + ᶠu³⁺ = p.ᶠu³ʲs.:1 ᶜu⁺ = @. (C123(Y.c.uₕ) + C123(ᶜinterp(ᶠu³⁺))) - ᶜts⁺ = @. ᶜtsʲs[1] - ᶜa⁺ = @. draft_area(ᶜρaʲs[1], TD.air_density(thermo_params, ᶜts⁺)) + ᶜts⁺ = @. ᶜtsʲs.:1 + ᶜa⁺ = @. draft_area(ᶜρaʲs.:1, TD.air_density(thermo_params, ᶜts⁺)) return (; ᶜu⁺, ᶠu³⁺, ᶜts⁺, ᶜa⁺) end diff --git a/src/callbacks/callbacks.jl b/src/callbacks/callbacks.jl index ff13a84e647..c229a11449c 100644 --- a/src/callbacks/callbacks.jl +++ b/src/callbacks/callbacks.jl @@ -49,28 +49,6 @@ function flux_accumulation!(integrator) return nothing end -NVTX.@annotate function turb_conv_affect_filter!(integrator) - p = integrator.p - (; edmf_cache) = p - (; edmf, param_set) = edmf_cache - t = integrator.t - Y = integrator.u - tc_params = CAP.turbconv_params(param_set) - - set_precomputed_quantities!(Y, p, t) # sets ᶜts for set_edmf_surface_bc - Fields.bycolumn(axes(Y.c)) do colidx - state = TC.tc_column_state(Y, p, nothing, colidx, t) - grid = TC.Grid(state) - TC.affect_filter!(edmf, grid, state, tc_params, t) - end - - # We're lying to OrdinaryDiffEq.jl, in order to avoid - # paying for an additional `∑tendencies!` call, which is required - # to support supplying a continuous representation of the - # solution. - SciMLBase.u_modified!(integrator, false) -end - NVTX.@annotate function rrtmgp_model_callback!(integrator) Y = integrator.u p = integrator.p @@ -351,60 +329,6 @@ NVTX.@annotate function compute_diagnostics(integrator) ᶜts, ) .+ ᶜa⁺ .* cloud_fraction.(thermo_params, ᶜts⁺, ᶜa⁺), ) - elseif p.atmos.turbconv_model isa TC.EDMFModel - tc_cent(p) = p.edmf_cache.aux.c.turbconv - tc_face(p) = p.edmf_cache.aux.f.turbconv - turbulence_convection_diagnostic = (; - bulk_up_area = tc_cent(p).bulk.area, - bulk_up_h_tot = tc_cent(p).bulk.h_tot, - bulk_up_q_tot = tc_cent(p).bulk.q_tot, - bulk_up_q_liq = tc_cent(p).bulk.q_liq, - bulk_up_q_ice = tc_cent(p).bulk.q_ice, - bulk_up_temperature = tc_cent(p).bulk.T, - bulk_up_cloud_fraction = tc_cent(p).bulk.cloud_fraction, - bulk_up_e_tot_tendency_precip_formation = tc_cent( - p, - ).bulk.e_tot_tendency_precip_formation, - bulk_up_qt_tendency_precip_formation = tc_cent( - p, - ).bulk.qt_tendency_precip_formation, - env_area = tc_cent(p).en.area, - env_q_tot = tc_cent(p).en.q_tot, - env_q_liq = tc_cent(p).en.q_liq, - env_q_ice = tc_cent(p).en.q_ice, - env_theta_liq_ice = tc_cent(p).en.θ_liq_ice, - env_theta_virt = tc_cent(p).en.θ_virt, - env_theta_dry = tc_cent(p).en.θ_dry, - env_e_tot = tc_cent(p).en.e_tot, - env_e_kin = tc_cent(p).en.e_kin, - env_h_tot = tc_cent(p).en.h_tot, - env_RH = tc_cent(p).en.RH, - env_temperature = tc_cent(p).en.T, - env_cloud_fraction = tc_cent(p).en.cloud_fraction, - env_TKE = tc_cent(p).en.tke, - env_e_tot_tendency_precip_formation = tc_cent( - p, - ).en.e_tot_tendency_precip_formation, - env_qt_tendency_precip_formation = tc_cent( - p, - ).en.qt_tendency_precip_formation, - face_env_buoyancy = tc_face(p).en.buoy, - face_up1_buoyancy = tc_face(p).bulk.buoy_up1, - face_bulk_w = tc_face(p).bulk.w, - face_env_w = tc_face(p).en.w, - bulk_up_filter_flag_1 = tc_cent(p).bulk.filter_flag_1, - bulk_up_filter_flag_2 = tc_cent(p).bulk.filter_flag_2, - bulk_up_filter_flag_3 = tc_cent(p).bulk.filter_flag_3, - bulk_up_filter_flag_4 = tc_cent(p).bulk.filter_flag_4, - env_q_vap = tc_cent(p).en.q_tot .- tc_cent(p).en.q_liq .- - tc_cent(p).en.q_ice, - draft_q_vap = tc_cent(p).bulk.q_tot .- tc_cent(p).bulk.q_liq .- - tc_cent(p).bulk.q_ice, - cloud_fraction = tc_cent(p).en.area .* - tc_cent(p).en.cloud_fraction .+ - tc_cent(p).bulk.area .* - tc_cent(p).bulk.cloud_fraction, - ) else turbulence_convection_diagnostic = NamedTuple() end diff --git a/src/dycore_equations_deprecated/sgs_flux_tendencies.jl b/src/dycore_equations_deprecated/sgs_flux_tendencies.jl deleted file mode 100644 index 9822e681c2b..00000000000 --- a/src/dycore_equations_deprecated/sgs_flux_tendencies.jl +++ /dev/null @@ -1,245 +0,0 @@ -using LinearAlgebra -import OrdinaryDiffEq as ODE -import Logging -import TerminalLoggers -import LinearAlgebra as LA -import LinearAlgebra: × -import Thermodynamics as TD -import CLIMAParameters as CP - -import ClimaCore.Fields as Fields -import ClimaCore.Spaces as Spaces -import ClimaCore.Operators as Operators -import ClimaCore.Geometry: ⊗ - -import ..Parameters as CAP -import ..TurbulenceConvection as TC -import ..InitialConditions as ICs -import ..TurbulenceConvection.Parameters as TCP -import ..TurbulenceConvection.Parameters.AbstractTurbulenceConvectionParameters as APS - -##### -##### TurbulenceConvection sgs flux tendencies and cache -##### - -function assign_thermo_aux!(state, moisture_model, thermo_params) - If = Operators.InterpolateC2F( - bottom = Operators.Extrapolate(), - top = Operators.Extrapolate(), - ) - aux_gm = TC.center_aux_grid_mean(state) - aux_gm_f = TC.face_aux_grid_mean(state) - prog_gm = TC.center_prog_grid_mean(state) - ᶜts = TC.center_aux_grid_mean_ts(state) - p_c = TC.center_aux_grid_mean_p(state) - ρ_c = prog_gm.ρ - ρ_f = aux_gm_f.ρ - @. ρ_f = If(ρ_c) - - @. aux_gm.q_tot = TD.total_specific_humidity(thermo_params, ᶜts) - @. aux_gm.q_liq = TD.liquid_specific_humidity(thermo_params, ᶜts) - @. aux_gm.q_ice = TD.ice_specific_humidity(thermo_params, ᶜts) - @. aux_gm.T = TD.air_temperature(thermo_params, ᶜts) - @. aux_gm.RH = TD.relative_humidity(thermo_params, ᶜts) - @. aux_gm.θ_liq_ice = TD.liquid_ice_pottemp(thermo_params, ᶜts) - @. aux_gm.h_tot = - TD.total_specific_enthalpy(thermo_params, ᶜts, prog_gm.ρe_tot / ρ_c) - @. p_c = TD.air_pressure(thermo_params, ᶜts) - @. aux_gm.θ_virt = TD.virtual_pottemp(thermo_params, ᶜts) - @. aux_gm.e_tot = prog_gm.ρe_tot / ρ_c - return -end - -function compute_implicit_gm_tendencies!( - edmf::TC.EDMFModel, - grid::TC.Grid, - state::TC.State, - param_set::APS, -) - tendencies_gm = TC.center_tendencies_grid_mean(state) - prog_gm = TC.center_prog_grid_mean(state) - aux_gm_f = TC.face_aux_grid_mean(state) - ρ_c = prog_gm.ρ - tendencies_gm_uₕ = TC.tendencies_grid_mean_uₕ(state) - - TC.compute_sgs_flux!(edmf, grid, state, param_set) - - ∇sgs = Operators.DivergenceF2C() - @. tendencies_gm.ρe_tot += -∇sgs(aux_gm_f.sgs_flux_h_tot) - @. tendencies_gm_uₕ += -∇sgs(aux_gm_f.sgs_flux_uₕ) / ρ_c - if hasproperty(tendencies_gm, :ρq_tot) - @. tendencies_gm.ρq_tot += -∇sgs(aux_gm_f.sgs_flux_q_tot) - end - - return nothing -end - -##### -##### No TurbulenceConvection scheme -##### - -turbconv_cache( - Y, - turbconv_model, - atmos, - param_set, - parsed_args, - initial_condition, -) = (; turbconv_model) - -implicit_sgs_flux_tendency!(Yₜ, Y, p, t, colidx, _) = nothing -explicit_sgs_flux_tendency!(Yₜ, Y, p, t, colidx, _) = nothing - -function turbconv_aux(atmos, edmf, Y, ::Type{FT}) where {FT} - face_aux_vars(FT, local_geometry, atmos, edmf) = (; - sgs_flux_h_tot = Geometry.Covariant3Vector(FT(0)), - sgs_flux_q_tot = Geometry.Covariant3Vector(FT(0)), - sgs_flux_uₕ = Geometry.Covariant3Vector(FT(0)) ⊗ - Geometry.Covariant12Vector(FT(0), FT(0)), - ρ = FT(0), - buoy = FT(0), - TC.face_aux_vars_edmf(FT, local_geometry, edmf)..., - ) - # Center only - cent_aux_vars(FT, local_geometry, atmos, edmf) = (; - tke = FT(0), - q_liq = FT(0), - q_ice = FT(0), - RH = FT(0), - T = FT(0), - buoy = FT(0), - cloud_fraction = FT(0), - θ_virt = FT(0), - Ri = FT(0), - θ_liq_ice = FT(0), - q_tot = FT(0), - h_tot = FT(0), - e_tot = FT(0), - TC.cent_aux_vars_edmf(FT, local_geometry, atmos)..., - ) - fspace = axes(Y.f) - cspace = axes(Y.c) - ᶜlocal_geometry = Fields.local_geometry_field(cspace) - ᶠlocal_geometry = Fields.local_geometry_field(fspace) - - aux_cent_fields = cent_aux_vars.(FT, ᶜlocal_geometry, atmos, edmf) - aux_face_fields = face_aux_vars.(FT, ᶠlocal_geometry, atmos, edmf) - aux = Fields.FieldVector(; c = aux_cent_fields, f = aux_face_fields) - return aux -end - -function turbconv_cache( - Y, - turbconv_model::TC.EDMFModel, - atmos, - param_set, - parsed_args, - initial_condition, -) - FT = Spaces.undertype(axes(Y.c)) - imex_edmf_turbconv = parsed_args["imex_edmf_turbconv"] - imex_edmf_gm = parsed_args["imex_edmf_gm"] - test_consistency = parsed_args["test_edmf_consistency"] - thermo_params = CAP.thermodynamics_params(param_set) - edmf = turbconv_model - @info "EDMFModel: \n$(summary(edmf))" - cache = (; - edmf, - turbconv_model, - imex_edmf_turbconv, - imex_edmf_gm, - test_consistency, - param_set, - aux = turbconv_aux(atmos, edmf, Y, FT), - Y_filtered = similar(Y), - ) - return (; edmf_cache = cache, turbconv_model) -end - -# TODO: Split update_aux! and other functions into implicit and explicit parts. - -function implicit_sgs_flux_tendency!(Yₜ, Y, p, t, colidx, ::TC.EDMFModel) - (; edmf_cache, Δt) = p - (; edmf, param_set, Y_filtered) = edmf_cache - (; imex_edmf_turbconv, imex_edmf_gm, test_consistency) = edmf_cache - thermo_params = CAP.thermodynamics_params(param_set) - tc_params = CAP.turbconv_params(param_set) - - imex_edmf_turbconv || imex_edmf_gm || return nothing - - Y_filtered.c[colidx] .= Y.c[colidx] - Y_filtered.f[colidx] .= Y.f[colidx] - - state = TC.tc_column_state(Y_filtered, p, Yₜ, colidx, t) - - grid = TC.Grid(state) - if test_consistency - parent(state.aux.f) .= NaN - parent(state.aux.c) .= NaN - end - - assign_thermo_aux!(state, edmf.moisture_model, thermo_params) - - TC.affect_filter!(edmf, grid, state, tc_params, t) - - TC.update_aux!(edmf, grid, state, tc_params, t, Δt) - - imex_edmf_turbconv && - TC.compute_implicit_turbconv_tendencies!(edmf, grid, state) - - imex_edmf_gm && - compute_implicit_gm_tendencies!(edmf, grid, state, tc_params) - - # Note: The "filter relaxation tendency" should not be included in the - # implicit tendency because its derivative with respect to Y is - # discontinuous, which means that including it would make the linear - # linear equation being solved by Newton's method ill-conditioned. - - return nothing -end - -function explicit_sgs_flux_tendency!(Yₜ, Y, p, t, colidx, ::TC.EDMFModel) - (; edmf_cache, Δt) = p - (; edmf, param_set, Y_filtered) = edmf_cache - (; imex_edmf_turbconv, imex_edmf_gm, test_consistency) = edmf_cache - thermo_params = CAP.thermodynamics_params(param_set) - tc_params = CAP.turbconv_params(param_set) - - # Note: We could also do Y_filtered .= Y further upstream if needed. - Y_filtered.c[colidx] .= Y.c[colidx] - Y_filtered.f[colidx] .= Y.f[colidx] - - state = TC.tc_column_state(Y_filtered, p, Yₜ, colidx, t) - - grid = TC.Grid(state) - if test_consistency - parent(state.aux.f) .= NaN - parent(state.aux.c) .= NaN - end - - assign_thermo_aux!(state, edmf.moisture_model, thermo_params) - - TC.affect_filter!(edmf, grid, state, tc_params, t) - - TC.update_aux!(edmf, grid, state, tc_params, t, Δt) - - # Ensure that, when a tendency is not computed with an IMEX formulation, - # both its implicit and its explicit components are computed here. - - TC.compute_explicit_turbconv_tendencies!(edmf, grid, state) - imex_edmf_turbconv || - TC.compute_implicit_turbconv_tendencies!(edmf, grid, state) - - # TODO: incrementally disable this and enable proper grid mean terms - imex_edmf_gm || - compute_implicit_gm_tendencies!(edmf, grid, state, tc_params) - - # Note: This "filter relaxation tendency" can be scaled down if needed, but - # it must be present in order to prevent Y and Y_filtered from diverging - # during each timestep. - Yₜ_turbconv = TC.field_vector_column(Yₜ, colidx) - Y_filtered_turbconv = TC.field_vector_column(Y_filtered, colidx) - Y_turbconv = TC.field_vector_column(Y, colidx) - Yₜ_turbconv .+= (Y_filtered_turbconv .- Y_turbconv) ./ Δt - return nothing -end diff --git a/src/initial_conditions/atmos_state.jl b/src/initial_conditions/atmos_state.jl index 84865dd0d24..cde0211e438 100644 --- a/src/initial_conditions/atmos_state.jl +++ b/src/initial_conditions/atmos_state.jl @@ -114,21 +114,6 @@ precip_variables(ls, _, ::PerfExperimental) = # or SGS type configurations (initial TKE needed), so we do not need to assert # that there is no TKE when there is no turbconv_model. turbconv_center_variables(ls, ::Nothing, gs_vars) = (;) -function turbconv_center_variables(ls, turbconv_model::TC.EDMFModel, gs_vars) - n = TC.n_updrafts(turbconv_model) - a_draft = max(ls.turbconv_state.draft_area, n * turbconv_model.minimum_area) - ρa = ls.ρ * a_draft / n - ρae_tot = - ρa * ( - TD.internal_energy(ls.thermo_params, ls.thermo_state) + - norm_sqr(ls.velocity) / 2 + - CAP.grav(ls.params) * ls.geometry.coordinates.z - ) - ρaq_tot = ρa * TD.total_specific_humidity(ls.thermo_params, ls.thermo_state) - en = (; ρatke = (ls.ρ - n * ρa) * ls.turbconv_state.tke) - up = ntuple(_ -> (; ρarea = ρa, ρae_tot, ρaq_tot), Val(n)) - return (; turbconv = (; en, up)) -end function turbconv_center_variables(ls, turbconv_model::EDMFX, gs_vars) n = n_mass_flux_subdomains(turbconv_model) a_draft = ls.turbconv_state.draft_area @@ -143,14 +128,6 @@ function turbconv_center_variables(ls, turbconv_model::DiagnosticEDMFX, gs_vars) end turbconv_face_variables(ls, ::Nothing) = (;) -turbconv_face_variables(ls, turbconv_model::TC.EDMFModel) = (; - turbconv = (; - up = ntuple( - _ -> (; w = C3(zero(eltype(ls)))), - Val(TC.n_updrafts(turbconv_model)), - ) - ) -) turbconv_face_variables(ls, turbconv_model::EDMFX) = (; sgsʲs = ntuple( _ -> (; u₃ = C3(ls.turbconv_state.velocity, ls.geometry)), diff --git a/src/parameterized_tendencies/microphysics/precipitation.jl b/src/parameterized_tendencies/microphysics/precipitation.jl index 68fd2e0ef47..b8a62694d71 100644 --- a/src/parameterized_tendencies/microphysics/precipitation.jl +++ b/src/parameterized_tendencies/microphysics/precipitation.jl @@ -44,24 +44,6 @@ function compute_precipitation_cache!(Y, p, colidx, ::Microphysics0Moment, _) ) end -function compute_precipitation_cache!( - Y, - p, - colidx, - ::Microphysics0Moment, - ::TC.EDMFModel, -) - (; ᶜS_ρq_tot) = p - qt_tendency_precip_formation_en = - p.edmf_cache.aux.c.turbconv.en.qt_tendency_precip_formation[colidx] - qt_tendency_precip_formation_bulk = - p.edmf_cache.aux.c.turbconv.bulk.qt_tendency_precip_formation[colidx] - - @. ᶜS_ρq_tot[colidx] = - Y.c.ρ[colidx] * - (qt_tendency_precip_formation_bulk + qt_tendency_precip_formation_en) -end - function compute_precipitation_cache!( Y, p, @@ -181,118 +163,6 @@ function precipitation_cache(Y, precip_model::Microphysics1Moment) ) end -function compute_precipitation_cache!( - Y, - p, - colidx, - ::Microphysics1Moment, - ::TC.EDMFModel, -) - (; ᶜq_rai, ᶜq_sno) = p - (; ᶜS_ρe_tot, ᶜS_ρq_tot, ᶜS_ρq_rai, ᶜS_ρq_sno) = p - (; ᶜS_q_rai_evap, ᶜS_q_sno_melt, ᶜS_q_sno_sub_dep) = p - (; ᶜts, ᶜΦ, ᶜT, params) = p - (; dt) = p.simulation - - FT = Spaces.undertype(axes(Y.c)) - - # Sources of precipitation from EDMF SGS sub-domains - e_tot_tendency_precip_formation_en = - p.edmf_cache.aux.c.turbconv.en.e_tot_tendency_precip_formation[colidx] - e_tot_tendency_precip_formation_bulk = - p.edmf_cache.aux.c.turbconv.bulk.e_tot_tendency_precip_formation[colidx] - qt_tendency_precip_formation_en = - p.edmf_cache.aux.c.turbconv.en.qt_tendency_precip_formation[colidx] - qt_tendency_precip_formation_bulk = - p.edmf_cache.aux.c.turbconv.bulk.qt_tendency_precip_formation[colidx] - qr_tendency_precip_formation_en = - p.edmf_cache.aux.c.turbconv.en.qr_tendency_precip_formation[colidx] - qr_tendency_precip_formation_bulk = - p.edmf_cache.aux.c.turbconv.bulk.qr_tendency_precip_formation[colidx] - qs_tendency_precip_formation_en = - p.edmf_cache.aux.c.turbconv.en.qs_tendency_precip_formation[colidx] - qs_tendency_precip_formation_bulk = - p.edmf_cache.aux.c.turbconv.bulk.qs_tendency_precip_formation[colidx] - - thermo_params = CAP.thermodynamics_params(params) - cm_params = CAP.microphysics_params(params) - rain_type = CM.CommonTypes.RainType() - snow_type = CM.CommonTypes.SnowType() - @. ᶜT[colidx] = TD.air_temperature(thermo_params, ᶜts[colidx]) - @. ᶜq_rai[colidx] = max(0, Y.c.ρq_rai[colidx] / Y.c.ρ[colidx]) - @. ᶜq_sno[colidx] = max(0, Y.c.ρq_sno[colidx] / Y.c.ρ[colidx]) - - # Sinks of precipitation (evaporation, melting, deposition/sublimation) - # Limiting the tendency by tracer/dt should be handeled in a better way - @. ᶜS_q_rai_evap[colidx] = - -min( - ᶜq_rai[colidx] / dt, - -CM1.evaporation_sublimation( - cm_params, - rain_type, - TD.PhasePartition(thermo_params, ᶜts[colidx]), - ᶜq_rai[colidx], - Y.c.ρ[colidx], - ᶜT[colidx], - ), - ) - @. ᶜS_q_sno_melt[colidx] = - -min( - ᶜq_sno[colidx] / dt, - CM1.snow_melt(cm_params, ᶜq_sno[colidx], Y.c.ρ[colidx], ᶜT[colidx]), - ) - @. ᶜS_q_sno_sub_dep[colidx] = CM1.evaporation_sublimation( - cm_params, - snow_type, - TD.PhasePartition(thermo_params, ᶜts[colidx]), - ᶜq_sno[colidx], - Y.c.ρ[colidx], - ᶜT[colidx], - ) - @. ᶜS_q_sno_sub_dep[colidx] = ifelse( - ᶜS_q_sno_sub_dep[colidx] > 0, - min( - TD.vapor_specific_humidity(thermo_params, ᶜts[colidx]) / dt, - ᶜS_q_sno_sub_dep[colidx], - ), - -min(ᶜq_sno[colidx] / dt, -(ᶜS_q_sno_sub_dep[colidx])), - ) - - # Combine all the sources - @. ᶜS_ρe_tot[colidx] = - Y.c.ρ[colidx] * ( - e_tot_tendency_precip_formation_bulk + - e_tot_tendency_precip_formation_en - - ᶜS_q_rai_evap[colidx] * ( - TD.internal_energy_liquid(thermo_params, ᶜts[colidx]) + - ᶜΦ[colidx] - ) - - ᶜS_q_sno_sub_dep[colidx] * - (TD.internal_energy_ice(thermo_params, ᶜts[colidx]) + ᶜΦ[colidx]) + - ᶜS_q_sno_melt[colidx] * - TD.latent_heat_fusion(thermo_params, ᶜts[colidx]) - ) - @. ᶜS_ρq_tot[colidx] = - Y.c.ρ[colidx] * ( - qt_tendency_precip_formation_bulk + - qt_tendency_precip_formation_en - ᶜS_q_rai_evap[colidx] - - ᶜS_q_sno_sub_dep[colidx] - ) - @. ᶜS_ρq_rai[colidx] = - Y.c.ρ[colidx] * ( - qr_tendency_precip_formation_bulk + - qr_tendency_precip_formation_en + - ᶜS_q_rai_evap[colidx] - ᶜS_q_sno_melt[colidx] - ) - @. ᶜS_ρq_sno[colidx] = - Y.c.ρ[colidx] * ( - qs_tendency_precip_formation_bulk + - qs_tendency_precip_formation_en + - ᶜS_q_sno_sub_dep[colidx] + - ᶜS_q_sno_melt[colidx] - ) -end - """ Computes the rain and snow advection (down) tendency """ diff --git a/src/prognostic_equations/implicit/implicit_tendency.jl b/src/prognostic_equations/implicit/implicit_tendency.jl index de6dadd30a0..039a0e96c47 100644 --- a/src/prognostic_equations/implicit/implicit_tendency.jl +++ b/src/prognostic_equations/implicit/implicit_tendency.jl @@ -10,9 +10,6 @@ NVTX.@annotate function implicit_tendency!(Yₜ, Y, p, t) set_precomputed_quantities!(Y, p, t) Fields.bycolumn(axes(Y.c)) do colidx implicit_vertical_advection_tendency!(Yₜ, Y, p, t, colidx) - if p.turbconv_model isa TurbulenceConvection.EDMFModel - implicit_sgs_flux_tendency!(Yₜ, Y, p, t, colidx, p.turbconv_model) - end # NOTE: All ρa tendencies should be applied before calling this function pressure_work_tendency!(Yₜ, Y, p, t, colidx, p.atmos.turbconv_model) diff --git a/src/prognostic_equations/remaining_tendency.jl b/src/prognostic_equations/remaining_tendency.jl index 6c428fd3890..0fcee804141 100644 --- a/src/prognostic_equations/remaining_tendency.jl +++ b/src/prognostic_equations/remaining_tendency.jl @@ -38,7 +38,6 @@ NVTX.@annotate function additional_tendency!(Yₜ, Y, p, t) edmfx_sgs_flux_tendency!(Yₜ, Y, p, t, colidx, p.atmos.turbconv_model) edmfx_nh_pressure_tendency!(Yₜ, Y, p, t, colidx, p.turbconv_model) edmfx_tke_tendency!(Yₜ, Y, p, t, colidx, p.turbconv_model) - explicit_sgs_flux_tendency!(Yₜ, Y, p, t, colidx, p.turbconv_model) precipitation_tendency!(Yₜ, Y, p, t, colidx, p.precip_model) # NOTE: All ρa tendencies should be applied before calling this function pressure_work_tendency!(Yₜ, Y, p, t, colidx, p.atmos.turbconv_model) diff --git a/src/solver/model_getters.jl b/src/solver/model_getters.jl index 3d7f1a3f1c7..039ed5dfb5d 100644 --- a/src/solver/model_getters.jl +++ b/src/solver/model_getters.jl @@ -343,15 +343,7 @@ function get_turbconv_model( turbconv = parsed_args["turbconv"] @assert turbconv in (nothing, "edmf", "edmfx", "diagnostic_edmfx") - return if turbconv == "edmf" - TC.EDMFModel( - FT, - moisture_model, - precip_model, - parsed_args, - turbconv_params, - ) - elseif turbconv == "edmfx" + return if turbconv == "edmfx" N = turbconv_params.updraft_number TKE = parsed_args["prognostic_tke"] EDMFX{N, TKE}(turbconv_params.min_area) diff --git a/src/solver/type_getters.jl b/src/solver/type_getters.jl index befa7660e6a..2daee76ac94 100644 --- a/src/solver/type_getters.jl +++ b/src/solver/type_getters.jl @@ -526,13 +526,6 @@ function get_callbacks(parsed_args, simulation, atmos, params) (callbacks..., call_every_dt(rrtmgp_model_callback!, dt_rad)) end - if atmos.turbconv_model isa TC.EDMFModel - callbacks = ( - callbacks..., - call_every_n_steps(turb_conv_affect_filter!; skip_first = true), - ) - end - return callbacks end diff --git a/src/surface_conditions/surface_conditions.jl b/src/surface_conditions/surface_conditions.jl index c469e59a65e..809a202fd32 100644 --- a/src/surface_conditions/surface_conditions.jl +++ b/src/surface_conditions/surface_conditions.jl @@ -389,18 +389,10 @@ function atmos_surface_conditions( energy_flux = if atmos.energy_form isa PotentialTemperature (; ρ_flux_θ = shf / TD.cp_m(thermo_params, ts) * surface_normal) elseif atmos.energy_form isa TotalEnergy - if atmos.turbconv_model isa TC.EDMFModel - (; - ρ_flux_h_tot = (shf + lhf) * surface_normal, - ρ_flux_θ = shf / TD.cp_m(thermo_params, ts) * surface_normal, - ) - else - (; ρ_flux_h_tot = (shf + lhf) * surface_normal) - end + (; ρ_flux_h_tot = (shf + lhf) * surface_normal) end moisture_flux = - atmos.moisture_model isa DryModel && - !(atmos.turbconv_model isa TC.EDMFModel) ? (;) : + atmos.moisture_model isa DryModel ? (;) : (; ρ_flux_q_tot = evaporation * surface_normal) return (; ts, @@ -433,15 +425,10 @@ function surface_conditions_type(atmos, ::Type{FT}) where {FT} energy_flux_names = if atmos.energy_form isa PotentialTemperature (:ρ_flux_θ,) elseif atmos.energy_form isa TotalEnergy - if atmos.turbconv_model isa TC.EDMFModel - (:ρ_flux_h_tot, :ρ_flux_θ) - else - (:ρ_flux_h_tot,) - end + (:ρ_flux_h_tot,) end moisture_flux_names = - atmos.moisture_model isa DryModel && - !(atmos.turbconv_model isa TC.EDMFModel) ? () : (:ρ_flux_q_tot,) + atmos.moisture_model isa DryModel ? () : (:ρ_flux_q_tot,) names = ( :ts, :ustar, diff --git a/src/utils/classify_case.jl b/src/utils/classify_case.jl index dcf0743030c..079a9761fa5 100644 --- a/src/utils/classify_case.jl +++ b/src/utils/classify_case.jl @@ -13,19 +13,9 @@ is_solid_body(parsed_args) = all(( parsed_args["perturb_initstate"] == false, )) -is_column_without_edmf(parsed_args) = all(( +is_column_without_edmfx(parsed_args) = all(( parsed_args["config"] == "column", parsed_args["turbconv"] == nothing, parsed_args["forcing"] == nothing, - parsed_args["turbconv"] != "edmf", -)) - -is_column_edmf(parsed_args) = all(( - parsed_args["config"] == "column", - parsed_args["energy_name"] == "rhoe", - parsed_args["forcing"] == nothing, - parsed_args["turbconv"] == "edmf", - parsed_args["rad"] == "DYCOMS_RF01" || - parsed_args["rad"] == "TRMM_LBA" || - parsed_args["rad"] == nothing, + parsed_args["turbconv"] == nothing, )) diff --git a/toml/edmf_trmm.toml b/toml/edmf_trmm.toml deleted file mode 100644 index a941cbf18ea..00000000000 --- a/toml/edmf_trmm.toml +++ /dev/null @@ -1,5 +0,0 @@ -[alpha_rayleigh_uh] -alias = "alpha_rayleigh_uh" -value = 0.0001 -type = "float" -