Fixes allocations and adds ci plot for gcm_driven_scm

post_processing/ci_plots.jl

@@ -42,7 +42,7 @@ import CairoMakie
 import CairoMakie.Makie
 import ClimaAnalysis
 import ClimaAnalysis: Visualize as viz
-import ClimaAnalysis: SimDir, slice, read_var
+import ClimaAnalysis: SimDir, slice, read_var, average_xy
 import ClimaAnalysis.Utils: kwargs as ca_kwargs
 
 import ClimaCoreSpectra: power_spectrum_2d
@@ -1310,3 +1310,50 @@ function make_plots(::EDMFSpherePlots, output_paths::Vector{<:AbstractString})
         MAX_NUM_ROWS = 4,
     )
 end
+
+
+function make_plots(
+    ::Val{:gcm_driven_scm},
+    output_paths::Vector{<:AbstractString},
+)
+    simdirs = SimDir.(output_paths)
+    short_names_2D = [
+        "rlut",
+        "rlutcs",
+        "rsut",
+        "rsutcs",
+        "clwvi",
+        "lwp",
+        "clivi",
+        "dsevi",
+        "clvi",
+        "prw",
+        "hurvi",
+    ]
+    short_names_3D = ["husv", "thetaa", "ta", "hur", "hus", "clw", "cl"]
+    reduction = "inst"
+    vars_2D = map_comparison(simdirs, short_names_2D) do simdir, short_name
+        average_xy(get(simdir; short_name, reduction))
+    end
+    vars_3D = map_comparison(simdirs, short_names_3D) do simdir, short_name
+        data = window(
+            get(simdir; short_name, reduction),
+            "z",
+            left = 0,
+            right = 4000,
+        )
+        return average_xy(data)
+    end
+    make_plots_generic(
+        output_paths,
+        vars_2D;
+        MAX_NUM_COLS = 2,
+        output_name = "summary_2D",
+    )
+    make_plots_generic(
+        output_paths,
+        vars_3D;
+        MAX_NUM_COLS = 2,
+        output_name = "summary_3D",
+    )
+end

src/diagnostics/core_diagnostics.jl

@@ -751,22 +751,22 @@ function compute_clwvi!(
 ) where {T <: Union{EquilMoistModel, NonEquilMoistModel}}
     if isnothing(out)
         out = zeros(axes(Fields.level(state.f, half)))
-        Operators.column_integral_definite!(
-            out,
-            state.c.ρ .* (
-                cache.precomputed.cloud_diagnostics_tuple.q_liq .+
+        clw = cache.scratch.ᶜtemp_scalar
+        @. clw =
+            state.c.ρ * (
+                cache.precomputed.cloud_diagnostics_tuple.q_liq +
                 cache.precomputed.cloud_diagnostics_tuple.q_ice
-            ),
-        )
+            )
+        Operators.column_integral_definite!(out, clw)
         return out
     else
-        Operators.column_integral_definite!(
-            out,
-            state.c.ρ .* (
-                cache.precomputed.cloud_diagnostics_tuple.q_liq .+
+        clw = cache.scratch.ᶜtemp_scalar
+        @. clw =
+            state.c.ρ * (
+                cache.precomputed.cloud_diagnostics_tuple.q_liq +
                 cache.precomputed.cloud_diagnostics_tuple.q_ice
-            ),
-        )
+            )
+        Operators.column_integral_definite!(out, clw)
     end
 end
 
@@ -799,16 +799,14 @@ function compute_lwp!(
 ) where {T <: Union{EquilMoistModel, NonEquilMoistModel}}
     if isnothing(out)
         out = zeros(axes(Fields.level(state.f, half)))
-        Operators.column_integral_definite!(
-            out,
-            state.c.ρ .* (cache.precomputed.cloud_diagnostics_tuple.q_liq),
-        )
+        lw = cache.scratch.ᶜtemp_scalar
+        @. lw = state.c.ρ * cache.precomputed.cloud_diagnostics_tuple.q_liq
+        Operators.column_integral_definite!(out, lw)
         return out
     else
-        Operators.column_integral_definite!(
-            out,
-            state.c.ρ .* (cache.precomputed.cloud_diagnostics_tuple.q_liq),
-        )
+        lw = cache.scratch.ᶜtemp_scalar
+        @. lw = state.c.ρ * cache.precomputed.cloud_diagnostics_tuple.q_liq
+        Operators.column_integral_definite!(out, lw)
     end
 end
 
@@ -841,16 +839,14 @@ function compute_clivi!(
 ) where {T <: Union{EquilMoistModel, NonEquilMoistModel}}
     if isnothing(out)
         out = zeros(axes(Fields.level(state.f, half)))
-        Operators.column_integral_definite!(
-            out,
-            state.c.ρ .* (cache.precomputed.cloud_diagnostics_tuple.q_ice),
-        )
+        cli = cache.scratch.ᶜtemp_scalar
+        @. cli = state.c.ρ * cache.precomputed.cloud_diagnostics_tuple.q_ice
+        Operators.column_integral_definite!(out, cli)
         return out
     else
-        Operators.column_integral_definite!(
-            out,
-            state.c.ρ .* (cache.precomputed.cloud_diagnostics_tuple.q_ice),
-        )
+        cli = cache.scratch.ᶜtemp_scalar
+        @. cli = state.c.ρ * cache.precomputed.cloud_diagnostics_tuple.q_ice
+        Operators.column_integral_definite!(out, cli)
     end
 end
 
@@ -878,7 +874,7 @@ function compute_dsevi!(out, state, cache, time)
         dse = cache.scratch.ᶜtemp_scalar
         @. dse =
             state.c.ρ * (
-                cp * TD.air_temperature.(thermo_params, cache.precomputed.ᶜts) +
+                cp * TD.air_temperature(thermo_params, cache.precomputed.ᶜts) +
                 cache.core.ᶜΦ
             )
         Operators.column_integral_definite!(out, dse)
@@ -888,7 +884,7 @@ function compute_dsevi!(out, state, cache, time)
         dse = cache.scratch.ᶜtemp_scalar
         @. dse =
             state.c.ρ * (
-                cp * TD.air_temperature.(thermo_params, cache.precomputed.ᶜts) +
+                cp * TD.air_temperature(thermo_params, cache.precomputed.ᶜts) +
                 cache.core.ᶜΦ
             )
         Operators.column_integral_definite!(out, dse)
@@ -923,23 +919,22 @@ function compute_clvi!(
     if isnothing(out)
         out = zeros(axes(Fields.level(state.f, half)))
         cloud_cover = cache.scratch.ᶜtemp_scalar
-        cloud_cover .=
-            FT = Spaces.undertype(axes(cloud_cover))
-            ifelse.(
-                cache.precomputed.cloud_diagnostics_tuple.cf .> zero(FT),
-                one(FT),
-                zero(FT),
-            )
+        FT = Spaces.undertype(axes(cloud_cover))
+        @. cloud_cover = ifelse(
+            cache.precomputed.cloud_diagnostics_tuple.cf > zero(FT),
+            one(FT),
+            zero(FT),
+        )
         Operators.column_integral_definite!(out, cloud_cover)
         return out
     else
         cloud_cover = cache.scratch.ᶜtemp_scalar
-        cloud_cover .=
-            ifelse.(
-                cache.precomputed.cloud_diagnostics_tuple.cf .> 0.0f0,
-                1.0f0,
-                0.0f0,
-            )
+        FT = Spaces.undertype(axes(cloud_cover))
+        @. cloud_cover = ifelse(
+            cache.precomputed.cloud_diagnostics_tuple.cf > zero(FT),
+            one(FT),
+            zero(FT),
+        )
         Operators.column_integral_definite!(out, cloud_cover)
     end
 end
@@ -1006,44 +1001,42 @@ function compute_hurvi!(
     thermo_params = CAP.thermodynamics_params(cache.params)
     if isnothing(out)
         out = zeros(axes(Fields.level(state.f, half)))
-        saturation_vi = zeros(axes(Fields.level(state.f, half)))
-        water_vapor_saturation = cache.scratch.ᶜtemp_scalar
-        @. water_vapor_saturation =
+        # compute vertical integral of saturation specific humidity
+        # note next line currently allocates; currently no correct scratch space
+        sat_vi = zeros(axes(Fields.level(state.f, half)))
+        sat = cache.scratch.ᶜtemp_scalar
+        @. sat =
             state.c.ρ *
-            TD.q_vap_saturation.(thermo_params, cache.precomputed.ᶜts)
-        Operators.column_integral_definite!(
-            saturation_vi,
-            water_vapor_saturation,
+            TD.q_vap_saturation(thermo_params, cache.precomputed.ᶜts)
+        Operators.column_integral_definite!(sat_vi, sat)
+        # compute saturation-weighted vertical integral of specific humidity
+        hur = cache.scratch.ᶜtemp_scalar
+        @. hur = TD.vapor_specific_humidity(
+            CAP.thermodynamics_params(cache.params),
+            cache.precomputed.ᶜts,
         )
-        q_vap = cache.scratch.ᶜtemp_scalar
-        q_vap .=
-            TD.vapor_specific_humidity.(
-                CAP.thermodynamics_params(cache.params),
-                cache.precomputed.ᶜts,
-            )
-        saturation = cache.scratch.ᶜtemp_scalar
-        @. saturation = state.c.ρ * q_vap / saturation_vi
-        Operators.column_integral_definite!(out, saturation)
+        hur_weighted = cache.scratch.ᶜtemp_scalar_2
+        @. hur_weighted = state.c.ρ * hur / sat_vi
+        Operators.column_integral_definite!(out, hur_weighted)
         return out
     else
-        saturation_vi = zeros(axes(Fields.level(state.f, half)))
-        water_vapor_saturation = cache.scratch.ᶜtemp_scalar
-        @. water_vapor_saturation =
+        # compute vertical integral of saturation specific humidity
+        # note next line currently allocates; currently no correct scratch space
+        sat_vi = zeros(axes(Fields.level(state.f, half)))
+        sat = cache.scratch.ᶜtemp_scalar
+        @. sat =
             state.c.ρ *
-            TD.q_vap_saturation.(thermo_params, cache.precomputed.ᶜts)
-        Operators.column_integral_definite!(
-            saturation_vi,
-            water_vapor_saturation,
+            TD.q_vap_saturation(thermo_params, cache.precomputed.ᶜts)
+        Operators.column_integral_definite!(sat_vi, sat)
+        # compute saturation-weighted vertical integral of specific humidity
+        hur = cache.scratch.ᶜtemp_scalar
+        @. hur = TD.vapor_specific_humidity(
+            CAP.thermodynamics_params(cache.params),
+            cache.precomputed.ᶜts,
         )
-        q_vap = cache.scratch.ᶜtemp_scalar
-        q_vap .=
-            TD.vapor_specific_humidity.(
-                CAP.thermodynamics_params(cache.params),
-                cache.precomputed.ᶜts,
-            )
-        saturation = cache.scratch.ᶜtemp_scalar
-        @. saturation = state.c.ρ * q_vap / saturation_vi
-        Operators.column_integral_definite!(out, saturation)
+        hur_weighted = cache.scratch.ᶜtemp_scalar_2
+        @. hur_weighted = state.c.ρ * hur / sat_vi
+        Operators.column_integral_definite!(out, hur_weighted)
     end
 end