From 76dc06d6dd439ad24bbeb35e7e3554e71b2bee56 Mon Sep 17 00:00:00 2001 From: Anna Jaruga Date: Wed, 13 Dec 2023 16:28:17 -0800 Subject: [PATCH] add mixing length to diagnostics --- ...uarez_rhoe_equilmoist_topography_dcmip.yml | 15 ++++++ src/cache/cloud_fraction.jl | 54 +++++++++++-------- src/diagnostics/Diagnostics.jl | 1 + src/diagnostics/core_diagnostics.jl | 21 +++++++- src/diagnostics/netcdf_writer.jl | 1 + 5 files changed, 68 insertions(+), 24 deletions(-) diff --git a/config/model_configs/sphere_held_suarez_rhoe_equilmoist_topography_dcmip.yml b/config/model_configs/sphere_held_suarez_rhoe_equilmoist_topography_dcmip.yml index d41ed96e0b0..793d9dfb9a3 100644 --- a/config/model_configs/sphere_held_suarez_rhoe_equilmoist_topography_dcmip.yml +++ b/config/model_configs/sphere_held_suarez_rhoe_equilmoist_topography_dcmip.yml @@ -11,3 +11,18 @@ topography: "DCMIP200" job_id: "sphere_held_suarez_rhoe_equilmoist_topography_dcmip" moist: "equil" netcdf_interpolate_z_over_msl: true +diagnostics: + - short_name: cl + period: 3hours + - short_name: cli + period: 3hours + - short_name: clw + period: 3hours + - short_name: ta + period: 3hours + - short_name: hus + period: 3hours + - short_name: hur + period: 3hours + - short_name: mixlgm + period: 3hours diff --git a/src/cache/cloud_fraction.jl b/src/cache/cloud_fraction.jl index 79d9e83123d..f08ae9d2614 100644 --- a/src/cache/cloud_fraction.jl +++ b/src/cache/cloud_fraction.jl @@ -18,7 +18,7 @@ end - Pr - Prandtl number - ϵ_st - strain rate norm """ -function compute_smagorinsky_length_scale(c_smag, N_eff, dz, Pr, ϵ_st) +function smagorinsky_length_scale(c_smag, N_eff, dz, Pr, ϵ_st) FT = eltype(c_smag) #return N_eff > FT(0) && N_eff < sqrt(2 * Pr * ϵ_st) ? return N_eff > FT(0) ? @@ -26,20 +26,14 @@ function compute_smagorinsky_length_scale(c_smag, N_eff, dz, Pr, ϵ_st) c_smag * dz end -""" - Compute the grid scale cloud fraction based on sub-grid scale properties -""" -function set_cloud_fraction!(Y, p, ::DryModel) - @. p.precomputed.ᶜcloud_fraction = 0 -end -function set_cloud_fraction!(Y, p, ::Union{EquilMoistModel, NonEquilMoistModel}) - (; SG_quad, params) = p +function compute_gm_mixing_length!(ᶜmixing_length, Y, p) + (; params) = p + thermo_params = CAP.thermodynamics_params(params) FT = eltype(params) - thermo_params = CAP.thermodynamics_params(params) ᶜdz = Fields.Δz_field(axes(Y.c)) ᶜlg = Fields.local_geometry_field(Y.c) - (; ᶜts, ᶜp, ᶠu³, ᶜcloud_fraction) = p.precomputed + (; ᶜts, ᶜp, ᶠu³) = p.precomputed (; obukhov_length) = p.precomputed.sfc_conditions ᶜlinear_buoygrad = p.scratch.ᶜtemp_scalar @@ -78,26 +72,40 @@ function set_cloud_fraction!(Y, p, ::Union{EquilMoistModel, NonEquilMoistModel}) ᶠu = p.scratch.ᶠtemp_C123 @. ᶠu = C123(ᶠinterp(Y.c.uₕ)) + C123(ᶠu³) - ᶜstrain_rate = p.scratch.ᶜtemp_UVWxUVW compute_strain_rate_center!(ᶜstrain_rate, ᶠu) - ᶜprandtl_nvec = p.scratch.ᶜtemp_scalar - @. ᶜprandtl_nvec = turbulent_prandtl_number( - params, - obukhov_length, - ᶜlinear_buoygrad, - norm_sqr(ᶜstrain_rate), - ) - - ᶜl_smag = p.scratch.ᶜtemp_scalar_2 - @. ᶜl_smag = compute_smagorinsky_length_scale( + @. ᶜmixing_length = smagorinsky_length_scale( CAP.c_smag(params), sqrt(max(ᶜlinear_buoygrad, 0)), #N_eff ᶜdz, - ᶜprandtl_nvec, + turbulent_prandtl_number( + params, + obukhov_length, + ᶜlinear_buoygrad, + norm_sqr(ᶜstrain_rate), + ), norm_sqr(ᶜstrain_rate), ) +end + +""" + Compute the grid scale cloud fraction based on sub-grid scale properties +""" +function set_cloud_fraction!(Y, p, ::DryModel) + @. p.precomputed.ᶜcloud_fractcompute_smagorinsky_length_scaleion = 0 +end +function set_cloud_fraction!(Y, p, ::Union{EquilMoistModel, NonEquilMoistModel}) + (; SG_quad, params) = p + + FT = eltype(params) + thermo_params = CAP.thermodynamics_params(params) + (; ᶜts, ᶜp, ᶜcloud_fraction) = p.precomputed + + # temp_scalar is already in use inside the compute_gm_mixing_length + # this is not a great pattern, but I don't have a more elegant idea + ᶜl_smag = p.scratch.ᶜtemp_scalar_2 + compute_gm_mixing_length!(ᶜl_smag, Y, p) coeff = FT(2.1) # TODO - move to parameters @. ᶜcloud_fraction = quad_loop( diff --git a/src/diagnostics/Diagnostics.jl b/src/diagnostics/Diagnostics.jl index 562778626f7..51d36696085 100644 --- a/src/diagnostics/Diagnostics.jl +++ b/src/diagnostics/Diagnostics.jl @@ -35,6 +35,7 @@ import ..DiagnosticEDMFX # functions used to calculate diagnostics import ..draft_area +import ..compute_gm_mixing_length! # We need the abbreviations for symbols like curl, grad, and so on include(joinpath("..", "utils", "abbreviations.jl")) diff --git a/src/diagnostics/core_diagnostics.jl b/src/diagnostics/core_diagnostics.jl index bdf7b9ce600..1add87444ac 100644 --- a/src/diagnostics/core_diagnostics.jl +++ b/src/diagnostics/core_diagnostics.jl @@ -4,7 +4,7 @@ # # In addition to the metadata (names, comments, ...), the most important step in adding a # new DiagnosticVariable is defining its compute! function. `compute!` has to take four -# arguments: (out, state, cache, time), and as to write the diagnostic in place into the +# arguments: (out, state, cache, time), and has to write the diagnostic in place into the # `out` variable. # # Often, it is possible to compute certain diagnostics only for specific models (e.g., @@ -226,6 +226,25 @@ add_diagnostic_variable!( end, ) +### +# Grid mean mixing length (3d) +### +function compute_mixlgm!(out, state, cache, time) + if isnothing(out) + out = similar(state.c.ρ) + end + compute_gm_mixing_length!(out, state, cache) +end + +add_diagnostic_variable!( + short_name = "mixlgm", + long_name = "Grid mean mixing length", + standard_name = "grid_mean_mixing_length", + units = "m", + comments = "Smagorinsky length scale, to be used as an approximation of mixing length in simulations without EDMF SGS model", + compute! = compute_mixlgm!, +) + ### # Relative humidity (3d) ### diff --git a/src/diagnostics/netcdf_writer.jl b/src/diagnostics/netcdf_writer.jl index 6c9592ca991..3b22375df43 100644 --- a/src/diagnostics/netcdf_writer.jl +++ b/src/diagnostics/netcdf_writer.jl @@ -597,6 +597,7 @@ function write_field!( deflatelevel = writer.compression_level, ) v.attrib["long_name"] = diagnostic.output_long_name + v.attrib["short_name"] = var.short_name v.attrib["units"] = var.units v.attrib["comments"] = var.comments