add edmf diagnostics

CliMA · Oct 11, 2023 · 08658c8 · 08658c8
1 parent 36429e9
commit 08658c8
Show file tree

Hide file tree

Showing 6 changed files with 1,001 additions and 8 deletions.
diff --git a/config/model_configs/diagnostic_edmfx_bomex_box.yml b/config/model_configs/diagnostic_edmfx_bomex_box.yml
@@ -26,4 +26,59 @@ z_max: 3e3
 z_stretch: false
 dt: "10secs" 
 t_end: "6hours" 
-dt_save_to_disk: "10mins"
+dt_save_to_disk: "10mins"
+diagnostics:
+  - short_name: ts
+    period: 10mins
+  - short_name: ta
+    period: 10mins
+  - short_name: thetaa
+    period: 10mins
+  - short_name: ha
+    period: 10mins
+  - short_name: pfull
+    period: 10mins
+  - short_name: rhoa
+    period: 10mins
+  - short_name: ua
+    period: 10mins
+  - short_name: va
+    period: 10mins
+  - short_name: wa
+    period: 10mins
+  - short_name: hur
+    period: 10mins
+  - short_name: hus
+    period: 10mins
+  - short_name: clw
+    period: 10mins
+  - short_name: cli
+    period: 10mins
+  - short_name: hussfc
+    period: 10mins
+  - short_name: evspsbl
+    period: 10mins
+  - short_name: arup
+    period: 10mins
+  - short_name: waup
+    period: 10mins
+  - short_name: taup
+    period: 10mins
+  - short_name: thetaaup
+    period: 10mins
+  - short_name: haup
+    period: 10mins
+  - short_name: husup
+    period: 10mins
+  - short_name: hurup
+    period: 10mins
+  - short_name: clwup
+    period: 10mins
+  - short_name: cliup
+    period: 10mins
+  - short_name: waen
+    period: 10mins
+  - short_name: tke
+    period: 10mins
+  - short_name: lmix
+    period: 10mins
diff --git a/src/diagnostics/Diagnostics.jl b/src/diagnostics/Diagnostics.jl
@@ -32,6 +32,9 @@ import ClimaAtmos.RRTMGPInterface as RRTMGPI
 import ..EDMFX
 import ..DiagnosticEDMFX
 
+# functions used to calculate diagnostics
+import ..draft_area
+
 # 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
@@ -47,14 +47,13 @@ error_diagnostic_variable(variable, model::T) where {T} =
     error_diagnostic_variable("Cannot compute $variable with model = $T")
 
 ###
-# Rho (3d)
+# Density (3d)
 ###
 add_diagnostic_variable!(
     short_name = "rhoa",
     long_name = "Air Density",
     standard_name = "air_density",
     units = "kg m^-3",
-    comments = "Density of air",
     compute! = (out, state, cache, time) -> begin
         if isnothing(out)
             return copy(state.c.ρ)
@@ -129,7 +128,6 @@ add_diagnostic_variable!(
     long_name = "Air Temperature",
     standard_name = "air_temperature",
     units = "K",
-    comments = "Temperature of air",
     compute! = (out, state, cache, time) -> begin
         thermo_params = CAP.thermodynamics_params(cache.params)
         if isnothing(out)
@@ -148,7 +146,6 @@ add_diagnostic_variable!(
     long_name = "Air Potential Temperature",
     standard_name = "air_potential_temperature",
     units = "K",
-    comments = "Potential temperature of air",
     compute! = (out, state, cache, time) -> begin
         thermo_params = CAP.thermodynamics_params(cache.params)
         if isnothing(out)
@@ -159,14 +156,30 @@ add_diagnostic_variable!(
     end,
 )
 
+###
+# Enthalpy (3d)
+###
+add_diagnostic_variable!(
+    short_name = "ha",
+    long_name = "Air Specific Enthalpy",
+    units = "m^2 s^-2",
+    compute! = (out, state, cache, time) -> begin
+        thermo_params = CAP.thermodynamics_params(cache.params)
+        if isnothing(out)
+            return TD.specific_enthalpy.(thermo_params, cache.ᶜts)
+        else
+            out .= TD.specific_enthalpy.(thermo_params, cache.ᶜts)
+        end
+    end,
+)
+
 ###
 # Air pressure (3d)
 ###
 add_diagnostic_variable!(
     short_name = "pfull",
     long_name = "Pressure at Model Full-Levels",
     units = "Pa",
-    comments = "Pressure of air",
     compute! = (out, state, cache, time) -> begin
         if isnothing(out)
             return copy(cache.ᶜp)
@@ -261,6 +274,70 @@ add_diagnostic_variable!(
     compute! = compute_hus!,
 )
 
+###
+# Liquid water specific humidity (3d)
+###
+compute_clw!(out, state, cache, time) =
+    compute_clw!(out, state, cache, time, cache.atmos.moisture_model)
+compute_clw!(_, _, _, _, model::T) where {T} =
+    error_diagnostic_variable("clw", model)
+
+function compute_clw!(
+    out,
+    state,
+    cache,
+    time,
+    moisture_model::T,
+) where {T <: Union{EquilMoistModel, NonEquilMoistModel}}
+    thermo_params = CAP.thermodynamics_params(cache.params)
+    if isnothing(out)
+        return TD.liquid_specific_humidity.(thermo_params, cache.ᶜts)
+    else
+        out .= TD.liquid_specific_humidity.(thermo_params, cache.ᶜts)
+    end
+end
+
+add_diagnostic_variable!(
+    short_name = "clw",
+    long_name = "Mass Fraction of Cloud Liquid Water",
+    standard_name = "mass_fraction_of_cloud_liquid_water_in_air",
+    units = "kg kg^-1",
+    comments = "Includes both large-scale and convective cloud. This is calculated as the mass of cloud liquid water in the grid cell divided by the mass of air (including the water in all phases) in the grid cells. ",
+    compute! = compute_clw!,
+)
+
+###
+# Ice water specific humidity (3d)
+###
+compute_cli!(out, state, cache, time) =
+    compute_cli!(out, state, cache, time, cache.atmos.moisture_model)
+compute_cli!(_, _, _, _, model::T) where {T} =
+    error_diagnostic_variable("cli", model)
+
+function compute_cli!(
+    out,
+    state,
+    cache,
+    time,
+    moisture_model::T,
+) where {T <: Union{EquilMoistModel, NonEquilMoistModel}}
+    thermo_params = CAP.thermodynamics_params(cache.params)
+    if isnothing(out)
+        return TD.ice_specific_humidity.(thermo_params, cache.ᶜts)
+    else
+        out .= TD.ice_specific_humidity.(thermo_params, cache.ᶜts)
+    end
+end
+
+add_diagnostic_variable!(
+    short_name = "cli",
+    long_name = "Mass Fraction of Cloud Ice",
+    standard_name = "mass_fraction_of_cloud_ice_in_air",
+    units = "kg kg^-1",
+    comments = "Includes both large-scale and convective cloud. This is calculated as the mass of cloud ice in the grid cell divided by the mass of air (including the water in all phases) in the grid cell.",
+    compute! = compute_cli!,
+)
+
 ###
 # Surface specific humidity (2d)
 ###

diff --git a/src/diagnostics/default_diagnostics.jl b/src/diagnostics/default_diagnostics.jl
@@ -192,7 +192,7 @@ hourly_average(short_names; output_writer = HDF5Writer()) =
 # Core #
 ########
 function core_default_diagnostics()
-    core_diagnostics = ["ts", "ta", "thetaa", "pfull", "rhoa", "ua", "va", "wa"]
+    core_diagnostics = ["ts", "ta", "thetaa", "ha", "pfull", "rhoa", "ua", "va", "wa"]
 
     return [
         daily_averages(core_diagnostics...)...,
@@ -217,7 +217,7 @@ end
 function default_diagnostics(
     ::T,
 ) where {T <: Union{EquilMoistModel, NonEquilMoistModel}}
-    moist_diagnostics = ["hur", "hus", "hussfc", "evspsbl"]
+    moist_diagnostics = ["hur", "hus", "clw", "cli", "hussfc", "evspsbl"]
 
     return [daily_averages(moist_diagnostics...)...]
 end
@@ -246,3 +246,55 @@ function default_diagnostics(::RRTMGPI.AllSkyRadiationWithClearSkyDiagnostics)
 
     return [daily_averages(clear_diagnostics...)...]
 end
+
+##################
+# Turbconv model #
+##################
+function default_diagnostics(::EDMFX)
+    edmfx_diagnostics = [
+        "arup",
+        "rhoaup",
+        "waup",
+        "taup",
+        "thetaaup",
+        "haup",
+        "husup",
+        "hurup",
+        "clwup",
+        "cliup",
+        "aren",
+        "rhoaen",
+        "waen",
+        "taen",
+        "thetaaen",
+        "husen",
+        "huren",
+        "clwen",
+        "clien",
+        "tke",
+        "lmix",
+    ]
+
+    return [daily_averages(edmfx_diagnostics...)...]
+end
+
+
+function default_diagnostics(::DiagnosticEDMFX)
+    diagnostic_edmfx_diagnostics = [
+        "arup",
+        "rhoaup",
+        "waup",
+        "taup",
+        "thetaaup",
+        "haup",
+        "husup",
+        "hurup",
+        "clwup",
+        "cliup",
+        "waen",
+        "tke",
+        "lmix",
+    ]
+
+    return [daily_averages(diagnostic_edmfx_diagnostics...)...]
+end
diff --git a/src/diagnostics/diagnostic.jl b/src/diagnostics/diagnostic.jl
@@ -182,6 +182,7 @@ end
 # Do you want to define more diagnostics? Add them here
 include("core_diagnostics.jl")
 include("radiation_diagnostics.jl")
+include("edmfx_diagnostics.jl")
 
 # Default diagnostics and higher level interfaces
 include("default_diagnostics.jl")