diff --git a/.buildkite/gpu_pipeline/pipeline.yml b/.buildkite/gpu_pipeline/pipeline.yml
index 2950b9ffdd1..6c773595dca 100644
--- a/.buildkite/gpu_pipeline/pipeline.yml
+++ b/.buildkite/gpu_pipeline/pipeline.yml
@@ -53,7 +53,7 @@ steps:
       - label: "gpu_aquaplanet_dyamond"
         command:
           - mkdir -p gpu_aquaplanet_dyamond
-          - > 
+          - >
             nsys profile --trace=nvtx,cuda --output=gpu_aquaplanet_dyamond/report
             julia --color=yes --project=examples examples/hybrid/driver.jl
             --config_file ${GPU_CONFIG_PATH}gpu_aquaplanet_dyamond.yml
@@ -63,7 +63,7 @@ steps:
 
       - label: "moist Held-Suarez"
         key: "gpu_hs_rhoe_equilmoist_nz63_0M_55km_rs35km"
-        command: 
+        command:
           - mkdir -p gpu_hs_rhoe_equilmoist_nz63_0M_55km_rs35km
           - >
             nsys profile --trace=nvtx,cuda --output=gpu_hs_rhoe_equilmoist_nz63_0M_55km_rs35km/report
diff --git a/.buildkite/longruns/pipeline.yml b/.buildkite/longruns/pipeline.yml
index 29540b4f184..a61921f0048 100644
--- a/.buildkite/longruns/pipeline.yml
+++ b/.buildkite/longruns/pipeline.yml
@@ -23,7 +23,7 @@ steps:
       - "julia --project -e 'using Pkg; Pkg.instantiate(;verbose=true)'"
       - "julia --project -e 'using Pkg; Pkg.precompile()'"
       - "julia --project -e 'using Pkg; Pkg.status()'"
-      
+
       - echo "--- Configure CUDA"
       # force the initialization of the CUDA runtime as it is lazily loaded by default
       - "julia --project=cuda_env -e 'using Pkg; Pkg.resolve(); Pkg.instantiate(;verbose=true);using CUDA; CUDA.precompile_runtime()'"
@@ -49,57 +49,47 @@ steps:
       - label: ":computer: baroclinic wave (ρe_tot) high resolution"
         command:
           - srun julia --project=examples examples/hybrid/driver.jl --config_file $CONFIG_PATH/$$JOB_NAME.yml
-          - julia --color=yes --project=examples post_processing/remap/remap_pipeline.jl --data_dir $$JOB_NAME --out_dir $$JOB_NAME
-          - julia --color=yes --project=examples post_processing/plot/plot_pipeline.jl --nc_dir $$JOB_NAME --fig_dir $$JOB_NAME --case_name dry_baroclinic_wave
         artifact_paths: "$$JOB_NAME/*"
         agents:
           slurm_ntasks: 32
           slurm_time: 24:00:00
-        env: 
+        env:
           JOB_NAME: "longrun_bw_rhoe_highres"
 
       - label: ":computer: no lim ARS baroclinic wave (ρe_tot) equilmoist high resolution"
         command:
           - srun julia --project=examples examples/hybrid/driver.jl --config_file $CONFIG_PATH/$$JOB_NAME.yml
-          - julia --color=yes --project=examples post_processing/remap/remap_pipeline.jl --data_dir $$JOB_NAME --out_dir $$JOB_NAME
-          - julia --color=yes --project=examples post_processing/plot/plot_pipeline.jl --nc_dir $$JOB_NAME --fig_dir $$JOB_NAME --case_name moist_baroclinic_wave
         artifact_paths: "$$JOB_NAME/*"
         agents:
           slurm_ntasks: 32
           slurm_time: 24:00:00
-        env: 
+        env:
           JOB_NAME: "longrun_bw_rhoe_equil_highres"
 
       - label: ":computer: lim ARS zalesak baroclinic wave (ρe_tot) equilmoist high resolution"
         command:
           - srun julia --project=examples examples/hybrid/driver.jl --config_file $CONFIG_PATH/$$JOB_NAME.yml
-          - julia --color=yes --project=examples post_processing/remap/remap_pipeline.jl --data_dir $$JOB_NAME --out_dir $$JOB_NAME
-          - julia --color=yes --project=examples post_processing/plot/plot_pipeline.jl --nc_dir $$JOB_NAME --fig_dir $$JOB_NAME --case_name moist_baroclinic_wave
         artifact_paths: "$$JOB_NAME/*"
         agents:
           slurm_ntasks: 32
           slurm_mem_per_cpu: 32GB
           slurm_time: 24:00:00
-        env: 
+        env:
           JOB_NAME: "longrun_zalesak_tracer_energy_bw_rhoe_equil_highres"
 
       - label: ":computer: SSP baroclinic wave (ρe_tot) equilmoist high resolution centered diff"
         command:
           - "srun julia --project=examples examples/hybrid/driver.jl --config_file $CONFIG_PATH/$$JOB_NAME.yml"
-          - julia --color=yes --project=examples post_processing/remap/remap_pipeline.jl --data_dir $$JOB_NAME --out_dir $$JOB_NAME
-          - julia --color=yes --project=examples post_processing/plot/plot_pipeline.jl --nc_dir $$JOB_NAME --fig_dir $$JOB_NAME --case_name moist_baroclinic_wave
         artifact_paths: "$$JOB_NAME/*"
         agents:
           slurm_ntasks: 32
           slurm_time: 24:00:00
         env:
           JOB_NAME: "longrun_ssp_bw_rhoe_equil_highres"
-      
+
       - label: ":computer: held-suarez, dry, high-topped (55km), high-sponge (35km), helem_16 np_3"
         command:
           - srun julia --project=examples examples/hybrid/driver.jl --config_file $CONFIG_PATH/$$JOB_NAME.yml
-          - julia --color=yes --project=examples post_processing/remap/remap_pipeline.jl --data_dir $$JOB_NAME --out_dir $$JOB_NAME
-          - julia --color=yes --project=examples post_processing/plot/plot_pipeline.jl --nc_dir $$JOB_NAME --fig_dir $$JOB_NAME --case_name dry_held_suarez
         artifact_paths: "$$JOB_NAME/*"
         env:
           CLIMACORE_DISTRIBUTED: "MPI"
@@ -112,8 +102,6 @@ steps:
       - label: ":computer: held-suarez, equilmoist, high-topped (55km), high-sponge (35km), helem_16 np_3"
         command:
           - srun julia --project=examples examples/hybrid/driver.jl --config_file $CONFIG_PATH/$$JOB_NAME.yml
-          - julia --color=yes --project=examples post_processing/remap/remap_pipeline.jl --data_dir $$JOB_NAME --out_dir $$JOB_NAME
-          - julia --color=yes --project=examples post_processing/plot/plot_pipeline.jl --nc_dir $$JOB_NAME --fig_dir $$JOB_NAME --case_name aquaplanet
         artifact_paths: "$$JOB_NAME/*"
         env:
           CLIMACORE_DISTRIBUTED: "MPI"
@@ -122,12 +110,10 @@ steps:
           slurm_ntasks: 64
           slurm_mem_per_cpu: 16GB
           slurm_time: 24:00:00
-      
+
       - label: ":computer: aquaplanet, equilmoist, high-topped (55km), gray-radiation, vertdiff, high-sponge (35km), helem_16 np_3"
         command:
           - srun julia --project=examples examples/hybrid/driver.jl --config_file $CONFIG_PATH/$$JOB_NAME.yml
-          - julia --color=yes --project=examples post_processing/remap/remap_pipeline.jl --data_dir $$JOB_NAME --out_dir $$JOB_NAME
-          - julia --color=yes --project=examples post_processing/plot/plot_pipeline.jl --nc_dir $$JOB_NAME --fig_dir $$JOB_NAME --case_name aquaplanet
         artifact_paths: "$$JOB_NAME/*"
         env:
           CLIMACORE_DISTRIBUTED: "MPI"
@@ -140,8 +126,6 @@ steps:
       - label: ":computer: aquaplanet (ρe_tot) equilmoist clearsky high resolution (nz63) hightop (55km) rayleigh sponge(35e3, 10) Float32 (time-varying insolation)"
         command:
           - srun julia --project=examples examples/hybrid/driver.jl --config_file $CONFIG_PATH/$$JOB_NAME.yml
-          - julia --color=yes --project=examples post_processing/remap/remap_pipeline.jl --data_dir $$JOB_NAME --out_dir $$JOB_NAME
-          - julia --color=yes --project=examples post_processing/plot/plot_pipeline.jl --nc_dir $$JOB_NAME --fig_dir $$JOB_NAME --case_name aquaplanet
         artifact_paths: "$$JOB_NAME/*"
         agents:
           slurm_ntasks: 64
@@ -152,50 +136,44 @@ steps:
       - label: ":computer: aquaplanet (ρe_tot) equilmoist clearsky high resolution (nz63) hightop (55km) rayleigh sponge(35e3, 10) Float32 (time-varying insolation) + earth topography"
         command:
           - srun julia --project=examples examples/hybrid/driver.jl --config_file $CONFIG_PATH/$$JOB_NAME.yml
-          - julia --color=yes --project=examples post_processing/remap/remap_pipeline.jl --data_dir $$JOB_NAME --out_dir $$JOB_NAME
-          - julia --color=yes --project=examples post_processing/plot/plot_pipeline.jl --nc_dir $$JOB_NAME --fig_dir $$JOB_NAME --case_name aquaplanet
         artifact_paths: "$$JOB_NAME/*"
         agents:
           slurm_ntasks: 64
           slurm_time: 24:00:00
         env:
           JOB_NAME: "longrun_aquaplanet_rhoe_equilmoist_nz63_0M_55km_rs35km_clearsky_tvinsolation_earth"
-      
+
       - label: ":computer: baroclinic wave (ρe_tot) equilmoist high resolution topography (earth)"
         command:
           - srun julia --project=examples examples/hybrid/driver.jl --config_file $CONFIG_PATH/$$JOB_NAME.yml
-          - julia --color=yes --project=examples post_processing/remap/remap_pipeline.jl --data_dir $$JOB_NAME --out_dir $$JOB_NAME
         artifact_paths: "$$JOB_NAME/*"
         agents:
           slurm_ntasks: 32
           slurm_mem_per_cpu: 16GB
           slurm_time: 24:00:00
-        env: 
+        env:
           JOB_NAME: "longrun_bw_rhoe_equil_highres_topography_earth"
 
       - label: ":computer: held suarez (ρe_tot) equilmoist high resolution topography (earth)"
         command:
           - srun julia --project=examples examples/hybrid/driver.jl --config_file $CONFIG_PATH/$$JOB_NAME.yml
-          - julia --color=yes --project=examples post_processing/remap/remap_pipeline.jl --data_dir $$JOB_NAME --out_dir $$JOB_NAME
         artifact_paths: "$$JOB_NAME/*"
         agents:
           slurm_ntasks: 64
           slurm_mem_per_cpu: 16GB
           slurm_time: 24:00:00
-        env: 
+        env:
           JOB_NAME: "longrun_hs_rhoe_equil_highres_topography_earth"
 
       - label: ":computer: no lim ARS aquaplanet (ρe_tot) equilmoist high resolution clearsky radiation Float32 (earth)"
         command:
           - srun julia --project=examples examples/hybrid/driver.jl --config_file $CONFIG_PATH/$$JOB_NAME.yml
-          - julia --color=yes --project=examples post_processing/remap/remap_pipeline.jl --data_dir $$JOB_NAME --out_dir $$JOB_NAME
-          - julia --color=yes --project=examples post_processing/plot/plot_pipeline.jl --nc_dir $$JOB_NAME --fig_dir $$JOB_NAME --case_name aquaplanet
         artifact_paths: "$$JOB_NAME/*"
         agents:
           slurm_ntasks: 64
           slurm_mem_per_cpu: 16GB
           slurm_time: 24:00:00
-        env: 
+        env:
           JOB_NAME: "longrun_aquaplanet_rhoe_equil_highres_clearsky_ft32_earth"
 
   - group: "Low resolution long runs"
@@ -205,15 +183,12 @@ steps:
       - label: ":computer: hydrostatic balance (ρe_tot)"
         command:
           - julia --color=yes --project=examples --threads=8 examples/hybrid/driver.jl --config_file $CONFIG_PATH/$$JOB_NAME.yml
-          - julia --color=yes --project=examples post_processing/remap/remap_pipeline.jl --data_dir $$JOB_NAME --out_dir $$JOB_NAME
-          - julia --color=yes --project=examples post_processing/plot/plot_pipeline.jl --nc_dir $$JOB_NAME --fig_dir $$JOB_NAME --case_name dry_held_suarez
         artifact_paths: "$$JOB_NAME/*"
         agents:
           slurm_cpus_per_task: 8
           slurm_time: 24:00:00
         env:
           JOB_NAME: "longrun_sphere_hydrostatic_balance_rhoe"
-      
 
   - group: "Experimental Long runs"
 
@@ -222,7 +197,6 @@ steps:
       - label: ":computer: aquaplanet (ρe_tot) equilmoist high resolution allsky radiation float64"
         command:
           - srun julia --project=examples examples/hybrid/driver.jl --config_file $CONFIG_PATH/$$JOB_NAME.yml
-          - julia --color=yes --project=examples post_processing/remap/remap_pipeline.jl --data_dir $$JOB_NAME --out_dir $$JOB_NAME
         artifact_paths: "$$JOB_NAME/*"
         agents:
           slurm_ntasks: 32
@@ -237,7 +211,6 @@ steps:
       - label: ":computer: aquaplanet dyamond"
         command:
           - srun julia --project=examples examples/hybrid/driver.jl --config_file $CONFIG_PATH/$$JOB_NAME.yml
-          - julia --color=yes --project=examples post_processing/remap/remap_pipeline.jl --data_dir $$JOB_NAME --out_dir $$JOB_NAME
         artifact_paths: "$$JOB_NAME/*"
         agents:
           slurm_ntasks: 64
@@ -245,7 +218,7 @@ steps:
           slurm_time: 24:00:00
         env:
           JOB_NAME: "longrun_aquaplanet_dyamond"
-  
+
   - group: "AMIP"
 
     steps:
@@ -253,7 +226,6 @@ steps:
       - label: ":computer: aquaplanet amip"
         command:
           - srun julia --project=examples examples/hybrid/driver.jl --config_file $CONFIG_PATH/$$JOB_NAME.yml
-          - julia --color=yes --project=examples post_processing/remap/remap_pipeline.jl --data_dir $$JOB_NAME --out_dir $$JOB_NAME
         artifact_paths: "$$JOB_NAME/*"
         agents:
           slurm_ntasks: 64
diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml
index 9df2b2ff7d0..49db0523b2d 100644
--- a/.buildkite/pipeline.yml
+++ b/.buildkite/pipeline.yml
@@ -196,15 +196,6 @@ steps:
         command: >
           julia --color=yes --project=examples examples/hybrid/driver.jl
           --config_file $CONFIG_PATH/sphere_hydrostatic_balance_rhoe_ft64.yml
-
-          julia --color=yes --project=examples post_processing/remap/remap_pipeline.jl
-          --data_dir sphere_hydrostatic_balance_rhoe_ft64
-          --out_dir sphere_hydrostatic_balance_rhoe_ft64
-
-          julia --color=yes --project=examples post_processing/plot/plot_pipeline.jl
-          --nc_dir sphere_hydrostatic_balance_rhoe_ft64
-          --fig_dir sphere_hydrostatic_balance_rhoe_ft64
-          --case_name dry_held_suarez
         artifact_paths: "sphere_hydrostatic_balance_rhoe_ft64/*"
         agents:
           slurm_mem: 20GB
@@ -214,45 +205,18 @@ steps:
         command: >
           julia --color=yes --project=examples examples/hybrid/driver.jl
           --config_file $CONFIG_PATH/sphere_baroclinic_wave_rhoe.yml
-
-          julia --color=yes --project=examples post_processing/remap/remap_pipeline.jl
-          --data_dir sphere_baroclinic_wave_rhoe --out_dir sphere_baroclinic_wave_rhoe
-
-          julia --color=yes --project=examples post_processing/plot/plot_pipeline.jl
-          --nc_dir sphere_baroclinic_wave_rhoe --fig_dir sphere_baroclinic_wave_rhoe
-          --case_name dry_baroclinic_wave
         artifact_paths: "sphere_baroclinic_wave_rhoe/*"
 
       - label: ":computer: no lim ARS baroclinic wave (ρe) equilmoist"
         command: >
           julia --color=yes --project=examples examples/hybrid/driver.jl
           --config_file $CONFIG_PATH/sphere_baroclinic_wave_rhoe_equilmoist.yml
-
-          julia --color=yes --project=examples post_processing/remap/remap_pipeline.jl
-          --data_dir sphere_baroclinic_wave_rhoe_equilmoist
-          --out_dir sphere_baroclinic_wave_rhoe_equilmoist
-
-          julia --color=yes --project=examples post_processing/plot/plot_pipeline.jl
-          --nc_dir sphere_baroclinic_wave_rhoe_equilmoist
-          --fig_dir sphere_baroclinic_wave_rhoe_equilmoist --case_name moist_baroclinic_wave
-
-          julia --color=yes --project=examples regression_tests/test_mse.jl
-          --job_id sphere_baroclinic_wave_rhoe_equilmoist
-          --out_dir sphere_baroclinic_wave_rhoe_equilmoist
         artifact_paths: "sphere_baroclinic_wave_rhoe_equilmoist/*"
-      
+
       - label: ":computer: no lim ARS baroclinic wave (ρe) equilmoist explicit vertdiff"
         command: >
           julia --color=yes --project=examples examples/hybrid/driver.jl
           --config_file $CONFIG_PATH/sphere_baroclinic_wave_rhoe_equilmoist_expvdiff.yml
-
-          julia --color=yes --project=examples post_processing/remap/remap_pipeline.jl
-          --data_dir sphere_baroclinic_wave_rhoe_equilmoist_expvdiff
-          --out_dir sphere_baroclinic_wave_rhoe_equilmoist_expvdiff
-
-          julia --color=yes --project=examples post_processing/plot/plot_pipeline.jl
-          --nc_dir sphere_baroclinic_wave_rhoe_equilmoist_expvdiff
-          --fig_dir sphere_baroclinic_wave_rhoe_equilmoist_expvdiff --case_name aquaplanet
         artifact_paths: "sphere_baroclinic_wave_rhoe_equilmoist_expvdiff/*"
 
       # Add this back when we figure out what to do with SSP and zalesak
@@ -260,15 +224,6 @@ steps:
       #   command: >
       #     julia --color=yes --project=examples examples/hybrid/driver.jl
       #     --config_file $CONFIG_PATH/$$JOB_NAME.yml
-
-      #     julia --color=yes --project=examples post_processing/remap/remap_pipeline.jl
-      #     --data_dir $$JOB_NAME
-      #     --out_dir $$JOB_NAME
-
-      #     julia --color=yes --project=examples post_processing/plot/plot_pipeline.jl
-      #     --nc_dir $$JOB_NAME
-      #     --fig_dir $$JOB_NAME --case_name moist_baroclinic_wave
-
       #   artifact_paths: "$$JOB_NAME/*"
       #   agents:
       #     slurm_mem: 64GB
@@ -283,13 +238,6 @@ steps:
         command: >
           julia --color=yes --project=examples examples/hybrid/driver.jl
           --config_file $CONFIG_PATH/sphere_held_suarez_rhotheta.yml
-
-          julia --color=yes --project=examples post_processing/remap/remap_pipeline.jl
-          --data_dir sphere_held_suarez_rhotheta --out_dir sphere_held_suarez_rhotheta
-
-          julia --color=yes --project=examples post_processing/plot/plot_pipeline.jl
-          --nc_dir sphere_held_suarez_rhotheta --fig_dir sphere_held_suarez_rhotheta
-          --case_name dry_held_suarez
         artifact_paths: "sphere_held_suarez_rhotheta/*"
 
       - label: ":computer: held suarez (ρe) hightop"
@@ -297,13 +245,6 @@ steps:
         command: >
           julia --color=yes --project=examples examples/hybrid/driver.jl
           --config_file $CONFIG_PATH/sphere_held_suarez_rhoe_hightop.yml
-
-          julia --color=yes --project=examples post_processing/remap/remap_pipeline.jl
-          --data_dir sphere_held_suarez_rhoe_hightop --out_dir sphere_held_suarez_rhoe_hightop
-
-          julia --color=yes --project=examples post_processing/plot/plot_pipeline.jl
-          --nc_dir sphere_held_suarez_rhoe_hightop --fig_dir sphere_held_suarez_rhoe_hightop
-          --case_name dry_held_suarez
         artifact_paths: "sphere_held_suarez_rhoe_hightop/*"
 
       - label: ":computer: no lim ARS held suarez (ρe) equilmoist hightop sponge"
@@ -311,14 +252,6 @@ steps:
           julia --color=yes --project=examples examples/hybrid/driver.jl
           --config_file $CONFIG_PATH/sphere_held_suarez_rhoe_equilmoist_hightop_sponge.yml
 
-          julia --color=yes --project=examples post_processing/remap/remap_pipeline.jl
-          --data_dir sphere_held_suarez_rhoe_equilmoist_hightop_sponge
-          --out_dir sphere_held_suarez_rhoe_equilmoist_hightop_sponge
-
-          julia --color=yes --project=examples post_processing/plot/plot_pipeline.jl
-          --nc_dir sphere_held_suarez_rhoe_equilmoist_hightop_sponge
-          --fig_dir sphere_held_suarez_rhoe_equilmoist_hightop_sponge --case_name moist_held_suarez
-
           julia --color=yes --project=examples regression_tests/test_mse.jl
           --job_id sphere_held_suarez_rhoe_equilmoist_hightop_sponge
           --out_dir sphere_held_suarez_rhoe_equilmoist_hightop_sponge
@@ -331,14 +264,6 @@ steps:
         command: >
           julia --color=yes --project=examples examples/hybrid/driver.jl
           --config_file $CONFIG_PATH/sphere_aquaplanet_rhoe_equilmoist_allsky_gw_res.yml
-
-          julia --color=yes --project=examples post_processing/remap/remap_pipeline.jl
-          --data_dir sphere_aquaplanet_rhoe_equilmoist_allsky_gw_res
-          --out_dir sphere_aquaplanet_rhoe_equilmoist_allsky_gw_res
-
-          julia --color=yes --project=examples post_processing/plot/plot_pipeline.jl
-          --nc_dir sphere_aquaplanet_rhoe_equilmoist_allsky_gw_res
-          --fig_dir sphere_aquaplanet_rhoe_equilmoist_allsky_gw_res --case_name aquaplanet
         artifact_paths: "sphere_aquaplanet_rhoe_equilmoist_allsky_gw_res/*"
         agents:
           slurm_mem: 20GB
@@ -348,15 +273,6 @@ steps:
           julia --color=yes --project=examples examples/hybrid/driver.jl
           --config_file $CONFIG_PATH/sphere_aquaplanet_rhoe_equilmoist_allsky_gw_raw_zonallyasymmetric.yml
 
-          julia --color=yes --project=examples post_processing/remap/remap_pipeline.jl
-          --data_dir sphere_aquaplanet_rhoe_equilmoist_allsky_gw_raw_zonallyasymmetric
-          --out_dir sphere_aquaplanet_rhoe_equilmoist_allsky_gw_raw_zonallyasymmetric
-
-          julia --color=yes --project=examples post_processing/plot/plot_pipeline.jl
-          --nc_dir sphere_aquaplanet_rhoe_equilmoist_allsky_gw_raw_zonallyasymmetric
-          --fig_dir sphere_aquaplanet_rhoe_equilmoist_allsky_gw_raw_zonallyasymmetric
-          --case_name aquaplanet
-
           julia --color=yes --project=examples regression_tests/test_mse.jl
           --job_id sphere_aquaplanet_rhoe_equilmoist_allsky_gw_raw_zonallyasymmetric
           --out_dir sphere_aquaplanet_rhoe_equilmoist_allsky_gw_raw_zonallyasymmetric
@@ -371,30 +287,18 @@ steps:
         command: >
           julia --color=yes --project=examples examples/hybrid/driver.jl
           --config_file $CONFIG_PATH/sphere_baroclinic_wave_rhoe_topography_dcmip_rs.yml
-
-          julia --color=yes --project=examples post_processing/remap/remap_pipeline.jl
-          --data_dir sphere_baroclinic_wave_rhoe_topography_dcmip_rs
-          --out_dir sphere_baroclinic_wave_rhoe_topography_dcmip_rs
         artifact_paths: "sphere_baroclinic_wave_rhoe_topography_dcmip_rs/*"
 
       - label: ":computer: held suarez (ρe) topography (dcmip)"
         command: >
           julia --color=yes --project=examples examples/hybrid/driver.jl
           --config_file $CONFIG_PATH/sphere_held_suarez_rhoe_topography_dcmip.yml
-
-          julia --color=yes --project=examples post_processing/remap/remap_pipeline.jl
-          --data_dir sphere_held_suarez_rhoe_topography_dcmip
-          --out_dir sphere_held_suarez_rhoe_topography_dcmip
         artifact_paths: "sphere_held_suarez_rhoe_topography_dcmip/*"
 
       - label: ":computer: held suarez (ρe) equilmoist topography (dcmip)"
         command: >
           julia --color=yes --project=examples examples/hybrid/driver.jl
           --config_file $CONFIG_PATH/sphere_held_suarez_rhoe_equilmoist_topography_dcmip.yml
-
-          julia --color=yes --project=examples post_processing/remap/remap_pipeline.jl
-          --data_dir sphere_held_suarez_rhoe_equilmoist_topography_dcmip
-          --out_dir sphere_held_suarez_rhoe_equilmoist_topography_dcmip
         artifact_paths: "sphere_held_suarez_rhoe_equilmoist_topography_dcmip/*"
         agents:
           slurm_mem: 20GB
@@ -403,30 +307,12 @@ steps:
         command: >
           julia --color=yes --project=examples examples/hybrid/driver.jl
           --config_file $CONFIG_PATH/sphere_ssp_baroclinic_wave_rhoe_equilmoist_dcmip200.yml
-
-          julia --color=yes --project=examples post_processing/remap/remap_pipeline.jl
-          --data_dir sphere_ssp_baroclinic_wave_rhoe_equilmoist_dcmip200
-          --out_dir sphere_ssp_baroclinic_wave_rhoe_equilmoist_dcmip200
-
-          julia --color=yes --project=examples post_processing/plot/plot_pipeline.jl
-          --nc_dir sphere_ssp_baroclinic_wave_rhoe_equilmoist_dcmip200
-          --fig_dir sphere_ssp_baroclinic_wave_rhoe_equilmoist_dcmip200
-          --case_name elevation_spectrum
         artifact_paths: "sphere_ssp_baroclinic_wave_rhoe_equilmoist_dcmip200/*"
 
       - label: ":computer: Diagnostic Earth surface elevation spectra"
         command: >
           julia --color=yes --project=examples examples/hybrid/driver.jl
           --config_file $CONFIG_PATH/sphere_ssp_baroclinic_wave_rhoe_equilmoist_earth.yml
-
-          julia --color=yes --project=examples post_processing/remap/remap_pipeline.jl
-          --data_dir sphere_ssp_baroclinic_wave_rhoe_equilmoist_earth
-          --out_dir sphere_ssp_baroclinic_wave_rhoe_equilmoist_earth
-
-          julia --color=yes --project=examples post_processing/plot/plot_pipeline.jl
-          --nc_dir sphere_ssp_baroclinic_wave_rhoe_equilmoist_earth
-          --fig_dir sphere_ssp_baroclinic_wave_rhoe_equilmoist_earth
-          --case_name elevation_spectrum
         artifact_paths: "sphere_ssp_baroclinic_wave_rhoe_equilmoist_earth/*"
 
   - group: "MPI Examples"
diff --git a/examples/Manifest.toml b/examples/Manifest.toml
index 00ea63173e0..8bd42e9e126 100644
--- a/examples/Manifest.toml
+++ b/examples/Manifest.toml
@@ -2,7 +2,7 @@
 
 julia_version = "1.9.3"
 manifest_format = "2.0"
-project_hash = "0c14077b57d8fdc7663f22496855895314b0dc22"
+project_hash = "73a755acd2133f41fe6635e6f44ecab666d097c8"
 
 [[deps.ADTypes]]
 git-tree-sha1 = "5d2e21d7b0d8c22f67483ef95ebdc39c0e6b6003"
@@ -301,6 +301,16 @@ weakdeps = ["SparseArrays"]
     [deps.ChainRulesCore.extensions]
     ChainRulesCoreSparseArraysExt = "SparseArrays"
 
+[[deps.ClimaAnalysis]]
+deps = ["NCDatasets", "OrderedCollections", "Statistics"]
+path = "../lib/ClimaAnalysis"
+uuid = "29b5916a-a76c-4e73-9657-3c8fd22e65e6"
+version = "0.1.0"
+weakdeps = ["CairoMakie"]
+
+    [deps.ClimaAnalysis.extensions]
+    CairoMakieExt = "CairoMakie"
+
 [[deps.ClimaAtmos]]
 deps = ["ArgParse", "ArtifactWrappers", "Artifacts", "AtmosphericProfilesLibrary", "CLIMAParameters", "ClimaComms", "ClimaCore", "ClimaTimeSteppers", "CloudMicrophysics", "Colors", "Dates", "Dierckx", "DiffEqBase", "DiffEqCallbacks", "DocStringExtensions", "FastGaussQuadrature", "HDF5_jll", "ImageFiltering", "Insolation", "Interpolations", "IntervalSets", "Krylov", "LinearAlgebra", "Logging", "NCDatasets", "NVTX", "Pkg", "Printf", "RRTMGP", "Random", "RootSolvers", "SciMLBase", "StaticArrays", "Statistics", "StatsBase", "SurfaceFluxes", "Test", "Thermodynamics", "YAML"]
 path = ".."
diff --git a/examples/Project.toml b/examples/Project.toml
index ed627719ed3..a0792404f56 100644
--- a/examples/Project.toml
+++ b/examples/Project.toml
@@ -5,6 +5,7 @@ AtmosphericProfilesLibrary = "86bc3604-9858-485a-bdbe-831ec50de11d"
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
 CLIMAParameters = "6eacf6c3-8458-43b9-ae03-caf5306d3d53"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
+ClimaAnalysis = "29b5916a-a76c-4e73-9657-3c8fd22e65e6"
 ClimaAtmos = "b2c96348-7fb7-4fe0-8da9-78d88439e717"
 ClimaComms = "3a4d1b5c-c61d-41fd-a00a-5873ba7a1b0d"
 ClimaCore = "d414da3d-4745-48bb-8d80-42e94e092884"
diff --git a/examples/hybrid/driver.jl b/examples/hybrid/driver.jl
index 85f42466a15..bf9d939bda4 100644
--- a/examples/hybrid/driver.jl
+++ b/examples/hybrid/driver.jl
@@ -12,6 +12,7 @@ sol_res = CA.solve_atmos!(integrator)
 (; p) = integrator
 
 import ClimaCore
+import ClimaCore: Spaces, Topologies
 import ClimaAtmos.InitialConditions as ICs
 using Statistics: mean
 import ClimaAtmos.Parameters as CAP
@@ -29,25 +30,11 @@ import OrderedCollections
 using ClimaCoreTempestRemap
 using ClimaCorePlots, Plots
 using ClimaCoreMakie, CairoMakie
-include(joinpath(pkgdir(CA), "post_processing", "common_utils.jl"))
-
-include(joinpath(pkgdir(CA), "post_processing", "contours_and_profiles.jl"))
-include(joinpath(pkgdir(CA), "post_processing", "post_processing_funcs.jl"))
-include(
-    joinpath(pkgdir(CA), "post_processing", "define_tc_quicklook_profiles.jl"),
-)
+include(joinpath(pkgdir(CA), "post_processing", "make_plots.jl"))
 
 ref_job_id = config.parsed_args["reference_job_id"]
 reference_job_id = isnothing(ref_job_id) ? simulation.job_id : ref_job_id
 
-is_edmfx =
-    atmos.turbconv_model isa CA.PrognosticEDMFX ||
-    atmos.turbconv_model isa CA.DiagnosticEDMFX
-if is_edmfx && config.parsed_args["post_process"]
-    contours_and_profiles(simulation.output_dir, reference_job_id)
-    zip_and_cleanup_output(simulation.output_dir, "hdf5files.zip")
-end
-
 if sol_res.ret_code == :simulation_crashed
     error(
         "The ClimaAtmos simulation has crashed. See the stack trace for details.",
@@ -58,31 +45,31 @@ end
 @assert last(sol.t) == simulation.t_end
 CA.verify_callbacks(sol.t)
 
-if CA.is_distributed(config.comms_ctx)
-    export_scaling_file(
-        sol,
-        simulation.output_dir,
-        walltime,
-        config.comms_ctx,
-        ClimaComms.nprocs(config.comms_ctx),
-    )
-end
+@info "Plotting"
+make_plots(Val{Symbol(reference_job_id)}, simulation.output_dir)
+@info "Plotting done"
 
-if !CA.is_distributed(config.comms_ctx) &&
-   config.parsed_args["post_process"] &&
-   !is_edmfx &&
-   !(atmos.model_config isa CA.SphericalModel)
-    ENV["GKSwstype"] = "nul" # avoid displaying plots
-    if is_column_without_edmfx(config.parsed_args)
-        custom_postprocessing(sol, simulation.output_dir, p)
-    elseif is_solid_body(config.parsed_args)
-        postprocessing(sol, simulation.output_dir, config.parsed_args["fps"])
-    elseif atmos.model_config isa CA.BoxModel
-        postprocessing_box(sol, simulation.output_dir)
-    elseif atmos.model_config isa CA.PlaneModel
-        postprocessing_plane(sol, simulation.output_dir, p)
-    else
-        error("Uncaught case")
+if CA.is_distributed(config.comms_ctx)
+    nprocs = ClimaComms.nprocs(config.comms_ctx)
+    comms_ctx = config.comms_ctx
+    output_dir = simulation.output_dir
+    # replace sol.u on the root processor with the global sol.u
+    if ClimaComms.iamroot(comms_ctx)
+        Y = sol.u[1]
+        center_space = axes(Y.c)
+        horz_space = Spaces.horizontal_space(center_space)
+        horz_topology = horz_space.topology
+        Nq = Spaces.Quadratures.degrees_of_freedom(horz_space.quadrature_style)
+        nlocalelems = Topologies.nlocalelems(horz_topology)
+        ncols_per_process = nlocalelems * Nq * Nq
+        scaling_file =
+            joinpath(output_dir, "scaling_data_$(nprocs)_processes.jld2")
+        @info(
+            "Writing scaling data",
+            "walltime (seconds)" = walltime,
+            scaling_file
+        )
+        JLD2.jldsave(scaling_file; nprocs, ncols_per_process, walltime)
     end
 end
 
diff --git a/perf/Manifest.toml b/perf/Manifest.toml
index af0f302a8ff..85be9e3af77 100644
--- a/perf/Manifest.toml
+++ b/perf/Manifest.toml
@@ -2,7 +2,7 @@
 
 julia_version = "1.9.3"
 manifest_format = "2.0"
-project_hash = "551bce7c2dee7d8c00ba5b736c87527020a6054b"
+project_hash = "b40a57c3a14df49f94c2fac1d141bb7693ebcff9"
 
 [[deps.ADTypes]]
 git-tree-sha1 = "5d2e21d7b0d8c22f67483ef95ebdc39c0e6b6003"
@@ -312,6 +312,16 @@ weakdeps = ["SparseArrays"]
     [deps.ChainRulesCore.extensions]
     ChainRulesCoreSparseArraysExt = "SparseArrays"
 
+[[deps.ClimaAnalysis]]
+deps = ["NCDatasets", "OrderedCollections", "Statistics"]
+path = "../lib/ClimaAnalysis"
+uuid = "29b5916a-a76c-4e73-9657-3c8fd22e65e6"
+version = "0.1.0"
+weakdeps = ["CairoMakie"]
+
+    [deps.ClimaAnalysis.extensions]
+    CairoMakieExt = "CairoMakie"
+
 [[deps.ClimaAtmos]]
 deps = ["ArgParse", "ArtifactWrappers", "Artifacts", "AtmosphericProfilesLibrary", "CLIMAParameters", "ClimaComms", "ClimaCore", "ClimaTimeSteppers", "CloudMicrophysics", "Colors", "Dates", "Dierckx", "DiffEqBase", "DiffEqCallbacks", "DocStringExtensions", "FastGaussQuadrature", "HDF5_jll", "ImageFiltering", "Insolation", "Interpolations", "IntervalSets", "Krylov", "LinearAlgebra", "Logging", "NCDatasets", "NVTX", "Pkg", "Printf", "RRTMGP", "Random", "RootSolvers", "SciMLBase", "StaticArrays", "Statistics", "StatsBase", "SurfaceFluxes", "Test", "Thermodynamics", "YAML"]
 path = ".."
diff --git a/perf/Project.toml b/perf/Project.toml
index cac2fafedde..5566b01542c 100644
--- a/perf/Project.toml
+++ b/perf/Project.toml
@@ -6,6 +6,7 @@ BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
 CLIMAParameters = "6eacf6c3-8458-43b9-ae03-caf5306d3d53"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
+ClimaAnalysis = "29b5916a-a76c-4e73-9657-3c8fd22e65e6"
 ClimaAtmos = "b2c96348-7fb7-4fe0-8da9-78d88439e717"
 ClimaComms = "3a4d1b5c-c61d-41fd-a00a-5873ba7a1b0d"
 ClimaCore = "d414da3d-4745-48bb-8d80-42e94e092884"
diff --git a/post_processing/common_utils.jl b/post_processing/common_utils.jl
deleted file mode 100644
index 6bcb940988b..00000000000
--- a/post_processing/common_utils.jl
+++ /dev/null
@@ -1,46 +0,0 @@
-import ClimaCore: Spaces, Topologies
-import JLD2
-
-function export_scaling_file(sol, output_dir, walltime, comms_ctx, nprocs)
-    # replace sol.u on the root processor with the global sol.u
-    if ClimaComms.iamroot(comms_ctx)
-        Y = sol.u[1]
-        center_space = axes(Y.c)
-        horz_space = Spaces.horizontal_space(center_space)
-        horz_topology = horz_space.topology
-        Nq = Spaces.Quadratures.degrees_of_freedom(horz_space.quadrature_style)
-        nlocalelems = Topologies.nlocalelems(horz_topology)
-        ncols_per_process = nlocalelems * Nq * Nq
-        scaling_file =
-            joinpath(output_dir, "scaling_data_$(nprocs)_processes.jld2")
-        @info(
-            "Writing scaling data",
-            "walltime (seconds)" = walltime,
-            scaling_file
-        )
-        JLD2.jldsave(scaling_file; nprocs, ncols_per_process, walltime)
-    end
-    return nothing
-end
-
-#TODO - do we want to change anything here now?
-is_baro_wave(parsed_args) = all((
-    parsed_args["config"] == "sphere",
-    parsed_args["forcing"] == nothing,
-    parsed_args["surface_setup"] == nothing,
-    parsed_args["perturb_initstate"] == true,
-))
-
-is_solid_body(parsed_args) = all((
-    parsed_args["config"] == "sphere",
-    parsed_args["forcing"] == nothing,
-    parsed_args["rad"] == nothing,
-    parsed_args["perturb_initstate"] == false,
-))
-
-is_column_without_edmfx(parsed_args) = all((
-    parsed_args["config"] == "column",
-    parsed_args["turbconv"] == nothing,
-    parsed_args["forcing"] == nothing,
-    parsed_args["turbconv"] == nothing,
-))
diff --git a/post_processing/compare_outputs.jl b/post_processing/compare_outputs.jl
index 30e67fbc25d..81810abe069 100644
--- a/post_processing/compare_outputs.jl
+++ b/post_processing/compare_outputs.jl
@@ -22,8 +22,6 @@ function parse_commandline()
         help = "A bool: compare the state between output_folder_1 and output_folder_2 [`true` [default], `false`]"
         arg_type = Bool
         default = true
-        "--compare_diagnostics"
-        help = "A bool: compare the diagnostics between output_folder_1 and output_folder_2 [`true` [default], `false`]"
         arg_type = Bool
         default = true
     end
@@ -45,8 +43,6 @@ function get_data(folder, name)
 end
 Y_1 = get_data(parsed_args["output_folder_1"], "Y");
 Y_2 = get_data(parsed_args["output_folder_2"], "Y");
-D_1 = get_data(parsed_args["output_folder_1"], "diagnostics");
-D_2 = get_data(parsed_args["output_folder_2"], "diagnostics");
 
 compare(a::FV, b::FV, pn0 = "") where {FV <: Fields.FieldVector} =
     compare(true, a, b, pn0)
@@ -83,11 +79,4 @@ end
     else
         @test_broken compare(Y_1, Y_2)
     end
-
-    @info "Comparing diagnostics:"
-    if parsed_args["compare_diagnostics"]
-        @test compare(D_1, D_2)
-    else
-        @test_broken compare(D_1, D_2)
-    end
 end
diff --git a/post_processing/comparison_plots.jl b/post_processing/comparison_plots.jl
deleted file mode 100644
index f9493941050..00000000000
--- a/post_processing/comparison_plots.jl
+++ /dev/null
@@ -1,203 +0,0 @@
-using Plots
-using Plots.PlotMeasures
-using JLD2
-
-job_id = ARGS[1]
-output_dir = "./"
-
-secs_per_hour = 60 * 60
-secs_per_day = 60 * 60 * 24
-days_per_year = 8760 / 24
-
-# tempest scaling data
-if occursin("low", job_id)
-    resolution = "low-resolution"
-    nprocs_tempest = [1, 2, 3, 6, 24]
-    walltime_tempest = [102.5892, 49.0212, 33.41088, 17.064, 5.77584]
-    t_int_days_tempest = 10 # integration time for Tempest
-    t_int_days_climaatmos = 10 # integration time for ClimaAtmos
-else
-    resolution = "high-resolution"
-    nprocs_tempest = [1, 2, 3, 6, 24, 54, 96, 216]
-    walltime_tempest = [
-        6186.595968,
-        3197.915424,
-        2160.269568,
-        1131.526368,
-        385.558272,
-        174.662784,
-        139.771872,
-        45.971712,
-        43.602624,
-    ]
-    t_int_days_tempest = 1 # integration time for Tempest
-    t_int_days_climaatmos = 1 # integration time for ClimaAtmos
-end
-t_int_days = t_int_days_climaatmos
-t_int = string(t_int_days_climaatmos) * " days"
-# normalalize Tempest walltime for comparison with ClimaAtmos (if necessary)
-walltime_tempest .*= (t_int_days / t_int_days_tempest)
-
-# read ClimaAtmos scaling data
-I, FT = Int, Float64
-nprocs_clima_atmos = I[]
-walltime_clima_atmos = FT[]
-for foldername in readdir(output_dir)
-    if occursin(job_id, foldername)
-        nprocs_string = split(foldername, "_")[end]
-        dict = load(
-            joinpath(
-                output_dir,
-                foldername,
-                "scaling_data_$(nprocs_string)_processes.jld2",
-            ),
-        )
-        push!(nprocs_clima_atmos, I(dict["nprocs"]))
-        push!(walltime_clima_atmos, FT(dict["walltime"]))
-    end
-end
-order = sortperm(nprocs_clima_atmos)
-nprocs_clima_atmos, walltime_clima_atmos =
-    nprocs_clima_atmos[order], walltime_clima_atmos[order]
-
-# update data if all data points are not available for ClimaAtmos
-minprocs = nprocs_clima_atmos[1]
-loc = findfirst(nprocs_tempest .== minprocs)
-nprocs_tempest = nprocs_tempest[loc:end]
-walltime_tempest = walltime_tempest[loc:end]
-
-@assert nprocs_clima_atmos == nprocs_tempest # needed for comparison Plots.plot
-
-# simulated years per day
-sypd_tempest = (secs_per_day ./ walltime_tempest) * t_int_days ./ days_per_year
-sypd_clima_atmos =
-    (secs_per_day ./ walltime_clima_atmos) * t_int_days ./ days_per_year
-# CPU hours
-cpu_hours_tempest = nprocs_tempest .* walltime_tempest / secs_per_hour
-cpu_hours_clima_atmos =
-    nprocs_clima_atmos .* walltime_clima_atmos / secs_per_hour
-# scaling efficiency
-single_proc_time_tempest = walltime_tempest[1] * nprocs_tempest[1]
-scaling_efficiency_tempest =
-    trunc.(
-        100 * (single_proc_time_tempest ./ nprocs_tempest) ./ walltime_tempest,
-        digits = 1,
-    )
-single_proc_time_clima_atmos = walltime_clima_atmos[1] * nprocs_clima_atmos[1]
-scaling_efficiency_clima_atmos =
-    trunc.(
-        100 * (single_proc_time_clima_atmos ./ nprocs_clima_atmos) ./
-        walltime_clima_atmos,
-        digits = 1,
-    )
-ENV["GKSwstype"] = "100"
-Plots.GRBackend()
-plt1 = Plots.plot(
-    [log2.(nprocs_clima_atmos) log2.(nprocs_tempest)],
-    [sypd_clima_atmos sypd_tempest],
-    markershape = [:circle :square],
-    markercolor = [:blue :orange],
-    xticks = (
-        log2.(nprocs_clima_atmos),
-        [string(i) for i in nprocs_clima_atmos],
-    ),
-    xlabel = "# of MPI processes",
-    ylabel = "SYPD",
-    yaxis = :log,
-    title = "Simulated years per day",
-    label = ["ClimaAtmos (Float64)" "Tempest (Float64)"],
-    legend = :topleft,
-    grid = :true,
-    left_margin = 10mm,
-    bottom_margin = 10mm,
-    top_margin = 10mm,
-)
-Plots.png(plt1, joinpath(output_dir, resolution * "_" * "sypd"))
-Plots.pdf(plt1, joinpath(output_dir, resolution * "_" * "sypd"))
-
-Plots.GRBackend()
-plt2 = Plots.plot(
-    [log2.(nprocs_clima_atmos) log2.(nprocs_tempest)],
-    [cpu_hours_clima_atmos cpu_hours_tempest],
-    markershape = [:circle :square],
-    markercolor = [:blue :orange],
-    xticks = (
-        log2.(nprocs_clima_atmos),
-        [string(i) for i in nprocs_clima_atmos],
-    ),
-    xlabel = "# of MPI processes",
-    ylabel = "CPU hours",
-    title = "Scaling data (T_int = $t_int)",
-    label = ["ClimaAtmos (Float64)" "Tempest (Float64)"],
-    ylims = (0.0, Inf),
-    legend = :topleft,
-    grid = :true,
-    left_margin = 10mm,
-    bottom_margin = 10mm,
-    top_margin = 10mm,
-)
-Plots.png(plt2, joinpath(output_dir, resolution * "_" * "Scaling"))
-Plots.pdf(plt2, joinpath(output_dir, resolution * "_" * "Scaling"))
-
-Plots.GRBackend()
-plt3 = Plots.plot(
-    log2.(nprocs_clima_atmos),
-    cpu_hours_clima_atmos ./ cpu_hours_tempest,
-    markershape = [:circle :square],
-    markercolor = [:blue :orange],
-    xticks = (
-        log2.(nprocs_clima_atmos),
-        [string(i) for i in nprocs_clima_atmos],
-    ),
-    xlabel = "# of MPI processes",
-    ylabel = "ratio",
-    ylims = (0.0, Inf),
-    label = "Ratio of CPU hours (ClimaAtmos/Tempest) (Float64)",
-    title = "Comparison with Tempest",
-    legend = :topleft,
-    grid = :true,
-    left_margin = 10mm,
-    bottom_margin = 10mm,
-    top_margin = 10mm,
-)
-Plots.png(plt3, joinpath(output_dir, resolution * "_" * "comparison"))
-Plots.pdf(plt3, joinpath(output_dir, resolution * "_" * "comparison"))
-
-
-Plots.GRBackend()
-plt4 = Plots.plot(
-    [log2.(nprocs_clima_atmos) log2.(nprocs_tempest)],
-    [scaling_efficiency_clima_atmos scaling_efficiency_tempest],
-    markershape = [:circle :square],
-    markercolor = [:blue :orange],
-    xticks = (
-        log2.(nprocs_clima_atmos),
-        [string(i) for i in nprocs_clima_atmos],
-    ),
-    xlabel = "# of MPI processes",
-    ylabel = "Efficiency (T1/N)/TN",
-    title = "Scaling efficiency (T_int = $t_int)",
-    label = ["ClimaAtmos (Float64)" "Tempest (Float64)"],
-    ylims = (0.0, Inf),
-    legend = :bottom,
-    grid = :true,
-    left_margin = 10mm,
-    bottom_margin = 10mm,
-    top_margin = 10mm,
-    annotations = [
-        (
-            log2.(nprocs_clima_atmos),
-            scaling_efficiency_clima_atmos .- 5,
-            [string(i) * "%" for i in scaling_efficiency_clima_atmos],
-            4,
-        ),
-        (
-            log2.(nprocs_tempest),
-            scaling_efficiency_tempest .- 5,
-            [string(i) * "%" for i in scaling_efficiency_tempest],
-            4,
-        ),
-    ],
-)
-Plots.png(plt4, joinpath(output_dir, resolution * "_" * "Scaling_efficiency"))
-Plots.pdf(plt4, joinpath(output_dir, resolution * "_" * "Scaling_efficiency"))
diff --git a/post_processing/contours_and_profiles.jl b/post_processing/contours_and_profiles.jl
deleted file mode 100644
index 7b01ddec4a8..00000000000
--- a/post_processing/contours_and_profiles.jl
+++ /dev/null
@@ -1,550 +0,0 @@
-import ClimaAtmos: time_from_filename
-import ClimaCore: Geometry, Spaces, Fields, InputOutput
-import ClimaComms
-import CairoMakie: Makie
-import Statistics: mean
-
-function best_time_unit(value_in_seconds)
-    unit_options = (
-        (60 * 60 * 24 * 365, "yr"),
-        (60 * 60 * 24, "d"),
-        (60 * 60, "h"),
-        (60, "min"),
-        (1, "s"),
-    )
-    minimum_value_in_units = 2
-    index = findfirst(unit_options) do (seconds_per_unit, _)
-        return value_in_seconds / seconds_per_unit >= minimum_value_in_units
-    end
-    return unit_options[isnothing(index) ? end : index]
-end
-
-function time_string_with_unit(value_in_seconds)
-    seconds_per_unit, abbreviation = best_time_unit(value_in_seconds)
-    value_in_units = round(value_in_seconds / seconds_per_unit; sigdigits = 6)
-    return "$value_in_units $abbreviation"
-end
-
-# This is an extended version of the time tick locator and formatter from the
-# Makie docs: https://docs.makie.org/stable/examples/blocks/axis/index.html.
-struct TimeTicks end
-function Makie.get_ticks(::TimeTicks, scale, formatter, vmin, vmax)
-    seconds_per_unit, abbreviation = best_time_unit(max(abs(vmin), abs(vmax)))
-    values_in_units = Makie.get_tickvalues(
-        Makie.automatic,
-        scale,
-        vmin / seconds_per_unit,
-        vmax / seconds_per_unit,
-    )
-    values_in_seconds = values_in_units .* seconds_per_unit
-    labels =
-        Makie.get_ticklabels(formatter, values_in_units) .* " $abbreviation"
-    return values_in_seconds, labels
-end
-
-function contour_plot!(
-    grid_position,
-    var_col_time_series,
-    ref_var_col_time_series,
-    category,
-    name,
-    draw_ref_plot,
-    label_ts,
-    label_zs,
-    negative_values_allowed,
-)
-    # There seems to be a bug in contourf! that occasionally causes the highest
-    # or lowest contour band to disappear. Adding a tiny amount of padding to
-    # the limits fixes this problem.
-    padding_fraction = 1e-6 # Increase this if any contour bands aren't plotted.
-
-    # There is a known issue in Makie that causes aliasing artifacts to appear
-    # between adjacent bands drawn by contourf! (and apparently also Colorbar):
-    # https://discourse.julialang.org/t/aliasing-with-makie-contourf/71783.
-    # Redrawing plots on top of themselves fixes this problem.
-    n_redraws = 2 # Increase this if any aliasing artifacts are observed.
-
-    if isnothing(var_col_time_series)
-        return
-    end
-    ts = getindex.(var_col_time_series, 1)
-    zs = var_col_time_series[1][2]
-
-    #values_adjoint = hcat(getindex.(var_col_time_series, 3)...)
-    #values = permutedims(values_adjoint, (2, 1))
-    values = hcat(getindex.(var_col_time_series, 3)...)'
-
-    if isnothing(ref_var_col_time_series)
-        limits = extrema(values)
-    else
-        ref_ts = getindex.(ref_var_col_time_series, 1)
-        ref_zs = ref_var_col_time_series[1][2]
-
-        #ref_values_adjoint = hcat(getindex.(ref_var_col_time_series, 3)...)
-        #ref_values = permutedims(ref_values_adjoint, (2, 1))
-        ref_values = hcat(getindex.(ref_var_col_time_series, 3)...)'
-
-        limits = extrema((extrema(values)..., extrema(ref_values)...))
-    end
-
-
-    limits = limits .+ (-1, 1) .* (padding_fraction * (limits[2] - limits[1]))
-    if negative_values_allowed
-        extendlow = nothing
-    else
-        limits = max.(limits, 0)
-        extendlow = :red
-    end
-
-    if limits[1] == limits[2]
-        colorbar_ticks_kwarg = (; ticks = [limits[1]])
-
-        # Since contourf! requires nonempty contours, the limits need to be
-        # expanded. If negative values are allowed, expand the limits so that
-        # the contour plot is filled with the color located in the middle of the
-        # color bar; otherwise, expand them so that it is filled with the color
-        # at the bottom of the color bar.
-        limits = limits[1] .+ (negative_values_allowed ? (-1, 1) : (0, 1))
-    else
-        colorbar_ticks_kwarg = (;)
-    end
-    n_bands = 20
-    levels = range(limits..., n_bands + 1)
-    kwargs = (; levels, extendlow, colormap = :viridis)
-
-    grid_layout = Makie.GridLayout(; default_rowgap = 10, default_colgap = 10)
-    grid_position[:, :] = grid_layout
-
-    axis = Makie.Axis(
-        grid_layout[1, :];
-        title = "$category $name",
-        xticks = TimeTicks(),
-        xlabel = label_ts && !draw_ref_plot ? "Time" : "",
-        ylabel = label_zs ?
-                 (draw_ref_plot ? "PR\n" : "") * "Elevation (km)" : "",
-        xticksvisible = !draw_ref_plot,
-        xticklabelsvisible = label_ts && !draw_ref_plot,
-        yticklabelsvisible = label_zs,
-    )
-    first_plot = Makie.contourf!(axis, ts, zs, values; kwargs...)
-    for _ in 1:n_redraws
-        Makie.contourf!(axis, ts, zs, values; kwargs...)
-    end
-
-    if draw_ref_plot
-        ref_axis = Makie.Axis(
-            grid_layout[2, :];
-            xticks = TimeTicks(),
-            xlabel = label_ts ? "Time" : "",
-            ylabel = label_zs ? "main\nElevation (km)" : "",
-            xticklabelsvisible = label_ts,
-            yticklabelsvisible = label_zs,
-        )
-        if isnothing(ref_var_col_time_series)
-            Makie.contourf!(ref_axis, ts, zs, fill(NaN, size(values)))
-        else
-            for _ in 1:(n_redraws + 1)
-                Makie.contourf!(ref_axis, ref_ts, ref_zs, ref_values; kwargs...)
-            end
-        end
-    end
-
-    for _ in 1:(n_redraws + 1)
-        Makie.Colorbar(grid_layout[:, 2], first_plot; colorbar_ticks_kwarg...)
-    end
-end
-
-function profile_plot!(
-    grid_position,
-    final_var_cols,
-    ref_final_var_cols,
-    categories,
-    name,
-    label_zs,
-)
-    # Use maximally distinguishable colors that are neither too bright nor too
-    # dark (restrict their luminance range from [0, 100] to [30, 40]).
-    colors = Makie.ColorSchemes.distinguishable_colors(6; lchoices = 30:40)
-    color_indices = Dict(:gm => 1, :draft => 2, :env => 3)
-    ref_color_indices = Dict(:gm => 4, :draft => 5, :env => 6)
-    alpha = 0.6 # Use transparency to see when lines are on top of each other.
-
-    axis = Makie.Axis(
-        grid_position;
-        xlabel = "$name",
-        ylabel = label_zs ? "Elevation (km)" : "",
-        yticksvisible = label_zs,
-        yticklabelsvisible = label_zs,
-    )
-    any_refs_available = any(!isnothing, getindex.(ref_final_var_cols, 2))
-    for (category, (zs, values)) in zip(categories, final_var_cols)
-        isnothing(values) && continue
-        color = (colors[color_indices[category]], alpha)
-        label = (any_refs_available ? "PR " : "") * "$category"
-        Makie.lines!(axis, values, zs; color, label)
-    end
-    for (category, (zs, values)) in zip(categories, ref_final_var_cols)
-        isnothing(values) && continue
-        color = (colors[ref_color_indices[category]], alpha)
-        label = "main $category"
-        Makie.lines!(axis, values, zs; color, label, linestyle = :dash)
-    end
-end
-
-function contours_and_profiles(output_dir, ref_job_id = nothing)
-    ##
-    ## Reading in time series data:
-    ##
-
-    function sorted_hdf5_files(dir)
-        file_paths = readdir(dir; join = true)
-        filter!(endswith(".hdf5"), file_paths)
-        sort!(file_paths; by = time_from_filename)
-        return file_paths
-    end
-
-    function read_hdf5_file(file_path)
-        reader = InputOutput.HDF5Reader(
-            file_path,
-            ClimaComms.SingletonCommsContext(ClimaComms.CPUSingleThreaded()),
-        )
-        diagnostics = InputOutput.read_field(reader, "diagnostics")
-        close(reader)
-        return time_from_filename(file_path), diagnostics
-    end
-
-    time_series = map(read_hdf5_file, sorted_hdf5_files(output_dir))
-    t_end = time_series[end][1]
-
-    ref_time_series = nothing
-    if !isnothing(ref_job_id) && haskey(ENV, "BUILDKITE_COMMIT")
-        main_dir = "/central/scratch/esm/slurm-buildkite/climaatmos-main"
-        isdir(main_dir) ||
-            error("Unable to find $main_dir when running on Buildkite")
-        commit_dirs = readdir(main_dir; join = true)
-
-        # Only commits that increment the "ref_counter" currently generate hdf5
-        # files. The files are zipped for TC cases, but not for other CI jobs.
-        # TODO: Unify how reference data is stored and simplify this code.
-        has_hdf5_files(commit_dir) =
-            isdir(joinpath(commit_dir, ref_job_id)) &&
-            any(endswith(".hdf5"), readdir(joinpath(commit_dir, ref_job_id)))
-        zip_file(commit_dir) = joinpath(commit_dir, ref_job_id, "hdf5files.zip")
-        filter!(commit_dirs) do commit_dir
-            return has_hdf5_files(commit_dir) || isfile(zip_file(commit_dir))
-        end
-        if isempty(commit_dirs)
-            @warn "Unable to find reference data for $ref_job_id in $main_dir"
-        else
-            latest_commit_dir =
-                argmax(commit_dir -> stat(commit_dir).mtime, commit_dirs)
-            commit_id = basename(latest_commit_dir)
-            @info "Loading reference data from main branch (commit $commit_id)"
-            if has_hdf5_files(latest_commit_dir)
-                ref_files_dir = joinpath(latest_commit_dir, ref_job_id)
-            else
-                ref_files_dir = mktempdir(output_dir; prefix = "unzip_dir_")
-                run(`unzip $(zip_file(latest_commit_dir)) -d $ref_files_dir`)
-            end
-            ref_file_paths = sorted_hdf5_files(ref_files_dir)
-            _, last_index = findmin(ref_file_paths) do ref_file_path
-                return abs(t_end - time_from_filename(ref_file_path))
-            end # End the reference time series as close to t_end as possible.
-            ref_time_series = map(read_hdf5_file, ref_file_paths[1:last_index])
-            ref_t_end = ref_time_series[end][1]
-        end
-    end
-
-    ##
-    ## Extracting a column of data from each time series:
-    ##
-
-    # Assume that all variables are on the same horizontal space.
-    horizontal_space(time_series) =
-        Spaces.horizontal_space(axes(time_series[1][2].temperature))
-
-    is_on_sphere(time_series) =
-        eltype(Fields.coordinate_field(horizontal_space(time_series))) <:
-        Geometry.LatLongPoint
-
-    function column_view_time_series(time_series, column)
-        ((i, j), h) = column
-        is_extruded_field(object) =
-            axes(object) isa Spaces.ExtrudedFiniteDifferenceSpace
-        column_view(field) = Fields.column(field, i, j, h)
-        column_zs(object) =
-            is_extruded_field(object) ?
-            vec(parent(Fields.coordinate_field(column_view(object)).z)) / 1000 :
-            nothing
-        column_values(object) =
-            is_extruded_field(object) ? vec(parent(column_view(object))) :
-            nothing
-        col_time_series = map(time_series) do (t, diagnostics)
-            objects = Fields._values(diagnostics)
-            return t, map(column_zs, objects), map(column_values, objects)
-        end
-
-        # Assume that all variables have the same horizontal coordinates.
-        coords =
-            Fields.coordinate_field(column_view(time_series[1][2].temperature))
-
-        coord_strings = map(filter(!=(:z), propertynames(coords))) do symbol
-            # Add 0 to every horizontal coordinate value so that -0.0 gets
-            # printed without the unnecessary negative sign.
-            value = round(mean(getproperty(coords, symbol)); sigdigits = 6) + 0
-            return "$symbol = $value"
-        end
-        col_string = "column data from $(join(coord_strings, ", "))"
-
-        return col_time_series, col_string
-    end
-
-    get_column_1(time_series) =
-        isnothing(time_series) ? (nothing, nothing) :
-        column_view_time_series(time_series, ((1, 1), 1))
-
-    column_at_latitude_getter(latitude) =
-        time_series -> begin
-            isnothing(time_series) && return (nothing, nothing)
-            column = if is_on_sphere(time_series)
-                horz_space = horizontal_space(time_series)
-                horz_coords = Fields.coordinate_field(horz_space)
-                FT = eltype(eltype(horz_coords))
-                target_column_coord =
-                    Geometry.LatLongPoint(FT(latitude), FT(0))
-                distance_to_target(((i, j), h)) =
-                    Geometry.great_circle_distance(
-                        Spaces.column(horz_coords, i, j, h)[],
-                        target_column_coord,
-                        horz_space.global_geometry,
-                    )
-                argmin(distance_to_target, Spaces.all_nodes(horz_space))
-            else
-                ((1, 1), 1) # If the data is not on a sphere, extract column 1.
-            end
-            return column_view_time_series(time_series, column)
-        end
-
-    # TODO: Add averaging over multiple columns.
-    column_info = if is_on_sphere(time_series)
-        map((0, 30, 60, 90)) do latitude
-            return column_at_latitude_getter(latitude), "_from_$(latitude)N_0E"
-        end
-    else
-        ((get_column_1, ""),)
-    end
-
-    ##
-    ## Extracting data for a specific variable from each column time series:
-    ##
-
-    function variable_column_time_series(col_time_series, category, name)
-        isnothing(col_time_series) && return nothing
-        symbol = category === :gm ? name : Symbol(category, :_, name)
-
-        # TODO: Remove this workaround when the old EDMF model is deleted.
-        deprecated_symbols = Dict(
-            :env_temperature => :env_temperature,
-            :draft_temperature => :bulk_up_temperature,
-            :env_potential_temperature => :env_theta_liq_ice,
-            :env_w_velocity => :face_env_w,
-            :draft_w_velocity => :face_bulk_w,
-            :draft_buoyancy => :bulk_up_buoyancy,
-            :draft_q_liq => :bulk_up_q_liq,
-            :draft_q_ice => :bulk_up_q_ice,
-            :env_relative_humidity => :env_RH,
-            :draft_cloud_fraction => :bulk_up_cloud_fraction,
-            :draft_area => :bulk_up_area,
-            :env_tke => :env_TKE,
-        )
-        if !hasproperty(col_time_series[1][2], symbol) &&
-           symbol in keys(deprecated_symbols)
-            symbol = deprecated_symbols[symbol]
-        end
-
-        hasproperty(col_time_series[1][2], symbol) || return nothing
-        return map(col_time_series) do (t, zs, diagnostics)
-            return t, getproperty(zs, symbol), getproperty(diagnostics, symbol)
-        end
-    end
-
-    final_variable_columns(col_time_series, categories, name) =
-        map(categories) do category
-            var_col_time_series =
-                variable_column_time_series(col_time_series, category, name)
-            return isnothing(var_col_time_series) ? (nothing, nothing) :
-                   var_col_time_series[end][2:3]
-        end
-
-    has_moisture = hasproperty(time_series[1][2], :q_vap)
-    has_precipitation = hasproperty(time_series[1][2], :q_rai)
-    has_sgs = hasproperty(time_series[1][2], :draft_area)
-
-    negative_values_allowed_names =
-        (:u_velocity, :v_velocity, :w_velocity, :buoyancy, :specific_enthalpy)
-
-    env_or_gm = has_sgs ? :env : :gm
-    draft_or_gm = has_sgs ? :draft : :gm
-    contour_variables = (
-        (draft_or_gm, :temperature),
-        (:gm, :potential_temperature),
-        (:gm, :u_velocity),
-        (:gm, :v_velocity),
-        (:gm, :w_velocity),
-        (draft_or_gm, :buoyancy),
-        (draft_or_gm, :specific_enthalpy),
-    )
-    if has_moisture
-        contour_variables = (
-            contour_variables...,
-            (env_or_gm, :relative_humidity),
-            (draft_or_gm, :q_tot),
-            (draft_or_gm, :q_vap),
-            (draft_or_gm, :q_liq),
-            (draft_or_gm, :q_ice),
-        )
-    end
-    if has_precipitation
-        contour_variables =
-            (contour_variables..., (draft_or_gm, :q_rai), (draft_or_gm, :q_sno))
-    end
-    if has_sgs
-        if has_moisture
-            contour_variables = (
-                contour_variables...,
-                (:gm, :cloud_fraction),
-                (:draft, :w_velocity),
-            )
-        end
-        contour_variables =
-            (contour_variables..., (:draft, :area), (:env, :tke))
-    end
-
-    all_categories = has_sgs ? (:gm, :env, :draft) : (:gm,)
-    profile_variables = (
-        (all_categories, :temperature),
-        (all_categories, :potential_temperature),
-        ((:gm,), :u_velocity),
-        ((:gm,), :v_velocity),
-        (all_categories, :w_velocity),
-        (all_categories, :buoyancy),
-        (all_categories, :specific_enthalpy),
-        (all_categories, :density),
-    )
-    if has_moisture
-        profile_variables = (
-            profile_variables...,
-            (all_categories, :relative_humidity),
-            (all_categories, :q_tot),
-            (all_categories, :q_vap),
-            (all_categories, :q_liq),
-            (all_categories, :q_ice),
-        )
-    end
-    if has_precipitation
-        profile_variables = (
-            profile_variables...,
-            (all_categories, :q_rai),
-            (all_categories, :q_sno),
-        )
-    end
-    if has_sgs
-        if has_moisture
-            profile_variables =
-                (profile_variables..., ((:gm,), :cloud_fraction))
-        end
-        profile_variables = (
-            profile_variables...,
-            ((:draft,), :area),
-            ((:env,), :tke),
-            ((:env,), :mixing_length),
-        )
-    end
-
-    ##
-    ## Generating the contour and profile plots:
-    ##
-
-    # Organize the variables into squares, or something close to squares.
-    # If they are not squares, the contour plots should have more rows than
-    # columns, whereas the profile plots should have more columns than rows.
-    sqrt_factors(n) = minmax(round(Int, sqrt(n)), cld(n, round(Int, sqrt(n))))
-    n_contour_cols, n_contour_rows = sqrt_factors(length(contour_variables))
-    n_profile_rows, n_profile_cols = sqrt_factors(length(profile_variables))
-
-    # Plot the variables by starting at the bottom of the first column and
-    # moving upward, then going to the bottom of the second column and moving
-    # upward, and so on. This ensures that the first column and the bottom row
-    # are always filled, which means that axis ticks do not need to be drawn in
-    # the remaining rows and columns.
-    function row_and_col(index, n_rows)
-        col, row = divrem(index - 1, n_rows, RoundDown) .+ 1
-        return n_rows + 1 - row, col
-    end
-
-    contour_resolution = (700 * n_contour_cols, 400 * n_contour_rows)
-    profile_resolution = (400 * n_profile_cols, 400 * n_profile_rows)
-
-    for (get_column, filename_suffix) in column_info
-        col_time_series, col_string = get_column(time_series)
-        ref_col_time_series, ref_col_string = get_column(ref_time_series)
-
-        contour_title =
-            isnothing(ref_col_time_series) || ref_col_string == col_string ?
-            col_string : "PR $col_string\nmain branch $ref_col_string"
-
-        profile_title =
-            isnothing(ref_col_time_series) ||
-            (ref_col_string == col_string && ref_t_end == t_end) ?
-            "$col_string after $(time_string_with_unit(t_end))" :
-            "PR $col_string after $(time_string_with_unit(t_end))\nmain branch \
-             $ref_col_string after $(time_string_with_unit(ref_t_end))"
-
-        # Assume that all variables have the same vertical coordinate limits.
-        z_limits = extrema(col_time_series[1][2].temperature)
-
-        figure = Makie.Figure()
-        for (index, (category, name)) in enumerate(contour_variables)
-            row, col = row_and_col(index, n_contour_rows)
-            contour_plot!(
-                figure[row, col],
-                variable_column_time_series(col_time_series, category, name),
-                variable_column_time_series(
-                    ref_col_time_series,
-                    category,
-                    name,
-                ),
-                category,
-                name,
-                !isnothing(ref_col_time_series),
-                row == n_contour_rows,
-                col == 1,
-                name in negative_values_allowed_names,
-            )
-        end
-        all_axes = filter(object -> object isa Makie.Axis, figure.content)
-        foreach(axis -> Makie.limits!(axis, (0, t_end), z_limits), all_axes)
-        Makie.Label(figure[0, :], contour_title; font = :bold, fontsize = 20)
-        file_path = joinpath(output_dir, "contours$(filename_suffix).png")
-        Makie.save(file_path, figure; resolution = contour_resolution)
-
-        figure = Makie.Figure()
-        for (index, (categories, name)) in enumerate(profile_variables)
-            row, col = row_and_col(index, n_profile_rows)
-            profile_plot!(
-                figure[row, col],
-                final_variable_columns(col_time_series, categories, name),
-                final_variable_columns(ref_col_time_series, categories, name),
-                categories,
-                name,
-                col == 1,
-            )
-        end
-        all_axes = filter(object -> object isa Makie.Axis, figure.content)
-        foreach(axis -> Makie.ylims!(axis, z_limits), all_axes)
-        Makie.Legend(figure[0, :], all_axes[1]; orientation = :horizontal)
-        Makie.Label(figure[-1, :], profile_title; font = :bold, fontsize = 20)
-        file_path = joinpath(output_dir, "final_profiles$(filename_suffix).png")
-        Makie.save(file_path, figure; resolution = profile_resolution)
-    end
-end
diff --git a/post_processing/define_tc_quicklook_profiles.jl b/post_processing/define_tc_quicklook_profiles.jl
deleted file mode 100644
index e607babb806..00000000000
--- a/post_processing/define_tc_quicklook_profiles.jl
+++ /dev/null
@@ -1,367 +0,0 @@
-import ClimaCore: Fields, InputOutput, Geometry
-ENV["GKSwstype"] = "nul";
-import Plots
-import ClimaAtmos as CA
-import ClimaComms
-
-function getfun(D, sym::Symbol)
-    if !hasproperty(D, sym)
-        @warn "Property $(sym) not in Diagnostics FieldVector."
-        return nothing
-    else
-        return (Y, D) -> getproperty(D, sym)
-    end
-end
-getfun(D, p::Function) = p
-
-function plot_tc_profiles(folder; hdf5_filename, main_branch_data_path)
-    args =
-        (; tickfontsize = 13, guidefontsize = 16, legendfontsize = 10, lw = 3)
-
-    # initialize all plot panes
-
-#! format: off
-    vars = Dict()
-    # vars[title] = (; fn = func to get variable(s), pfx = label prefix(es))
-    vars["area fraction"] =     (;fn=(:bulk_up_area, :env_area), pfx = ("up", "en"))
-    vars["buoy"] =              (;fn=(:face_up1_buoyancy, :face_env_buoyancy), pfx = ("up", "en"))
-    vars["T"] =                 (;fn=(:bulk_up_temperature, :env_temperature), pfx = ("up", "en"))
-    vars["CF"] =                (;fn=(:bulk_up_cloud_fraction, :env_cloud_fraction), pfx = ("up", "env"))
-
-    vars["up w"]  =             (;fn=(Y, D) -> Geometry.WVector.(D.face_bulk_w), pfx="")
-    vars["en w"] =              (;fn=(Y, D) -> Geometry.WVector.(D.face_env_w), pfx="")
-    vars["gm u"] =              (;fn=(Y, D) -> Geometry.UVector.(Y.c.uₕ), pfx = "")
-    vars["gm v"] =              (;fn=(Y, D) -> Geometry.VVector.(Y.c.uₕ), pfx = "")
-
-    vars["up qt"] =             (;fn=:bulk_up_q_tot, pfx = "")
-    vars["up ql"] =             (;fn=:bulk_up_q_liq, pfx = "")
-    vars["up qi"] =             (;fn=:bulk_up_q_ice, pfx = "")
-    vars["en qt"] =             (;fn=:env_q_tot, pfx = "")
-    vars["en ql"] =             (;fn=:env_q_liq, pfx = "")
-    vars["en qi"] =             (;fn=:env_q_ice, pfx = "")
-    vars["gm theta"] =          (;fn=:potential_temperature, pfx = "")
-    vars["en RH"] =             (;fn=:env_RH, pfx = "")
-    vars["en TKE"] =            (;fn=:env_TKE, pfx = "")
-    vars["up filter flag 1"] =  (;fn=:bulk_up_filter_flag_1, pfx = "")
-    vars["up filter flag 2"] =  (;fn=:bulk_up_filter_flag_2, pfx = "")
-    vars["up filter flag 3"] =  (;fn=:bulk_up_filter_flag_3, pfx = "")
-    vars["up filter flag 4"] =  (;fn=:bulk_up_filter_flag_4, pfx = "")
-#! format: on
-
-    p_all = map(k -> Plots.plot(; title = k, args...), collect(keys(vars)))
-
-    function add_to_plots!(input_filename; data_source)
-        if !isfile(input_filename)
-            @info "Data file `$input_filename` not found for data source `$data_source`."
-            return
-        end
-
-        reader = InputOutput.HDF5Reader(
-            input_filename,
-            ClimaComms.SingletonCommsContext(ClimaComms.CPUSingleThreaded()),
-        )
-        Y = InputOutput.read_field(reader, "Y")
-        D = InputOutput.read_field(reader, "diagnostics")
-
-        zc = parent(Fields.coordinate_field(Y.c).z)[:]
-        zf = parent(Fields.coordinate_field(Y.f).z)[:]
-
-        for (i, (k, v)) in enumerate(vars)
-            if v.pfx isa Tuple
-                for j in 1:length(v.pfx)
-                    fn = getfun(D, v.fn[j])
-                    isnothing(fn) && continue
-                    var = fn(Y, D)
-                    z = parent(Fields.coordinate_field(axes(var)).z)[:]
-                    vardata = parent(var)[:]
-                    prefix = isempty(v.pfx[j]) ? "" : "$(v.pfx[j]) "
-                    Plots.plot!(
-                        p_all[i],
-                        vardata,
-                        z;
-                        label = "$prefix$data_source",
-                    )
-                end
-            elseif v.pfx isa String
-                fn = getfun(D, v.fn)
-                isnothing(fn) && continue
-                var = fn(Y, D)
-                z = parent(Fields.coordinate_field(axes(var)).z)[:]
-                vardata = parent(var)[:]
-                prefix = isempty(v.pfx) ? "" : "$(v.pfx) "
-                Plots.plot!(p_all[i], vardata, z; label = "$prefix$data_source")
-            end
-        end
-    end
-
-    PR_filename = joinpath(folder, hdf5_filename)
-    add_to_plots!(PR_filename; data_source = "PR")
-    if ispath(main_branch_data_path)
-        main_filename = joinpath(main_branch_data_path, hdf5_filename)
-        add_to_plots!(main_filename; data_source = "main")
-    end
-
-    n_cols = 5
-    more_args = (;
-        size = (2400.0, 1500.0),
-        bottom_margin = 20.0 * Plots.PlotMeasures.px,
-        left_margin = 20.0 * Plots.PlotMeasures.px,
-        layout = (ceil(Int, length(p_all) / n_cols), n_cols),
-    )
-    p = Plots.plot(p_all...; more_args...)
-
-    # Save output
-    output_filename = joinpath(
-        folder,
-        "____________________________________final_profiles.png",
-    )
-    Plots.png(p, output_filename)
-end
-
-"""
-    get_contours(vars, input_filenames; data_source, have_main)
-
-An array of contour plots, given
- - `vars::Tuple{String,Function}` a Tuple of variable name and function
-    that obtains the field given the diagnostics
- - `input_filenames::Array{String}` an array of input filenames
- - `data_source::String` a string containing the data source
- - `have_main::Bool` indicates whether the main branch data is available
-    in order to adjust the margins
-"""
-function get_contours(vars, input_filenames; data_source, have_main)
-    files_found = map(x -> isfile(x), input_filenames)
-    if !all(files_found)
-        @warn "Some data files were missing in data source `$data_source`."
-        @warn "$(count(files_found)) files found out of $(length(input_filenames))"
-        @warn "Filtering out only found files."
-        input_filenames = filter(isfile, input_filenames)
-    end
-
-    data = map(input_filenames) do input_filename
-        reader = InputOutput.HDF5Reader(
-            input_filename,
-            ClimaComms.SingletonCommsContext(ClimaComms.CPUSingleThreaded()),
-        )
-        Y = InputOutput.read_field(reader, "Y")
-        D = InputOutput.read_field(reader, "diagnostics")
-        (Y, D)
-    end
-    (Ys, Ds) = first.(data), last.(data)
-    t = CA.time_from_filename.(input_filenames)
-
-    clims = Dict()
-    contours = get_empty_contours(vars)
-    K = collect(keys(contours))
-    n = length(K)
-    fig_width = 4500
-    fig_height = 5000
-    left_side = data_source == "main"
-    if have_main
-        l_margin = left_side ? 100 : 0
-        r_margin = left_side ? 0 : 0 # else case doesn't seem to be working
-    else
-        l_margin = 100
-        r_margin = -100
-    end
-    parent_data(data) = parent(data)[:]
-
-    for (i, name) in enumerate(K)
-        fn = getfun(first(Ds), contours[name].fn)
-        isnothing(fn) && continue
-        space = axes(fn(first(Ys), first(Ds)))
-        z = parent(Fields.coordinate_field(space).z)[:]
-        Ds_parent = parent_data.(fn.(Ys, Ds))
-        cdata = hcat(Ds_parent...)
-        clims[name] = (minimum(cdata), maximum(cdata))
-        @info "clims[$name] ($data_source) = $(clims[name])"
-        Plots.contourf!(
-            contours[name].plot,
-            t ./ 3600,
-            z ./ 10^3,
-            cdata;
-            xticks = i == n,
-            yticks = !have_main || left_side,
-            colorbar = !have_main || !left_side,
-            c = :viridis,
-            xlabel = i == n ? "Time (hours)" : "",
-            ylabel = left_side || !have_main ? "Height (km)" : "",
-            left_margin = l_margin * Plots.PlotMeasures.px,
-            right_margin = r_margin * Plots.PlotMeasures.px,
-            bottom_margin = i == n ? 50 * Plots.PlotMeasures.px :
-                            0 * Plots.PlotMeasures.px,
-            top_margin = 0 * Plots.PlotMeasures.px,
-            size = (fig_width, fig_height), # extra space for colorbar
-            title = "$name ($data_source)",
-        )
-    end
-    return contours, clims
-end
-
-hdf5_files(path) = filter(x -> endswith(x, ".hdf5"), readdir(path, join = true))
-
-function zip_and_cleanup_output(path, zip_file)
-    files = basename.(hdf5_files(path))
-    cd(path) do
-        if !isfile(zip_file)
-            zipname = first(splitext(basename(zip_file)))
-            run(pipeline(Cmd(["zip", zipname, files...]), stdout = IOBuffer()))
-        end
-        for f in files
-            rm(f)
-        end
-    end
-end
-
-function unzip_file_in_path(path, zip_file, unzip_path)
-    if !ispath(path)
-        @warn "Path $path not found."
-        return nothing
-    end
-    if !isfile(joinpath(path, zip_file))
-        @warn "Zip file does not exist"
-        @show path
-        @show zip_file
-        @show readdir(path, join = true)
-        return nothing
-    end
-
-    cp(joinpath(path, zip_file), joinpath(unzip_path, zip_file))
-    @info "Unzipping files in `$path`"
-    cd(unzip_path) do
-        zipname = first(splitext(basename(zip_file)))
-        run(pipeline(Cmd(["unzip", zipname]); stdout = IOBuffer()))
-    end
-    files = hdf5_files(unzip_path)
-    @assert !isempty(files)
-end
-
-function get_main_filenames(main_branch_data_path)
-    if ispath(main_branch_data_path)
-        files = hdf5_files(main_branch_data_path)
-        if any(isfile.(files))
-            CA.sort_files_by_time(files)
-        else
-            nothing
-        end
-    else
-        nothing
-    end
-end
-
-function debug_filenames(filenames, data_source)
-    filenames == nothing && return
-    println("---------------- Files for $data_source")
-    for f in filenames
-        println("    File: $f")
-    end
-    println("---- Additional info")
-    @show allunique(filenames)
-    @show allunique(hdf5_files(dirname(first(filenames))))
-    println("----------------")
-end
-
-function plot_tc_contours(folder; main_branch_data_path)
-    PR_filenames = CA.sort_files_by_time(hdf5_files(folder))
-    main_filenames = get_main_filenames(main_branch_data_path)
-    # debug_filenames(PR_filenames, "PR")
-    # debug_filenames(main_filenames, "main")
-    _plot_tc_contours(folder; PR_filenames, main_filenames)
-end
-
-function get_empty_contours(vars)
-    plots = map(vars) do name_fn
-        name = first(name_fn)
-        fn = last(name_fn)
-        cplot = Plots.contourf()
-        Pair(name, (; plot = cplot, fn))
-    end
-    plots = Dict(plots...)
-    return plots
-end
-
-union_clims(a::Tuple{T, T}, b::Tuple{T, T}) where {T} =
-    (min(first(a), first(b)), max(last(a), last(b)))
-
-is_trivial_clims(tup::Tuple{T, T}) where {T} = tup[1] == tup[2]
-
-function union_clims(a::Dict, b::Dict)
-    clims_dict = Dict()
-    for k in union(keys(a), keys(b))
-        clims_tup = if haskey(a, k) && haskey(b, k)
-            union_clims(a[k], b[k])
-        elseif haskey(a, k)
-            a[k]
-        else
-            b[k]
-        end
-        # Avoid https://github.com/JuliaPlots/Plots.jl/issues/3924
-        clims_dict[k] =
-            is_trivial_clims(clims_tup) ? NamedTuple() : (; clims = clims_tup)
-        @info "union_clims[$k] = $(clims_dict[k])"
-    end
-    return clims_dict
-end
-
-function _plot_tc_contours(folder; PR_filenames, main_filenames)
-
-    vars = [
-        ("area fraction", :bulk_up_area),
-        ("up qt", :bulk_up_q_tot),
-        ("up ql", :bulk_up_q_liq),
-        ("up qi", :bulk_up_q_ice),
-        ("en qt", :env_q_tot),
-        ("en TKE", :env_TKE),
-        ("gm theta", :potential_temperature),
-        ("up filter flag 1", :bulk_up_filter_flag_1),
-        ("up filter flag 2", :bulk_up_filter_flag_2),
-        ("up filter flag 3", :bulk_up_filter_flag_3),
-        ("up filter flag 4", :bulk_up_filter_flag_4),
-        # ("up qr", D-> parent(D.)[:])
-    ]
-
-    have_main = main_filenames ≠ nothing
-    contours_PR, clims_PR =
-        get_contours(vars, PR_filenames; data_source = "PR", have_main)
-    if have_main
-        contours_main, clims_main =
-            get_contours(vars, main_filenames; data_source = "main", have_main)
-        clims = union_clims(clims_main, clims_PR)
-        for k in keys(contours_main)
-            haskey(clims, k) || continue # need to skip when adding new variables
-            Plots.contourf!(contours_main[k].plot; clims[k]...)
-        end
-        for k in keys(contours_PR)
-            haskey(clims, k) || continue # need to skip when adding new variables
-            Plots.contourf!(contours_PR[k].plot; clims[k]...)
-        end
-        # Get array of plots from dict of NamedTuples
-        contours_main =
-            map(k -> contours_main[k].plot, collect(keys(contours_main)))
-        contours_PR = map(k -> contours_PR[k].plot, collect(keys(contours_PR)))
-        P = Iterators.flatten(collect(zip(contours_main, contours_PR)))
-        p = Plots.plot(
-            P...;
-            layout = Plots.grid(
-                length(contours_main),
-                2;
-                widths = [0.45, 0.55],
-            ),
-        )
-    else
-        @warn "No main branch to compare against"
-        # Get array of plots from dict of NamedTuples
-        contours_PR = map(k -> contours_PR[k].plot, collect(keys(contours_PR)))
-        p = Plots.plot(
-            contours_PR...;
-            layout = Plots.grid(length(contours_PR), 1),
-        )
-    end
-
-    # Save output
-    output_filename = joinpath(
-        folder,
-        "____________________________________final_contours.png",
-    )
-    Plots.png(p, output_filename)
-end
diff --git a/post_processing/make_plots.jl b/post_processing/make_plots.jl
new file mode 100644
index 00000000000..cbcfd9f08ae
--- /dev/null
+++ b/post_processing/make_plots.jl
@@ -0,0 +1,93 @@
+import CairoMakie
+import ClimaAnalysis
+import ClimaAnalysis: Visualize as viz
+
+# swap_2_z plots z as second dimension and computes contours on the third
+function dim1_dim2_z_contour_plot!(
+    fig,
+    simdir,
+    p_loc;
+    short_name,
+    reduction = "average",
+    period = "1d",
+    swap_2_z = false,
+)
+    var = get(simdir; short_name, reduction, period)
+
+    t_name, dim1_name, dim2_name, z_name = var.index2dim
+    if swap_2_z
+        dim2_name, z_name = z_name, dim2_name
+    end
+    time = var.dims[t_name][end]
+    dim1 = var.dims[dim1_name]
+    dim2 = var.dims[dim2_name]
+    ZIDX = 3
+    z = var.dims[z_name][ZIDX]
+
+    Z = swap_2_z ? var.var[end, :, ZIDX, :] : var.var[end, :, :, ZIDX]
+
+    units = var.attributes["units"]
+    dim1_units = var.dim_attributes[dim1_name]["units"]
+    dim2_units = var.dim_attributes[dim2_name]["units"]
+    z_units = var.dim_attributes[z_name]["units"]
+    t_units = var.dim_attributes[t_name]["units"]
+
+    title =
+        var.attributes["long_name"] *
+        " $z_name = $z $z_units, $t_name = $time $t_units"
+    viz.contour_plot!(
+        fig;
+        X = dim1,
+        Y = dim2,
+        Z = Z,
+        title,
+        xlabel = "$dim1_name [$dim1_units]",
+        ylabel = "$dim2_name [$dim2_units]",
+        colorbar_label = "$short_name [$units]",
+        p_loc = Tuple(p_loc),
+    )
+end
+
+function make_plots(::Val{:sphere_baroclinic_wave_rhoe}, simulation_path)
+    simdir = ClimaAnalysis.SimDir(simulation_path)
+
+    short_names = ["pfull", "ta"]
+
+    fig = CairoMakie.Figure(resolution = (300 * length(short_names), 600))
+    p_loc = [1, 1]
+
+    for short_name in short_names
+        dim1_dim2_z_contour_plot!(fig, simdir, p_loc; short_name)
+        p_loc[1] = p_loc[1] + 1
+    end
+
+    CairoMakie.save(joinpath(simulation_path, "summary.png"), fig)
+end
+
+function make_plots(::Val{:bomex_box_rhoe}, simulation_path)
+    simdir = ClimaAnalysis.SimDir(simulation_path)
+
+    short_names = ["wa", "hur"]
+
+    fig = CairoMakie.Figure(resolution = (300 * length(short_names), 600))
+    p_loc = [1, 1]
+
+    for short_name in short_names
+        dim1_dim2_z_contour_plot!(
+            fig,
+            simdir,
+            p_loc;
+            short_name,
+            reduction = "inst",
+            period = nothing,
+            swap_2_z = true,
+        )
+        p_loc[1] = p_loc[1] + 1
+    end
+
+    CairoMakie.save(joinpath(simulation_path, "summary.png"), fig)
+end
+
+function make_plots(sim, simulation_path)
+    @warn "No plots for $sim"
+end
diff --git a/post_processing/plot/plot_helpers.jl b/post_processing/plot/plot_helpers.jl
deleted file mode 100644
index 1732afd7cd5..00000000000
--- a/post_processing/plot/plot_helpers.jl
+++ /dev/null
@@ -1,710 +0,0 @@
-using NCDatasets
-import ClimaCoreSpectra: power_spectrum_1d, power_spectrum_2d
-using Statistics
-
-"""
-    create_plot!(fig; X,Y,Z, args...)
-
-Create and save a Makie plot object, which consists 
-of either a line-plot (X,Y) or a surface plot ((X,Y),Z),
-Basic functionality is included, this can be extended with
-`args` compatible with the Makie plot package.  
-"""
-function create_plot!(
-    fig::Figure;
-    X = lon,
-    Y = lat,
-    Z = nothing,
-    p_loc::Tuple = (1, 1),
-    title = "",
-    xlabel = "Longitude",
-    ylabel = "Latitude",
-    xscale = identity,
-    yscale = identity,
-    linewidth = 6,
-    level::Int = 1,
-    timeind::Int = 1,
-)
-    if Z == nothing
-        generic_axis = fig[p_loc[1], p_loc[2]] = GridLayout()
-        Axis(generic_axis[1, 1]; title, xlabel, ylabel, xscale, yscale)
-        CairoMakie.lines!(X, Y; title, linewidth)
-    else
-        generic_axis = fig[p_loc[1], p_loc[2]] = GridLayout() # Generic Axis Layout
-        Axis(generic_axis[1, 1]; title, xlabel, ylabel, yscale) # Plot specific attributes
-        # custom_levels is a workaround for plotting constant fields with CairoMakie
-        custom_levels =
-            minimum(Z) ≈ maximum(Z) ? (minimum(Z):0.1:(minimum(Z) + 0.2)) : 25
-        generic_plot = CairoMakie.contourf!(X, Y, Z, levels = custom_levels) # Add plot contents
-        Colorbar(generic_axis[1, 2], generic_plot)
-    end
-end
-
-"""
-    generate_empty_figure(;resolution, bgcolor, fontsize)
-
-Generate an empty `Makie.Figure` object, with specified
-`resolution`, background color `bgcolor` and `fontsize`.
-Returns the empty figure which can then be populated with 
-the `create_plot!(...)`. 
-"""
-function generate_empty_figure(;
-    resolution::Tuple = (4098, 2048),
-    bgcolor::Tuple = (0.98, 0.98, 0.98),
-    fontsize = 36,
-)
-    fig = Figure(;
-        backgroundcolor = RGBf(bgcolor[1], bgcolor[2], bgcolor[3]),
-        resolution,
-        fontsize,
-    )
-    return fig
-end
-
-function generate_paperplots_dry_baro_wave(fig_dir, nc_files)
-    for day in [8, 10, 100]
-        ncfile = filter(x -> endswith(x, "day$day.0.nc"), nc_files)
-        if !isempty(ncfile)
-            nt = NCDataset(ncfile[1], "r") do nc
-                lon = Array(nc["lon"])
-                lat = Array(nc["lat"])
-                z = Array(nc["z"])
-                p = Array(nc["pressure"])
-                T = Array(nc["temperature"])
-                vort = Array(nc["vorticity"])
-                (; lon, lat, z, p, T, vort)
-            end
-            (; lon, lat, z, p, T, vort) = nt
-
-            latidx = findall(x -> x >= 0, lat)
-            lonidx = findall(x -> 0 <= x <= 240, lon)
-
-            FT = eltype(vort)
-
-            fig = generate_empty_figure()
-            create_plot!(
-                fig;
-                X = lon[lonidx],
-                Y = lat[latidx],
-                Z = p[lonidx, latidx, 1, 1],
-                title = "Pressure: Day $day; z $(round(z[1])) m",
-            )
-            create_plot!(
-                fig;
-                p_loc = (2, 1),
-                X = lon[lonidx],
-                Y = lat[latidx],
-                Z = T[lonidx, latidx, 1, 1],
-                title = "Temperature: Day $day; z $(round(z[1])) m",
-            )
-            create_plot!(
-                fig;
-                p_loc = (3, 1),
-                X = lon[lonidx],
-                Y = lat[latidx],
-                Z = vort[lonidx, latidx, 1, 1] .* FT(1e5),
-                title = "Vorticity(× 10⁵): Day $day; z $(round(z[1])) m",
-            )
-            CairoMakie.save(fig_dir * "/dbw_day$day.png", fig)
-        else
-            @warn "day$day.0.nc DOES NOT EXIST!!!"
-        end
-    end
-end
-
-function generate_paperplots_moist_baro_wave(fig_dir, nc_files)
-    for day in [8, 10, 100]
-        ncfile = filter(x -> endswith(x, "day$day.0.nc"), nc_files)
-        if !isempty(ncfile)
-            nt = NCDataset(ncfile[1], "r") do nc
-                lon = Array(nc["lon"])
-                lat = Array(nc["lat"])
-                z = Array(nc["z"])
-                p = Array(nc["pressure"])
-                T = Array(nc["temperature"])
-                vort = Array(nc["vorticity"])
-                qt = Array(nc["qt"])
-                cloud_water =
-                    Array(nc["cloud_liquid"]) + Array(nc["cloud_ice"])
-                water_vapor = Array(nc["water_vapor"])
-                rho = Array(nc["rho"])
-                w = Array(nc["w"])
-                u = Array(nc["u"])
-                v = Array(nc["v"])
-                z_sfc = Array(nc["sfc_elevation"])
-                (;
-                    lon,
-                    lat,
-                    z,
-                    p,
-                    T,
-                    vort,
-                    qt,
-                    cloud_water,
-                    water_vapor,
-                    rho,
-                    w,
-                    u,
-                    v,
-                    z_sfc,
-                )
-            end
-            (;
-                lon,
-                lat,
-                z,
-                p,
-                T,
-                vort,
-                qt,
-                cloud_water,
-                water_vapor,
-                rho,
-                w,
-                u,
-                v,
-                z_sfc,
-            ) = nt
-
-            vert_intg_cloud_water =
-                sum(cloud_water .* rho, dims = 3) ./ sum(rho, dims = 3)
-            vert_intg_water_vapor =
-                sum(water_vapor .* rho, dims = 3) ./ sum(rho, dims = 3)
-
-            latidx = findall(x -> x >= 0, lat)
-            lonidx = findall(x -> 0 <= x <= 240, lon)
-
-            FT = eltype(vort)
-
-            # Basic
-            fig = generate_empty_figure()
-
-            create_plot!(
-                fig;
-                p_loc = (1, 1),
-                X = lon[lonidx],
-                Y = lat[latidx],
-                Z = p[lonidx, latidx, 1, 1],
-                title = "Pressure: Day $day; z $(round(z[1])) m",
-            )
-            create_plot!(
-                fig;
-                p_loc = (2, 1),
-                X = lon[lonidx],
-                Y = lat[latidx],
-                Z = T[lonidx, latidx, 1, 1],
-                title = "Temperature: Day $day; z $(round(z[1])) m",
-            )
-            create_plot!(
-                fig;
-                p_loc = (3, 1),
-                X = lon[lonidx],
-                Y = lat[latidx],
-                Z = vort[lonidx, latidx, 1, 1] .* FT(1e5),
-                title = "Vorticity(×10⁵): Day $day; z $(round(z[1])) m",
-            )
-            CairoMakie.save(fig_dir * "/mbw_day$day.png", fig)
-
-            fig = generate_empty_figure()
-            create_plot!(
-                fig;
-                p_loc = (1, 1),
-                X = lon[lonidx],
-                Y = lat[latidx],
-                Z = w[lonidx, latidx, 1, 1],
-                title = "Vertical Velocity: Day $day; z $(round(z[1])) m",
-            )
-            create_plot!(
-                fig;
-                p_loc = (2, 1),
-                X = lon[lonidx],
-                Y = lat[latidx],
-                Z = w[lonidx, latidx, 3, 1],
-                title = "Vertical Velocity: Day $day; z $(round(z[3])) m",
-            )
-            create_plot!(
-                fig;
-                p_loc = (1, 2),
-                X = lon[lonidx],
-                Y = lat[latidx],
-                Z = u[lonidx, latidx, 1, 1],
-                title = "Zonal Velocity: Day $day; z $(round(z[1])) m",
-            )
-            create_plot!(
-                fig;
-                p_loc = (2, 2),
-                X = lon[lonidx],
-                Y = lat[latidx],
-                Z = u[lonidx, latidx, 3, 1],
-                title = "Zonal Velocity: Day $day; z $(round(z[3])) m",
-            )
-            create_plot!(
-                fig;
-                p_loc = (1, 3),
-                X = lon[lonidx],
-                Y = lat[latidx],
-                Z = v[lonidx, latidx, 1, 1],
-                title = "Meridional Velocity: Day $day; z $(round(z[1])) m",
-            )
-            create_plot!(
-                fig;
-                p_loc = (2, 3),
-                X = lon[lonidx],
-                Y = lat[latidx],
-                Z = v[lonidx, latidx, 3, 1],
-                title = "Meridional Velocity: Day $day; z $(round(z[3])) m",
-            )
-            CairoMakie.save(fig_dir * "/mbw_velocity_day$day.png", fig)
-
-            # Moisture 
-            fig = generate_empty_figure()
-            create_plot!(
-                fig;
-                p_loc = (1, 1),
-                X = lon[lonidx],
-                Y = lat[latidx],
-                Z = qt[lonidx, latidx, 1, 1] .* 1e3,
-                title = "Tot. Specific Moisture (qt×10³): Day $day; z $(round(z[1])) m",
-            )
-            create_plot!(
-                fig;
-                p_loc = (2, 1),
-                X = lon[lonidx],
-                Y = lat[latidx],
-                Z = vert_intg_cloud_water[lonidx, latidx, 1, 1] .* 1e3,
-                title = "Vert Int Cloud Water [g/kg]: Day $day; z $(round(z[1])) m",
-            )
-            create_plot!(
-                fig;
-                p_loc = (3, 1),
-                X = lon[lonidx],
-                Y = lat[latidx],
-                Z = vert_intg_water_vapor[lonidx, latidx, 1, 1] .* 1e3,
-                title = "Vert Int Cloud Vapor [g/kg]: Day $day; z $(round(z[1])) m",
-            )
-            CairoMakie.save(fig_dir * "/mbw_moisture_day$day.png", fig)
-
-            # Spectrum calculation
-            mass_weight = ones(FT, 1) # only 1 level
-
-            # Compute 2d spectrum of u-component
-            u_2d_spectrum, wave_numbers, spherical, mesh_info =
-                power_spectrum_2d(FT, u[:, :, 1, 1], mass_weight) # use first level, for default CI config
-            # size(u) = (180, 90, 10, 1) (10 vert levels for CI default config, 1 component)
-
-            # Compute 2d spectrum of v-component
-            v_2d_spectrum, wave_numbers, spherical, mesh_info =
-                power_spectrum_2d(FT, v[:, :, 1, 1], mass_weight) # use first level
-
-            orography_spectrum, wave_numbers, spherical, mesh_info =
-                power_spectrum_2d(FT, z_sfc[:, :, 1, 1], mass_weight) # use first level, for default CI config
-
-            spectrum_2d = 0.5 .* (u_2d_spectrum + v_2d_spectrum)
-
-            # Spectra
-            fig = generate_empty_figure()
-            create_plot!(
-                fig;
-                p_loc = (1, 1),
-                X = collect(0:1:(mesh_info.num_fourier))[:],
-                Y = collect(0:1:(mesh_info.num_spherical))[:],
-                Z = (spectrum_2d[:, :, 1]),
-                xlabel = "",
-                ylabel = "",
-                title = "2D KE (horz) Spectrum: Day $day; z $(round(z[1])) m",
-            )
-            create_plot!(
-                fig;
-                p_loc = (2, 1),
-                X = collect(0:1:(mesh_info.num_fourier))[2:(end - 1)], # plot against the zonal wavenumber, m
-                Y = sum(spectrum_2d[:, :, 1], dims = 2)[2:(end - 1), 1], # sum along the total wavenumber, n
-                xlabel = "",
-                ylabel = "",
-                title = "Σₙ KE (horz) Spectrum: Day $day; z $(round(z[1])) m",
-            )
-            create_plot!(
-                fig;
-                p_loc = (3, 1),
-                X = collect(0:1:(mesh_info.num_spherical))[2:(end - 1)], # plot against the total wavenumber, n
-                Y = log.(sum(spectrum_2d[:, :, 1], dims = 1))'[2:(end - 1), 1], # sum along the zonal wavenumber, m
-                xlabel = "",
-                ylabel = "",
-                title = "Zonal Wavenumber (m) KE (horz) Spectrum: Day $day; z $(round(z[1])) m",
-            )
-            CairoMakie.save(fig_dir * "/mbw_spectrum_day$day.png", fig)
-        else
-            @warn "day$day.0.nc DOES NOT EXIST!!!"
-        end
-    end
-end
-
-function generate_elevation_spectra(fig_dir, nc_files)
-    ncfile = filter(x -> endswith(x, "day0.0.nc"), nc_files)
-    if !isempty(ncfile)
-        nt = NCDataset(ncfile[1], "r") do nc
-            z_sfc = Array(nc["sfc_elevation"])
-            lon = Array(nc["lon"])
-            lat = Array(nc["lat"])
-            (; lon, lat, z_sfc)
-        end
-        (; lon, lat, z_sfc) = nt
-
-        FT = eltype(lat)
-        # Spectrum calculation
-        mass_weight = ones(FT, 1) # only 1 level
-
-        orography_spectrum, wave_numbers, spherical, mesh_info =
-            power_spectrum_2d(FT, z_sfc[:, :, 1, 1], mass_weight) # use first level, for default CI config
-
-        fig = generate_empty_figure()
-        create_plot!(
-            fig;
-            p_loc = (1, 1),
-            X = collect(0:1:(mesh_info.num_fourier))[2:end], # plot against the zonal wavenumber, m
-            Y = log.(sum(orography_spectrum[:, :, 1], dims = 2))[2:end], # sum along the total wavenumber, n
-            title = "Diagnostic: Surface Elevation Spectrum",
-            xlabel = "Zonal wavenumber",
-            ylabel = "Log surface elevation spectrum",
-        )
-
-        CairoMakie.save(fig_dir * "/surface_elev_spectrum.png", fig)
-
-    else
-        @warn "day0.0.nc DOES NOT EXIST!!!"
-    end
-end
-
-# plots in the held-suarez paper: https://journals.ametsoc.org/view/journals/bams/75/10/1520-0477_1994_075_1825_apftio_2_0_co_2.xml?tab_body=pdf
-calc_zonalave_timeave(x) =
-    dropdims(dropdims(mean(mean(x, dims = 1), dims = 4), dims = 4), dims = 1)
-calc_zonalave_variance(x) = calc_zonalave_timeave((x .- mean(x, dims = 1)) .^ 2)
-calc_zonalave_covariance(x, y) =
-    calc_zonalave_timeave((x .- mean(x, dims = 1)) .* (y .- mean(y, dims = 1)))
-
-function generate_paperplots_held_suarez(fig_dir, nc_files; moist)
-    days = map(x -> parse(Int, split(basename(x), ".")[1][4:end]), nc_files)
-    # filter days for plotting:
-    # if simulated time is more than 200 days, days after day 200 will be collected and used for longterm average
-    # if simulated time is less than 200 days, the last day will be used to create a napshot
-    if maximum(days) > 200
-        filter!(x -> x > 200, days)
-    else
-        filter!(x -> x == maximum(days), days)
-    end
-
-    # identify nc files for plotting 
-    nc_file = []
-    for day in unique(days)
-        push!(nc_file, filter(x -> occursin(string(day), x), nc_files)...)
-    end
-
-    println("files to generate the plots:")
-    println.(nc_file)
-
-    # Fluxes
-    sfc_flux_energy = sfc_evaporation = lw_flux_net = sw_flux_net = nothing
-
-    # collect data from all nc files and combine data from different time step into one array
-    for (i, ifile) in enumerate(nc_file)
-        nc = NCDataset(ifile, "r")
-        if i == 1
-            global lat = Array(nc["lat"])
-            global lon = Array(nc["lon"])
-            global z = Array(nc["z"])
-            global u = Array(nc["u"])
-            global v = Array(nc["v"])
-            global w = Array(nc["w"])
-            global T = Array(nc["temperature"])
-            global θ = Array(nc["potential_temperature"])
-            if moist
-                global qt = Array(nc["qt"])
-            end
-            if haskey(nc, "sfc_flux_energy")
-                sfc_flux_energy = Array(nc["sfc_flux_energy"])
-            end
-            if haskey(nc, "sfc_evaporation")
-                sfc_evaporation = Array(nc["sfc_evaporation"])
-            end
-            if haskey(nc, "lw_flux_up")
-                # NOTE: We define net as up - down
-                lw_flux_net =
-                    Array(nc["lw_flux_up"]) - Array(nc["lw_flux_down"])
-                sw_flux_net =
-                    Array(nc["sw_flux_up"]) - Array(nc["sw_flux_down"])
-            end
-        else
-            u = cat(u, Array(nc["u"]), dims = 4)
-            v = cat(v, Array(nc["v"]), dims = 4)
-            w = cat(w, Array(nc["w"]), dims = 4)
-            T = cat(T, Array(nc["temperature"]), dims = 4)
-            θ = cat(θ, Array(nc["potential_temperature"]), dims = 4)
-            if moist
-                qt = cat(qt, Array(nc["qt"]), dims = 4)
-            end
-            if haskey(nc, "sfc_flux_energy")
-                sfc_flux_energy =
-                    cat(sfc_flux_energy, Array(nc["sfc_flux_energy"]), dims = 3)
-            end
-            if haskey(nc, "sfc_evaporation")
-                sfc_evaporation =
-                    cat(sfc_evaporation, Array(nc["sfc_evaporation"]), dims = 3)
-            end
-            if haskey(nc, "lw_flux_up")
-                lw_flux_net = cat(
-                    lw_flux_net,
-                    Array(nc["lw_flux_up"]) - Array(nc["lw_flux_down"]),
-                    dims = 4,
-                )
-                sw_flux_net = cat(
-                    sw_flux_net,
-                    Array(nc["sw_flux_up"]) - Array(nc["sw_flux_down"]),
-                    dims = 4,
-                )
-            end
-        end
-    end
-
-    # compute longterm zonal mean statistics
-    u_timeave_zonalave = calc_zonalave_timeave(u)
-    v_timeave_zonalave = calc_zonalave_timeave(v)
-    w_timeave_zonalave = calc_zonalave_timeave(w)
-    T_timeave_zonalave = calc_zonalave_timeave(T)
-    θ_timeave_zonalave = calc_zonalave_timeave(θ)
-    T2_timeave_zonalave = calc_zonalave_variance(T)
-    u2_timeave_zonalave = calc_zonalave_variance(u)
-    v2_timeave_zonalave = calc_zonalave_variance(v)
-    w2_timeave_zonalave = calc_zonalave_variance(w)
-    uT_cov_timeave_zonalave = calc_zonalave_covariance(u, T)
-    wT_cov_timeave_zonalave = calc_zonalave_covariance(w, T)
-
-    if moist
-        qt_timeave_zonalave = calc_zonalave_timeave(qt)
-        uqt_cov_timeave_zonalave = calc_zonalave_covariance(u, qt)
-        wqt_cov_timeave_zonalave = calc_zonalave_covariance(w, qt)
-    end
-
-    # plot!!
-    fig = generate_empty_figure(; resolution = (4098, 4098))
-    create_plot!(
-        fig;
-        p_loc = (1, 1),
-        X = lat,
-        Y = z,
-        Z = u_timeave_zonalave,
-        yscale = log10,
-        title = "Zonal Velocity",
-        xlabel = "Latitude",
-        ylabel = "Altitude[m]",
-    ) # TODO Make lims (-30:30)
-    create_plot!(
-        fig;
-        p_loc = (1, 2),
-        X = lat,
-        Y = z,
-        Z = T_timeave_zonalave,
-        yscale = log10,
-        title = "Temperature T",
-        xlabel = "Latitude",
-        ylabel = "Altitude[m]",
-    ) # TODO Make lims (190:10:310)
-    create_plot!(
-        fig;
-        p_loc = (1, 3),
-        X = lat,
-        Y = z,
-        Z = θ_timeave_zonalave,
-        yscale = log10,
-        title = "Potential Temperature θ",
-        xlabel = "Latitude",
-        ylabel = "Altitude[m]",
-    ) # TODO Make lims (190:10:310)
-    create_plot!(
-        fig;
-        p_loc = (2, 1),
-        X = lat,
-        Y = z,
-        Z = u2_timeave_zonalave,
-        yscale = log10,
-        title = "Zonal Vel Variance u′²",
-        xlabel = "Latitude",
-        ylabel = "Altitude[m]",
-    ) # TODO Make lims (0:40)
-    create_plot!(
-        fig;
-        p_loc = (2, 2),
-        X = lat,
-        Y = z,
-        Z = T2_timeave_zonalave,
-        yscale = log10,
-        title = "Temp Variance T′²",
-        xlabel = "Latitude",
-        ylabel = "Altitude[m]",
-    ) # TODO Make lims (0:40)
-    create_plot!(
-        fig;
-        p_loc = (3, 1),
-        X = lat,
-        Y = z,
-        Z = uT_cov_timeave_zonalave,
-        yscale = log10,
-        title = "Covariance u′T′",
-        xlabel = "Latitude",
-        ylabel = "Altitude[m]",
-    ) # TODO Make lims (0:40)
-    create_plot!(
-        fig;
-        p_loc = (3, 2),
-        X = lon,
-        Y = lat,
-        Z = T[:, :, 1, end],
-        title = "Temperature: Day $(days[end]); z $(round(z[1])) m",
-        xlabel = "Longitude",
-        ylabel = "Latitude",
-    ) # TODO Make lims (0:40)
-    create_plot!(
-        fig;
-        p_loc = (3, 3),
-        X = lat,
-        Y = z,
-        Z = wT_cov_timeave_zonalave,
-        yscale = log10,
-        title = "Covariance w′T′",
-        xlabel = "Latitude",
-        ylabel = "Altitude[m]",
-    ) # TODO Make lims (0:40)
-    if moist
-        create_plot!(
-            fig;
-            p_loc = (4, 1),
-            X = lat,
-            Y = z,
-            Z = qt_timeave_zonalave .* 1e3,
-            yscale = log10,
-            title = "qt × 10³ (Zonal-Time-Averaged)",
-            xlabel = "Latitude",
-            ylabel = "Altitude[m]",
-        ) # TODO Make lims (0:40)
-        create_plot!(
-            fig;
-            p_loc = (4, 2),
-            X = lon,
-            Y = lat,
-            Z = qt[:, :, 1, end] .* 1e3,
-            title = "qt × 10³: Day $(days[end]); z $(round(z[1])) m",
-            xlabel = "Longitude",
-            ylabel = "Latitude",
-        ) # TODO Make lims (0:40)
-        create_plot!(
-            fig;
-            p_loc = (4, 3),
-            X = lat,
-            Y = z,
-            Z = uqt_cov_timeave_zonalave,
-            yscale = log10,
-            title = "Covariance u′q′ (Zonal-Time-Averaged)",
-            xlabel = "Latitude",
-            ylabel = "Altitude[m]",
-        ) # TODO Make lims (0:40)
-        create_plot!(
-            fig;
-            p_loc = (2, 3),
-            X = lat,
-            Y = z,
-            Z = wqt_cov_timeave_zonalave,
-            yscale = log10,
-            title = "Covariance w′qt′",
-            xlabel = "Latitude",
-            ylabel = "Altitude[m]",
-        ) # TODO Make lims (0:40)
-    end
-
-    CairoMakie.save(fig_dir * "/diagnostics.png", fig)
-
-    # sfc_flux_energy, sfc_evaporation, lw_flux_net, sw_flux_net are set to
-    # nothing if there's no flux information
-    has_flux = any(
-        map(
-            !isnothing,
-            [sfc_flux_energy, sfc_evaporation, lw_flux_net, sw_flux_net],
-        ),
-    )
-
-    if has_flux
-        fig_flux = generate_empty_figure(; resolution = (2048, 2048))
-    end
-
-    if !isnothing(sfc_flux_energy)
-        create_plot!(
-            fig_flux;
-            p_loc = (1, 1),
-            X = lon,
-            Y = lat,
-            Z = sfc_flux_energy[:, :, end],
-            title = "Surface flux of energy: Day $(days[end]);",
-            xlabel = "Longitude",
-            ylabel = "Latitude",
-        )
-    end
-    if !isnothing(sfc_evaporation)
-        create_plot!(
-            fig_flux;
-            p_loc = (1, 2),
-            X = lon,
-            Y = lat,
-            Z = sfc_evaporation[:, :, end],
-            title = "Surface evaporation: Day $(days[end]);",
-            xlabel = "Longitude",
-            ylabel = "Latitude",
-        )
-    end
-    if !isnothing(lw_flux_net)
-        create_plot!(
-            fig_flux;
-            p_loc = (2, 1),
-            X = lon,
-            Y = lat,
-            # There's a minus because we define net = upward - downward
-            Z = -lw_flux_net[:, :, end, end],
-            title = "Net (down - up) longwave flux at TOA: Day $(days[end]);",
-            xlabel = "Longitude",
-            ylabel = "Latitude",
-        )
-        create_plot!(
-            fig_flux;
-            p_loc = (2, 2),
-            X = lon,
-            Y = lat,
-            Z = lw_flux_net[:, :, 1, end],
-            title = "Net (up - down) longwave flux at surface: Day $(days[end]);",
-            xlabel = "Longitude",
-            ylabel = "Latitude",
-        )
-    end
-    if !isnothing(sw_flux_net)
-        create_plot!(
-            fig_flux;
-            p_loc = (3, 1),
-            X = lon,
-            Y = lat,
-            # There's a minus because we define net = upward - downward
-            Z = -sw_flux_net[:, :, end, end],
-            title = "Net (down - up) shortwave flux at TOA: Day $(days[end]);",
-            xlabel = "Longitude",
-            ylabel = "Latitude",
-        )
-        create_plot!(
-            fig_flux;
-            p_loc = (3, 2),
-            X = lon,
-            Y = lat,
-            Z = sw_flux_net[:, :, 1, end],
-            title = "Net (up - down) shortwave flux at surface: Day $(days[end]);",
-            xlabel = "Longitude",
-            ylabel = "Latitude",
-        )
-    end
-
-    if has_flux
-        CairoMakie.save(fig_dir * "/diagnostics_flux.png", fig_flux)
-    end
-end
diff --git a/post_processing/plot/plot_pipeline.jl b/post_processing/plot/plot_pipeline.jl
deleted file mode 100644
index 8d3a5a7629f..00000000000
--- a/post_processing/plot/plot_pipeline.jl
+++ /dev/null
@@ -1,65 +0,0 @@
-using CairoMakie
-
-include(joinpath(@__DIR__, "plot_helpers.jl"))
-
-import ArgParse
-function parse_commandline()
-    s = ArgParse.ArgParseSettings()
-    ArgParse.@add_arg_table s begin
-        "--nc_dir"
-        help = "NetCDF directory"
-        arg_type = String
-        "--fig_dir"
-        help = "Figure directory"
-        arg_type = String
-        "--case_name"
-        help = "Experiment Name"
-        arg_type = String
-    end
-    parsed_args = ArgParse.parse_args(ARGS, s)
-    return (s, parsed_args)
-end
-
-function get_params()
-    (s, parsed_args) = parse_commandline()
-    nc_dir = parsed_args["nc_dir"]
-    fig_dir = parsed_args["fig_dir"]
-    case_name = parsed_args["case_name"]
-    if isnothing(fig_dir)
-        fig_dir = joinpath(nc_dir, "fig")
-    end
-    if isnothing(case_name)
-        case_name = "aquaplanet"
-    end
-    mkpath(fig_dir)
-
-    nc_files = filter(
-        x -> startswith(basename(x), "day") && endswith(x, ".nc"),
-        readdir(nc_dir, join = true),
-    )
-
-    return (; nc_files, fig_dir, case_name)
-end
-
-(; nc_files, fig_dir, case_name) = get_params()
-
-if case_name == "dry_baroclinic_wave"
-    generate_paperplots(case_name, args...) =
-        generate_paperplots_dry_baro_wave(args...)
-elseif case_name == "moist_baroclinic_wave"
-    generate_paperplots(case_name, args...) =
-        generate_paperplots_moist_baro_wave(args...)
-elseif case_name == "dry_held_suarez"
-    generate_paperplots(case_name, args...) =
-        generate_paperplots_held_suarez(args...; moist = false)
-elseif case_name == "moist_held_suarez" || case_name == "aquaplanet"
-    generate_paperplots(case_name, args...) =
-        generate_paperplots_held_suarez(args...; moist = true)
-elseif case_name == "elevation_spectrum"
-    generate_paperplots(case_name, args...) =
-        generate_elevation_spectra(args...)
-else
-    error("Unknown `case_name`: $(case_name)")
-end
-
-generate_paperplots(case_name, fig_dir, nc_files)
diff --git a/post_processing/post_processing_funcs.jl b/post_processing/post_processing_funcs.jl
deleted file mode 100644
index 1b1eae488b5..00000000000
--- a/post_processing/post_processing_funcs.jl
+++ /dev/null
@@ -1,319 +0,0 @@
-include(joinpath(@__DIR__, "remap", "remap_helpers.jl"))
-
-space_string(::Spaces.FaceExtrudedFiniteDifferenceSpace) = "(Face field)"
-space_string(::Spaces.CenterExtrudedFiniteDifferenceSpace) = "(Center field)"
-
-import ClimaCoreTempestRemap: def_space_coord
-import ClimaCoreSpectra: power_spectrum_1d, power_spectrum_2d
-import ClimaCore: Geometry, Fields, Spaces
-using LinearAlgebra: norm
-import ClimaAtmos:
-    DryModel,
-    EquilMoistModel,
-    NonEquilMoistModel,
-    PotentialTemperature,
-    TotalEnergy
-
-function process_name(s::AbstractString)
-    # "c_ρ", "c_ρe", "c_uₕ_1", "c_uₕ_2", "f_w_1"
-    s = replace(s, "components_data_" => "")
-    s = replace(s, "ₕ" => "_h")
-    s = replace(s, "ρ" => "rho")
-    return s
-end
-processed_varname(pc::Tuple) = process_name(join(pc, "_"))
-
-# TODO: Make this a RecipesBase.@recipe
-function profile_animation(sol, output_dir, fps)
-    # Column animations
-    Y0 = first(sol.u)
-    for prop_chain in Fields.property_chains(Y0)
-        var_name = processed_varname(prop_chain)
-        var_space = axes(Fields.single_field(Y0, prop_chain))
-        Ni, Nj, _, _, Nh = size(ClimaCore.Spaces.local_geometry_data(var_space))
-        n_columns = Nh * Nj * Ni # TODO: is this correct?
-        @info(
-            "Creating profile animation",
-            n_columns,
-            var_name,
-            n_timesteps = length(sol.u)
-        )
-        anim = Plots.@animate for Y in sol.u
-            var = Fields.single_field(Y, prop_chain)
-            temporary = ClimaCore.column(var, 1, 1, 1)
-            ϕ_col_ave = deepcopy(vec(temporary))
-            ϕ_col_std = deepcopy(vec(temporary))
-            ϕ_col_ave .= 0
-            ϕ_col_std .= 0
-            local_geom = Fields.local_geometry_field(axes(var))
-            z_f = ClimaCore.column(local_geom, 1, 1, 1)
-            z_f = z_f.coordinates.z
-            z = vec(z_f)
-            for h in 1:Nh, j in 1:Nj, i in 1:Ni
-                ϕ_col = ClimaCore.column(var, i, j, h)
-                ϕ_col_ave .+= vec(ϕ_col) ./ n_columns
-            end
-            for h in 1:Nh, j in 1:Nj, i in 1:Ni
-                ϕ_col = ClimaCore.column(var, i, j, h)
-                ϕ_col_std .+=
-                    sqrt.((vec(ϕ_col) .- ϕ_col_ave) .^ 2 ./ n_columns)
-            end
-
-            # TODO: use xribbon when supported: https://github.com/JuliaPlots/Plots.jl/issues/2702
-            Plots.plot(
-                ϕ_col_ave,
-                z ./ 1000;
-                label = "Mean & Std",
-                grid = false,
-                xerror = ϕ_col_std,
-                fillalpha = 0.5,
-            )
-            Plots.plot!(;
-                xlabel = "$var_name",
-                ylabel = "z [km]",
-                markershape = :circle,
-            )
-            Plots.title!("$(space_string(var_space))")
-        end
-        Plots.mp4(
-            anim,
-            joinpath(output_dir, "profile_$var_name.mp4"),
-            fps = fps,
-        )
-    end
-end
-
-function contour_animations(sol, output_dir, fps)
-    for prop_chain in Fields.property_chains(sol.u[end])
-        var_name = processed_varname(prop_chain)
-        @info "Creating contour animation" var_name n_timesteps = length(sol.u)
-        anim = Plots.@animate for Y in sol.u
-            var = Fields.single_field(Y, prop_chain)
-            level = 3
-            # TODO: do not use ClimaCore internals
-            if axes(var) isa Spaces.FaceExtrudedFiniteDifferenceSpace
-                level = ClimaCore.Utilities.PlusHalf(level)
-            end
-            clim = (minimum(var), maximum(var))
-            Plots.plot(var, level = level, clim = clim)
-        end
-        Plots.mp4(
-            anim,
-            joinpath(output_dir, "contour_$var_name.mp4"),
-            fps = fps,
-        )
-    end
-end
-
-function postprocessing_box(sol, output_dir)
-    for prop_chain in Fields.property_chains(sol.u[1])
-        var_name = processed_varname(prop_chain)
-        t_start = sol.t[1]
-        var_start = Fields.single_field(sol.u[1], prop_chain)
-        t_end = sol.t[end]
-        var_end = Fields.single_field(sol.u[end], prop_chain)
-        @info(
-            "L₂ norm",
-            var_name,
-            t_start,
-            norm(var_start),
-            t_end,
-            norm(var_end)
-        )
-    end
-
-    Y = sol.u[end]
-    FT = Spaces.undertype(axes(Y.c))
-    y_max = maximum(Fields.coordinate_field(axes(Y.f)).y)
-    ᶠw = Geometry.WVector.(Y.f.u₃).components.data.:1
-    p = Plots.plot(ᶠw, slice = (:, FT(y_max / 2), :), clim = (-0.1, 0.1))
-    Plots.png(p, joinpath(output_dir, "w.png"))
-
-    if :ρq_tot in propertynames(Y.c)
-        qt = Y.c.ρq_tot ./ Y.c.ρ
-        pq =
-            Plots.plot(qt, slice = (:, FT(y_max / 2), :), clim = (-0.02, 0.020))
-        Plots.png(pq, joinpath(output_dir, "qt.png"))
-    end
-end
-
-function postprocessing(sol, output_dir, fps)
-    for prop_chain in Fields.property_chains(sol.u[1])
-        var_name = processed_varname(prop_chain)
-        t_start = sol.t[1]
-        var_start = Fields.single_field(sol.u[1], prop_chain)
-        t_end = sol.t[end]
-        var_end = Fields.single_field(sol.u[end], prop_chain)
-        @info(
-            "L₂ norm",
-            var_name,
-            t_start,
-            norm(var_start),
-            t_end,
-            norm(var_end)
-        )
-    end
-
-    ᶠw_max = maximum(
-        map(u -> maximum(parent(ClimaCore.Geometry.WVector.(u.f.w))), sol.u),
-    )
-    ᶠw_min = minimum(
-        map(u -> minimum(parent(ClimaCore.Geometry.WVector.(u.f.w))), sol.u),
-    )
-    @info "maximum vertical velocity" ᶠw_max
-    @info "maximum vertical velocity" ᶠw_min
-
-    # contour_animations(sol, output_dir, fps) # For generic contours:
-
-    anim = Plots.@animate for Y in sol.u
-        ᶜv = Geometry.UVVector.(Y.c.uₕ).components.data.:2
-        Plots.plot(ᶜv, level = 3, clim = (-6, 6))
-    end
-    Plots.mp4(anim, joinpath(output_dir, "v.mp4"), fps = fps)
-
-    anim = Plots.@animate for Y in sol.u
-        ᶠw = Geometry.WVector.(Y.f.u₃).components.data.:1
-        Plots.plot(
-            ᶠw,
-            level = ClimaCore.Utilities.PlusHalf(3),
-            clim = (-0.02, 0.02),
-        )
-    end
-    Plots.mp4(anim, joinpath(output_dir, "w.mp4"), fps = fps)
-
-    prop_chains = Fields.property_chains(sol.u[1])
-    FT = Spaces.undertype(axes(Y.c))
-    if any(pc -> pc == (:c, :ρq_tot), prop_chains)
-        anim = Plots.@animate for Y in sol.u
-            ᶜq_tot = Y.c.ρq_tot ./ Y.c.ρ
-            Plots.plot(ᶜq_tot .* FT(1e3), level = 3, clim = (0, 1))
-        end
-        Plots.mp4(anim, joinpath(output_dir, "contour_q_tot.mp4"), fps = fps)
-    else
-        @info "Moisture not found" prop_chains
-    end
-
-    profile_animation(sol, output_dir, fps)
-end
-
-function safe_index(ius, t)
-    iu = if isempty(ius)
-        @warn "Could not find desired time for plotting, falling back on last day."
-        length(t)
-    else
-        first(ius)
-    end
-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)
-    n = length(sol.u)
-    #! format: off
-    get_row(var) = [
-        "Y.$(join(var, '.'))";;
-        "$(norm(get_var(1, var), 2)) → $(norm(get_var(n, var), 2))";;
-        "$(mean(get_var(1, var))) → $(mean(get_var(n, var)))";;
-        "$(maximum(abs, get_var(1, var))) → $(maximum(abs, get_var(n, var)))";;
-        "$(minimum(abs, get_var(1, var))) → $(minimum(abs, get_var(n, var)))";;
-    ]
-    #! format: on
-    pretty_table(
-        vcat(map(get_row, Fields.property_chains(sol.u[1]))...);
-        title = "Change in Y from t = $(sol.t[1]) to t = $(sol.t[n]):",
-        header = ["var", "‖var‖₂", "mean(var)", "max(∣var∣)", "min(∣var∣)"],
-        alignment = :c,
-    )
-
-    anim = @animate for (Y, t) in zip(sol.u, sol.t)
-        CA.set_precomputed_quantities!(Y, p, t) # sets ᶜts
-        Plots.plot(
-            vec(TD.air_temperature.(thermo_params, p.precomputed.ᶜts)),
-            vec(Fields.coordinate_field(Y.c).z ./ 1000);
-            xlabel = "T [K]",
-            ylabel = "z [km]",
-            xlims = (190, 310),
-            legend = false,
-        )
-    end
-    Plots.mp4(anim, joinpath(output_dir, "T.mp4"), fps = 10)
-
-    anim = @animate for Y in sol.u
-        w = Geometry.WVector.(Y.f.u₃).components.data.:1
-        Plots.plot(
-            vec(w),
-            vec(Fields.coordinate_field(Y.f).z ./ 1000);
-            xlabel = "w [m/s]",
-            ylabel = "z [km]",
-            xlims = (-0.02, 0.02),
-            legend = false,
-        )
-    end
-    Plots.mp4(anim, joinpath(output_dir, "w.mp4"), fps = 10)
-end
-
-function postprocessing_plane(sol, output_dir, p)
-    thermo_params = CAP.thermodynamics_params(p.params)
-    get_var(i, var) = Fields.single_field(sol.u[i], var)
-    n = length(sol.u)
-    #! format: off
-    get_row(var) = [
-        "Y.$(join(var, '.'))";;
-        "$(norm(get_var(1, var), 2)) → $(norm(get_var(n, var), 2))";;
-        "$(mean(get_var(1, var))) → $(mean(get_var(n, var)))";;
-        "$(maximum(abs, get_var(1, var))) → $(maximum(abs, get_var(n, var)))";;
-        "$(minimum(abs, get_var(1, var))) → $(minimum(abs, get_var(n, var)))";;
-    ]
-    #! format: on
-    pretty_table(
-        vcat(map(get_row, Fields.property_chains(sol.u[1]))...);
-        title = "Change in Y from t = $(sol.t[1]) to t = $(sol.t[n]):",
-        header = ["var", "‖var‖₂", "mean(var)", "max(∣var∣)", "min(∣var∣)"],
-        alignment = :c,
-    )
-
-    Y = sol.u[end]
-
-    ## Plots for last timestep
-    function gen_plot_plane(
-        variable::Fields.Field,
-        filename::String,
-        title::String,
-        xlabel::String,
-        ylabel::String;
-        output_dir = output_dir,
-    )
-        # Set up Figure and Axes
-        f = CairoMakie.Figure(; font = "CMU Serif")
-        gaa = f[1, 1] = GridLayout()
-        Axis(gaa[1, 1], aspect = 2, title = title)
-        paa = fieldcontourf!(variable)
-        Colorbar(gaa[1, 2], paa, label = xlabel)
-        fig_png = joinpath(output_dir, filename)
-        CairoMakie.save(fig_png, f)
-    end
-
-    gen_plot_plane(
-        Geometry.UVector.(p.precomputed.ᶜu),
-        "horz_velocity.png",
-        "Horizontal Velocity",
-        "u[m/s]",
-        "z[m]",
-    )
-
-    gen_plot_plane(
-        Geometry.WVector.(p.precomputed.ᶜu),
-        "vert_velocity.png",
-        "Vertical Velocity",
-        "w[m/s]",
-        "z[m]",
-    )
-
-    gen_plot_plane(
-        TD.virtual_pottemp.(thermo_params, p.precomputed.ᶜts),
-        "virtual_pottemp.png",
-        "Virtual Pottemp",
-        "Theta[K]",
-        "z[m]",
-    )
-end
diff --git a/post_processing/remap/remap_helpers.jl b/post_processing/remap/remap_helpers.jl
index 2a17dd6b73f..55d54a0cbd9 100644
--- a/post_processing/remap/remap_helpers.jl
+++ b/post_processing/remap/remap_helpers.jl
@@ -79,7 +79,6 @@ function remap2latlon(filein, data_dir, remap_tmpdir, weightfile, nlat, nlon)
         ClimaComms.SingletonCommsContext(ClimaComms.CPUSingleThreaded()),
     )
     Y = InputOutput.read_field(reader, "Y")
-    diag = InputOutput.read_field(reader, "diagnostics")
     t_now = InputOutput.HDF5.read_attribute(reader.file, "time")
 
     remap_tmpsubdir = if @isdefined pmap
@@ -121,76 +120,6 @@ function remap2latlon(filein, data_dir, remap_tmpdir, weightfile, nlat, nlon)
     nc_u = defVar(nc, "u", FT, cspace, ("time",))
     nc_v = defVar(nc, "v", FT, cspace, ("time",))
     nc_w = defVar(nc, "w", FT, cspace, ("time",))
-    # define variables for the diagnostic states
-    nc_pres = defVar(nc, "pressure", FT, cspace, ("time",))
-    nc_T = defVar(nc, "temperature", FT, cspace, ("time",))
-    nc_θ = defVar(nc, "potential_temperature", FT, cspace, ("time",))
-    nc_K = defVar(nc, "kinetic_energy", FT, cspace, ("time",))
-    nc_vort = defVar(nc, "vorticity", FT, cspace, ("time",))
-    nc_T_sfc = defVar(nc, "sfc_temperature", FT, hspace, ("time",))
-    nc_z_sfc = defVar(nc, "sfc_elevation", FT, hspace, ("time",))
-    nc_qt_sfc = defVar(nc, "sfc_qt", FT, hspace, ("time",))
-    # define moist variables
-    if :ρq_tot in propertynames(Y.c)
-        nc_qt = defVar(nc, "qt", FT, cspace, ("time",))
-        nc_RH = defVar(nc, "RH", FT, cspace, ("time",))
-        nc_cloudliq = defVar(nc, "cloud_liquid", FT, cspace, ("time",))
-        nc_cloudice = defVar(nc, "cloud_ice", FT, cspace, ("time",))
-        nc_watervapor = defVar(nc, "water_vapor", FT, cspace, ("time",))
-        nc_cloudfrac_gm = defVar(nc, "cloud_fraction_gm", FT, cspace, ("time",))
-    end
-    # define precip variables
-    if :precipitation_removal in propertynames(diag)
-        nc_precipitation_removal =
-            defVar(nc, "precipitation_removal", FT, cspace, ("time",))
-        nc_column_integrated_rain =
-            defVar(nc, "column_integrated_rain", FT, hspace, ("time",))
-        nc_column_integrated_snow =
-            defVar(nc, "column_integrated_snow", FT, hspace, ("time",))
-    end
-    # define surface flux variables
-    if :sfc_flux_energy in propertynames(diag)
-        nc_sfc_flux_energy =
-            defVar(nc, "sfc_flux_energy", FT, hspace, ("time",))
-        nc_sfc_flux_u = defVar(nc, "sfc_flux_u", FT, hspace, ("time",))
-        nc_sfc_flux_v = defVar(nc, "sfc_flux_v", FT, hspace, ("time",))
-        if :sfc_evaporation in propertynames(diag)
-            nc_sfc_evaporation =
-                defVar(nc, "sfc_evaporation", FT, hspace, ("time",))
-        end
-    end
-    # define radiative flux variables
-    if :lw_flux_down in propertynames(diag)
-        nc_lw_flux_down = defVar(nc, "lw_flux_down", FT, fspace, ("time",))
-        nc_lw_flux_up = defVar(nc, "lw_flux_up", FT, fspace, ("time",))
-        nc_sw_flux_down = defVar(nc, "sw_flux_down", FT, fspace, ("time",))
-        nc_sw_flux_up = defVar(nc, "sw_flux_up", FT, fspace, ("time",))
-    end
-    if :clear_lw_flux_down in propertynames(diag)
-        nc_clear_lw_flux_down =
-            defVar(nc, "clear_lw_flux_down", FT, fspace, ("time",))
-        nc_clear_lw_flux_up =
-            defVar(nc, "clear_lw_flux_up", FT, fspace, ("time",))
-        nc_clear_sw_flux_down =
-            defVar(nc, "clear_sw_flux_down", FT, fspace, ("time",))
-        nc_clear_sw_flux_up =
-            defVar(nc, "clear_sw_flux_up", FT, fspace, ("time",))
-    end
-    if :draft_q_tot in propertynames(diag)
-        nc_draft_qt = defVar(nc, "draft_qt", FT, cspace, ("time",))
-        nc_draft_RH = defVar(nc, "draft_RH", FT, cspace, ("time",))
-        nc_draft_cloudliq =
-            defVar(nc, "draft_cloud_liquid", FT, cspace, ("time",))
-        nc_draft_cloudice = defVar(nc, "draft_cloud_ice", FT, cspace, ("time",))
-        nc_draft_watervapor =
-            defVar(nc, "draft_water_vapor", FT, cspace, ("time",))
-        nc_draft_cloudfrac =
-            defVar(nc, "draft_cloud_fraction", FT, cspace, ("time",))
-        nc_draft_area = defVar(nc, "draft_area_fraction", FT, cspace, ("time",))
-        nc_cloudfrac = defVar(nc, "cloud_fraction", FT, cspace, ("time",))
-        nc_tke = defVar(nc, "tke", FT, cspace, ("time",))
-        nc_mixing_length = defVar(nc, "mixing_length", FT, cspace, ("time",))
-    end
 
     # time
     nc_time[1] = t_now
@@ -199,11 +128,8 @@ function remap2latlon(filein, data_dir, remap_tmpdir, weightfile, nlat, nlon)
     # density
     nc_rho[:, 1] = Y.c.ρ
     # thermodynamics
-    if :ρe_tot in propertynames(Y.c)
-        nc_thermo[:, 1] = Y.c.ρe_tot ./ Y.c.ρ
-    elseif :ρθ in propertynames(Y.c)
-        nc_thermo[:, 1] = Y.c.ρθ ./ Y.c.ρ
-    end
+    nc_thermo[:, 1] = Y.c.ρe_tot ./ Y.c.ρ
+
     # physical horizontal velocity
     uh_phy = Geometry.transform.(tuple(Geometry.UVAxis()), Y.c.uₕ)
     nc_u[:, 1] = uh_phy.components.data.:1
@@ -212,69 +138,9 @@ function remap2latlon(filein, data_dir, remap_tmpdir, weightfile, nlat, nlon)
     ᶠw = Geometry.WVector.(Y.f.u₃)
     ᶜw = ᶜinterp.(ᶠw)
     nc_w[:, 1] = ᶜw
-    # diagnostic variables
-    nc_pres[:, 1] = diag.pressure
-    nc_T[:, 1] = diag.temperature
-    nc_θ[:, 1] = diag.potential_temperature
-    nc_K[:, 1] = diag.kinetic_energy
-    nc_vort[:, 1] = diag.vorticity
-    nc_T_sfc[:, 1] = diag.sfc_temperature
-    nc_z_sfc[:, 1] = Fields.level(Fields.coordinate_field(Y.f.u₃).z, half)
-    nc_qt_sfc[:, 1] = diag.sfc_qt
-
-    if :ρq_tot in propertynames(Y.c)
-        nc_qt[:, 1] = Y.c.ρq_tot ./ Y.c.ρ
-        nc_RH[:, 1] = diag.relative_humidity
-        nc_cloudliq[:, 1] = diag.q_liq
-        nc_cloudice[:, 1] = diag.q_ice
-        nc_watervapor[:, 1] = diag.q_vap
-        nc_cloudfrac_gm[:, 1] = diag.cloud_fraction_gm
-    end
-
-    if :precipitation_removal in propertynames(diag)
-        nc_precipitation_removal[:, 1] = diag.precipitation_removal
-        nc_column_integrated_rain[:, 1] = diag.column_integrated_rain
-        nc_column_integrated_snow[:, 1] = diag.column_integrated_snow
-    end
-
-    if :sfc_flux_energy in propertynames(diag)
-        nc_sfc_flux_energy[:, 1] = diag.sfc_flux_energy
-        nc_sfc_flux_u[:, 1] = diag.sfc_flux_u
-        nc_sfc_flux_v[:, 1] = diag.sfc_flux_v
-        if :sfc_evaporation in propertynames(diag)
-            nc_sfc_evaporation[:, 1] = diag.sfc_evaporation
-        end
-    end
-
-    if :lw_flux_down in propertynames(diag)
-        nc_lw_flux_down[:, 1] = diag.lw_flux_down
-        nc_lw_flux_up[:, 1] = diag.lw_flux_up
-        nc_sw_flux_down[:, 1] = diag.sw_flux_down
-        nc_sw_flux_up[:, 1] = diag.sw_flux_up
-    end
 
-    if :clear_lw_flux_down in propertynames(diag)
-        nc_clear_lw_flux_down[:, 1] = diag.clear_lw_flux_down
-        nc_clear_lw_flux_up[:, 1] = diag.clear_lw_flux_up
-        nc_clear_sw_flux_down[:, 1] = diag.clear_sw_flux_down
-        nc_clear_sw_flux_up[:, 1] = diag.clear_sw_flux_up
-    end
-
-    if :draft_q_tot in propertynames(diag)
-        nc_draft_qt[:, 1] = diag.draft_q_tot
-        nc_draft_RH[:, 1] = diag.draft_relative_humidity
-        nc_draft_cloudliq[:, 1] = diag.draft_q_liq
-        nc_draft_cloudice[:, 1] = diag.draft_q_ice
-        nc_draft_watervapor[:, 1] = diag.draft_q_vap
-        nc_draft_cloudfrac[:, 1] = diag.draft_cloud_fraction
-        nc_draft_area[:, 1] = diag.draft_area
-        nc_cloudfrac[:, 1] = diag.cloud_fraction
-        nc_tke[:, 1] = diag.env_tke
-        nc_mixing_length[:, 1] = diag.env_mixing_length
-    end
     close(nc)
 
-
     datafile_latlon =
         joinpath(out_dir, first(splitext(basename(filein))) * ".nc")
     dry_variables = [
@@ -283,85 +149,10 @@ function remap2latlon(filein, data_dir, remap_tmpdir, weightfile, nlat, nlon)
         "u",
         "v",
         "w",
-        "pressure",
-        "temperature",
-        "potential_temperature",
-        "kinetic_energy",
-        "vorticity",
-        "sfc_temperature",
-        "sfc_elevation",
-        "sfc_qt",
     ]
-    if :ρq_tot in propertynames(Y.c)
-        moist_variables = [
-            "qt",
-            "RH",
-            "cloud_ice",
-            "cloud_liquid",
-            "water_vapor",
-            "cloud_fraction_gm",
-        ]
-    else
-        moist_variables = String[]
-    end
-    if :precipitation_removal in propertynames(diag)
-        precip_variables = [
-            "precipitation_removal",
-            "column_integrated_rain",
-            "column_integrated_snow",
-        ]
-    else
-        precip_variables = String[]
-    end
-    if :sfc_flux_energy in propertynames(diag)
-        sfc_flux_variables = ["sfc_flux_energy", "sfc_flux_u", "sfc_flux_v"]
-        if :sfc_evaporation in propertynames(diag)
-            sfc_flux_variables = [sfc_flux_variables..., "sfc_evaporation"]
-        end
-    else
-        sfc_flux_variables = String[]
-    end
-    if :lw_flux_down in propertynames(diag)
-        rad_flux_variables =
-            ["lw_flux_down", "lw_flux_up", "sw_flux_down", "sw_flux_up"]
-    else
-        rad_flux_variables = String[]
-    end
-    if :clear_lw_flux_down in propertynames(diag)
-        rad_flux_clear_variables = [
-            "clear_lw_flux_down",
-            "clear_lw_flux_up",
-            "clear_sw_flux_down",
-            "clear_sw_flux_up",
-        ]
-    else
-        rad_flux_clear_variables = String[]
-    end
-    if :draft_q_tot in propertynames(diag)
-        edmf_variables = [
-            "draft_qt",
-            "draft_RH",
-            "draft_cloud_ice",
-            "draft_cloud_liquid",
-            "draft_water_vapor",
-            "draft_cloud_fraction",
-            "draft_area_fraction",
-            "cloud_fraction",
-            "tke",
-            "mixing_length",
-        ]
-    else
-        edmf_variables = String[]
-    end
 
     netcdf_variables = vcat(
         dry_variables,
-        moist_variables,
-        precip_variables,
-        sfc_flux_variables,
-        rad_flux_variables,
-        rad_flux_clear_variables,
-        edmf_variables,
     )
     apply_remap(datafile_latlon, datafile_cc, weightfile, netcdf_variables)
     rm(datafile_cc)
diff --git a/post_processing/remap/remap_pipeline.jl b/post_processing/remap/remap_pipeline.jl
index c2f69248868..36803b1e942 100644
--- a/post_processing/remap/remap_pipeline.jl
+++ b/post_processing/remap/remap_pipeline.jl
@@ -9,7 +9,7 @@ Sbatch script example:
 #SBATCH --time=01:00:00
 #SBATCH --output=remap.log
 module purge
-module load julia/1.8.5
+module load julia/1.9.3
 export JULIA_CUDA_USE_BINARYBUILDER=false
 export JULIA_NUM_THREADS=${SLURM_CPUS_PER_TASK:=1}
 export HDF5_DIR=/central/groups/esm/{username}/ClimaAtmos_remap/
diff --git a/src/callbacks/callbacks.jl b/src/callbacks/callbacks.jl
index f0af2081d0d..af409f875a4 100644
--- a/src/callbacks/callbacks.jl
+++ b/src/callbacks/callbacks.jl
@@ -172,267 +172,19 @@ NVTX.@annotate function rrtmgp_model_callback!(integrator)
     return nothing
 end
 
-function common_diagnostics(p, ᶜu, ᶜts)
-    (; env_thermo_quad, params) = p
-    thermo_params = CAP.thermodynamics_params(params)
-    ᶜρ = TD.air_density.(thermo_params, ᶜts)
-    diagnostics = (;
-        u_velocity = Geometry.UVector.(ᶜu),
-        v_velocity = Geometry.VVector.(ᶜu),
-        w_velocity = Geometry.WVector.(ᶜu),
-        temperature = TD.air_temperature.(thermo_params, ᶜts),
-        potential_temperature = TD.dry_pottemp.(thermo_params, ᶜts),
-        specific_enthalpy = TD.specific_enthalpy.(thermo_params, ᶜts),
-        buoyancy = CAP.grav(p.params) .* (p.core.ᶜρ_ref .- ᶜρ) ./ ᶜρ,
-        density = TD.air_density.(thermo_params, ᶜts),
-    )
-    if !(p.atmos.moisture_model isa DryModel)
-        diagnostics = (;
-            diagnostics...,
-            q_vap = TD.vapor_specific_humidity.(thermo_params, ᶜts),
-            q_liq = TD.liquid_specific_humidity.(thermo_params, ᶜts),
-            q_ice = TD.ice_specific_humidity.(thermo_params, ᶜts),
-            q_tot = TD.total_specific_humidity.(thermo_params, ᶜts),
-            relative_humidity = TD.relative_humidity.(thermo_params, ᶜts),
-            cloud_fraction_gm = get_cloud_fraction.(
-                thermo_params,
-                env_thermo_quad,
-                p.precomputed.ᶜp,
-                ᶜts,
-            ),
-        )
-    end
-    return diagnostics
-end
-
-# Adds a prefix to the front of each name in the named tuple. This function
-# is not type stable, but that probably doesn't matter for diagnostics.
-add_prefix(diagnostics::NamedTuple{names}, prefix) where {names} =
-    NamedTuple{Symbol.(prefix, names)}(values(diagnostics))
-
-cloud_fraction(thermo_params, ts, area::FT) where {FT} =
-    TD.has_condensate(thermo_params, ts) && area > 1e-3 ? FT(1) : FT(0)
-
-NVTX.@annotate function compute_diagnostics(integrator)
-    (; t, u, p) = integrator
-    Y = u
-    (; params, env_thermo_quad) = p
-    FT = eltype(params)
-    thermo_params = CAP.thermodynamics_params(params)
-
-    (; ᶜu, ᶜK, ᶜts, ᶜp, sfc_conditions) = p.precomputed
-    dycore_diagnostic = (;
-        common_diagnostics(p, ᶜu, ᶜts)...,
-        pressure = ᶜp,
-        kinetic_energy = ᶜK,
-        sfc_temperature = TD.air_temperature.(thermo_params, sfc_conditions.ts),
-        sfc_qt = TD.total_specific_humidity.(thermo_params, sfc_conditions.ts),
-    )
-
-    if eltype(Fields.coordinate_field(axes(Y.c))) <: Geometry.Abstract3DPoint
-        ᶜvort = @. Geometry.WVector(curlₕ(Y.c.uₕ))
-        if p.do_dss
-            Spaces.weighted_dss!(ᶜvort)
-        end
-        dycore_diagnostic = (; dycore_diagnostic..., vorticity = ᶜvort)
-    end
-
-    if p.atmos.precip_model isa Microphysics0Moment
-        (; ᶜS_ρq_tot, col_integrated_rain, col_integrated_snow) =
-            p.precipitation
-        Fields.bycolumn(axes(Y.c)) do colidx
-            precipitation_tendency!(
-                nothing,
-                Y,
-                p,
-                t,
-                colidx,
-                p.atmos.precip_model,
-            )
-        end # TODO: Set the diagnostics without computing the tendency.
-        precip_diagnostic = (;
-            precipitation_removal = ᶜS_ρq_tot,
-            column_integrated_rain = col_integrated_rain,
-            column_integrated_snow = col_integrated_snow,
-        )
-    elseif p.atmos.precip_model isa Microphysics1Moment
-        # TODO: Get column integrals for the land model.
-        precip_diagnostic =
-            (; q_rai = Y.c.ρq_rai ./ Y.c.ρ, q_sno = Y.c.ρq_sno ./ Y.c.ρ)
-    else
-        precip_diagnostic = NamedTuple()
-    end
-
-    if p.atmos.turbconv_model isa PrognosticEDMFX
-        (; ᶜtke⁰, ᶜu⁰, ᶜts⁰, ᶜmixing_length) = p.precomputed
-        (; ᶜu⁺, ᶜts⁺, ᶜa⁺, ᶜa⁰) = output_prognostic_sgs_quantities(Y, p, t)
-        env_diagnostics = (;
-            common_diagnostics(p, ᶜu⁰, ᶜts⁰)...,
-            area = ᶜa⁰,
-            cloud_fraction = get_cloud_fraction.(
-                thermo_params,
-                env_thermo_quad,
-                ᶜp,
-                ᶜts⁰,
-            ),
-            tke = ᶜtke⁰,
-            mixing_length = ᶜmixing_length,
-        )
-        draft_diagnostics = (;
-            common_diagnostics(p, ᶜu⁺, ᶜts⁺)...,
-            area = ᶜa⁺,
-            cloud_fraction = cloud_fraction.(thermo_params, ᶜts⁺, ᶜa⁺),
-        )
-        turbulence_convection_diagnostic = (;
-            add_prefix(env_diagnostics, :env_)...,
-            add_prefix(draft_diagnostics, :draft_)...,
-            cloud_fraction = ᶜa⁰ .*
-                             get_cloud_fraction.(
-                thermo_params,
-                env_thermo_quad,
-                ᶜp,
-                ᶜts⁰,
-            ) .+ ᶜa⁺ .* cloud_fraction.(thermo_params, ᶜts⁺, ᶜa⁺),
-        )
-    elseif p.atmos.turbconv_model isa DiagnosticEDMFX
-        (; ᶜtke⁰, ᶜmixing_length) = p.precomputed
-        (; ᶜu⁺, ᶜts⁺, ᶜa⁺) = output_diagnostic_sgs_quantities(Y, p, t)
-        env_diagnostics = (;
-            cloud_fraction = get_cloud_fraction.(
-                thermo_params,
-                env_thermo_quad,
-                ᶜp,
-                ᶜts,
-            ),
-            tke = ᶜtke⁰,
-            mixing_length = ᶜmixing_length,
-        )
-        draft_diagnostics = (;
-            common_diagnostics(p, ᶜu⁺, ᶜts⁺)...,
-            area = ᶜa⁺,
-            cloud_fraction = cloud_fraction.(thermo_params, ᶜts⁺, ᶜa⁺),
-        )
-        turbulence_convection_diagnostic = (;
-            add_prefix(env_diagnostics, :env_)...,
-            add_prefix(draft_diagnostics, :draft_)...,
-            cloud_fraction = get_cloud_fraction.(
-                thermo_params,
-                env_thermo_quad,
-                ᶜp,
-                ᶜts,
-            ) .+ ᶜa⁺ .* cloud_fraction.(thermo_params, ᶜts⁺, ᶜa⁺),
-        )
-    else
-        turbulence_convection_diagnostic = NamedTuple()
-    end
-
-    if p.atmos.energy_form isa TotalEnergy
-        sfc_local_geometry =
-            Fields.level(Fields.local_geometry_field(Y.f), Fields.half)
-        surface_ct3_unit =
-            CT3.(unit_basis_vector_data.(CT3, sfc_local_geometry))
-        (; ρ_flux_uₕ, ρ_flux_h_tot) = p.precomputed.sfc_conditions
-        sfc_flux_momentum =
-            Geometry.UVVector.(
-                adjoint.(ρ_flux_uₕ ./ Spaces.level(ᶠinterp.(Y.c.ρ), half)) .*
-                surface_ct3_unit
-            )
-        vert_diff_diagnostic = (;
-            sfc_flux_u = sfc_flux_momentum.components.data.:1,
-            sfc_flux_v = sfc_flux_momentum.components.data.:2,
-            sfc_flux_energy = dot.(ρ_flux_h_tot, surface_ct3_unit),
-        )
-        if :ρq_tot in propertynames(Y.c)
-            (; ρ_flux_q_tot) = p.precomputed.sfc_conditions
-            vert_diff_diagnostic = (;
-                vert_diff_diagnostic...,
-                sfc_evaporation = dot.(ρ_flux_q_tot, surface_ct3_unit),
-            )
-        end
-    else
-        vert_diff_diagnostic = NamedTuple()
-    end
-
-    if p.atmos.radiation_mode isa RRTMGPI.AbstractRRTMGPMode
-        (; face_lw_flux_dn, face_lw_flux_up, face_sw_flux_dn, face_sw_flux_up) =
-            p.radiation.radiation_model
-        rad_diagnostic = (;
-            lw_flux_down = RRTMGPI.array2field(FT.(face_lw_flux_dn), axes(Y.f)),
-            lw_flux_up = RRTMGPI.array2field(FT.(face_lw_flux_up), axes(Y.f)),
-            sw_flux_down = RRTMGPI.array2field(FT.(face_sw_flux_dn), axes(Y.f)),
-            sw_flux_up = RRTMGPI.array2field(FT.(face_sw_flux_up), axes(Y.f)),
-        )
-        if p.atmos.radiation_mode isa
-           RRTMGPI.AllSkyRadiationWithClearSkyDiagnostics
-            (;
-                face_clear_lw_flux_dn,
-                face_clear_lw_flux_up,
-                face_clear_sw_flux_dn,
-                face_clear_sw_flux_up,
-            ) = p.radiation.radiation_model
-            rad_clear_diagnostic = (;
-                clear_lw_flux_down = RRTMGPI.array2field(
-                    FT.(face_clear_lw_flux_dn),
-                    axes(Y.f),
-                ),
-                clear_lw_flux_up = RRTMGPI.array2field(
-                    FT.(face_clear_lw_flux_up),
-                    axes(Y.f),
-                ),
-                clear_sw_flux_down = RRTMGPI.array2field(
-                    FT.(face_clear_sw_flux_dn),
-                    axes(Y.f),
-                ),
-                clear_sw_flux_up = RRTMGPI.array2field(
-                    FT.(face_clear_sw_flux_up),
-                    axes(Y.f),
-                ),
-            )
-        else
-            rad_clear_diagnostic = NamedTuple()
-        end
-    elseif p.atmos.radiation_mode isa RadiationDYCOMS_RF01
-        # TODO: add radiation diagnostics
-        rad_diagnostic = NamedTuple()
-        rad_clear_diagnostic = NamedTuple()
-    elseif p.atmos.radiation_mode isa RadiationTRMM_LBA
-        # TODO: add radiation diagnostics
-        rad_diagnostic = NamedTuple()
-        rad_clear_diagnostic = NamedTuple()
-    else
-        rad_diagnostic = NamedTuple()
-        rad_clear_diagnostic = NamedTuple()
-    end
-
-    diagnostic = merge(
-        dycore_diagnostic,
-        precip_diagnostic,
-        vert_diff_diagnostic,
-        rad_diagnostic,
-        rad_clear_diagnostic,
-        turbulence_convection_diagnostic,
-    )
-    return diagnostic
-end
-
 function save_to_disk_func(integrator)
     (; t, u, p) = integrator
     Y = u
-    diagnostic = compute_diagnostics(integrator)
 
     (; output_dir) = p.simulation
     day = floor(Int, t / (60 * 60 * 24))
     sec = floor(Int, t % (60 * 60 * 24))
-    @info "Saving diagnostics to HDF5 file on day $day second $sec"
+    @info "Saving state to HDF5 file on day $day second $sec"
     output_file = joinpath(output_dir, "day$day.$sec.hdf5")
     comms_ctx = ClimaComms.context(integrator.u.c)
     hdfwriter = InputOutput.HDF5Writer(output_file, comms_ctx)
     InputOutput.HDF5.write_attribute(hdfwriter.file, "time", t) # TODO: a better way to write metadata
     InputOutput.write!(hdfwriter, Y, "Y")
-    FT = Spaces.undertype(axes(Y.c))
-    values = map(Fields.wrap, diagnostic)
-    fv = Fields.FieldVector{FT}(values)
-    InputOutput.write!(hdfwriter, fv, "diagnostics")
     Base.close(hdfwriter)
     return nothing
 end