revamp ConservationChecker

add to driver

add to driver

conserves

clean up

Author: LenkaNovak

config/model_configs/interactive_debug.yml

@@ -4,9 +4,9 @@ rad: "gray"
 energy_check: true
 mode_name: "slabplanet"
 t_end: "10days"
-dt_save_to_sol: "3600secs"
-dt_cpl: 200
-dt: "200secs"
+dt_save_to_sol: "2days"
+dt_cpl: 400
+dt: "400secs"
 mono_surface: true
 turb_flux_partition: "CombinedStateFluxes"
 h_elem: 4

experiments/AMIP/modular/components/atmosphere/climaatmos_init.jl

@@ -71,9 +71,14 @@ end
 # extensions required by the Interfacer
 get_field(sim::ClimaAtmosSimulation, ::Val{:radiative_energy_flux}) =
     Fields.level(sim.integrator.p.ᶠradiation_flux, half)
-get_field(sim::ClimaAtmosSimulation, ::Val{:liquid_precipitation}) =
-    sim.integrator.p.col_integrated_rain .+ sim.integrator.p.col_integrated_snow # all fallen snow melts for now
-get_field(sim::ClimaAtmosSimulation, ::Val{:snow_precipitation}) = sim.integrator.p.col_integrated_snow .* FT(0)
+function get_field(sim::ClimaAtmosSimulation, ::Val{:liquid_precipitation})
+    ρ_liq = CAP.ρ_cloud_liq(sim.integrator.p.params)
+    sim.integrator.p.col_integrated_rain .* ρ_liq #.+ sim.integrator.p.col_integrated_snow .* ρ_liq # kg/m^2/s ; all fallen snow melts for now
+end
+function get_field(sim::ClimaAtmosSimulation, ::Val{:snow_precipitation}) # kg/m^2/s
+    ρ_liq = CAP.ρ_cloud_liq(sim.integrator.p.params)
+    sim.integrator.p.col_integrated_snow .* ρ_liq #.* FT(0)
+end
 
 get_field(sim::ClimaAtmosSimulation, ::Val{:turbulent_energy_flux}) =
     Geometry.WVector.(sim.integrator.p.sfc_conditions.ρ_flux_h_tot)
@@ -252,3 +257,40 @@ Extension of Checkpointer.get_model_state_vector to get the model state.
 function get_model_state_vector(sim::ClimaAtmosSimulation)
     return sim.integrator.u
 end
+
+function get_field(atmos_sim::ClimaAtmosSimulation, ::Val{:F_radiative_TOA})
+    radiation = atmos_sim.integrator.p.radiation_model
+    FT = eltype(atmos_sim.integrator.u)
+    # save radiation source
+    if radiation != nothing
+        face_space = axes(atmos_sim.integrator.u.f)
+        z = parent(Fields.coordinate_field(face_space).z)
+        Δz_top = round(FT(0.5) * (z[end, 1, 1, 1, 1] - z[end - 1, 1, 1, 1, 1]))
+        n_faces = length(z[:, 1, 1, 1, 1])
+
+        LWd_TOA = Fields.level(
+            RRTMGPI.array2field(FT.(atmos_sim.integrator.p.radiation_model.face_lw_flux_dn), face_space),
+            n_faces - half,
+        )
+        LWu_TOA = Fields.level(
+            RRTMGPI.array2field(FT.(atmos_sim.integrator.p.radiation_model.face_lw_flux_up), face_space),
+            n_faces - half,
+        )
+        SWd_TOA = Fields.level(
+            RRTMGPI.array2field(FT.(atmos_sim.integrator.p.radiation_model.face_sw_flux_dn), face_space),
+            n_faces - half,
+        )
+        SWu_TOA = Fields.level(
+            RRTMGPI.array2field(FT.(atmos_sim.integrator.p.radiation_model.face_sw_flux_up), face_space),
+            n_faces - half,
+        )
+
+        return @. -(LWd_TOA + SWd_TOA - LWu_TOA - SWu_TOA) # [W/m^2]
+    else
+        return FT(0)
+    end
+end
+
+get_field(atmos_sim::ClimaAtmosSimulation, ::Val{:energy}) = atmos_sim.integrator.u.c.ρe_tot
+
+get_field(atmos_sim::ClimaAtmosSimulation, ::Val{:water}) = atmos_sim.integrator.u.c.ρq_tot

experiments/AMIP/modular/components/land/bucket_utils.jl

@@ -1,13 +1,7 @@
 
-"""
-    get_bucket_energy(bucket_sim)
-
-Returns the volumetric internal energy of the bucket land model.
-"""
-function get_land_energy(bucket_sim::BucketSimulation, e_per_area)
+function get_field(bucket_sim::BucketSimulation, ::Val{:energy})
     # required by ConservationChecker
-    e_per_area .= zeros(axes(bucket_sim.integrator.u.bucket.W))
-    soil_depth = FT = eltype(bucket_sim.integrator.u.bucket.W)
+    e_per_area = zeros(axes(bucket_sim.integrator.u.bucket.W))
     ClimaCore.Fields.bycolumn(axes(bucket_sim.integrator.u.bucket.T)) do colidx
         e_per_area[colidx] .=
             bucket_sim.model.parameters.ρc_soil .* mean(bucket_sim.integrator.u.bucket.T[colidx]) .*
@@ -20,6 +14,13 @@ function get_land_energy(bucket_sim::BucketSimulation, e_per_area)
     return e_per_area
 end
 
+function get_field(bucket_sim::BucketSimulation, ::Val{:water})
+    ρ_cloud_liq = ClimaLSM.LSMP.ρ_cloud_liq(bucket_sim.model.parameters.earth_param_set)
+    return
+    @. (bucket_sim.integrator.u.bucket.σS + bucket_sim.integrator.u.bucket.W + bucket_sim.integrator.u.bucket.Ws) *
+       ρ_cloud_liq  # kg water / m2
+end
+
 """
     get_land_temp_from_state(land_sim, u)
 Returns the surface temperature of the earth, computed
@@ -105,10 +106,12 @@ function update_field!(sim::BucketSimulation, ::Val{:radiative_energy_flux}, fie
     parent(sim.integrator.p.bucket.R_n) .= parent(field)
 end
 function update_field!(sim::BucketSimulation, ::Val{:liquid_precipitation}, field)
-    parent(sim.integrator.p.bucket.P_liq) .= parent(field)
+    ρ_liq = (LSMP.ρ_cloud_liq(sim.model.parameters.earth_param_set))
+    parent(sim.integrator.p.bucket.P_liq) .= parent(field ./ ρ_liq)
 end
 function update_field!(sim::BucketSimulation, ::Val{:snow_precipitation}, field)
-    parent(sim.integrator.p.bucket.P_snow) .= parent(field)
+    ρ_liq = (LSMP.ρ_cloud_liq(sim.model.parameters.earth_param_set))
+    parent(sim.integrator.p.bucket.P_snow) .= parent(field ./ ρ_liq)
 end
 
 function update_field!(sim::BucketSimulation, ::Val{:air_density}, field)

experiments/AMIP/modular/components/ocean/slab_ocean_init.jl

@@ -157,3 +157,8 @@ Extension of Checkpointer.get_model_state_vector to get the model state.
 function get_model_state_vector(sim::SlabOceanSimulation)
     return sim.integrator.u
 end
+
+get_field(sim::SlabOceanSimulation, ::Val{:energy}) =
+    sim.integrator.p.params.ρ .* sim.integrator.p.params.c .* sim.integrator.u.T_sfc .* sim.integrator.p.params.h
+
+get_field(sim::SlabOceanSimulation, ::Val{:water}) = nothing

experiments/AMIP/modular/components/ocean/slab_seaice_init.jl

@@ -167,3 +167,8 @@ Extension of Checkpointer.get_model_state_vector to get the model state.
 function get_model_state_vector(sim::PrescribedIceSimulation)
     return sim.integrator.u
 end
+
+get_field(sim::PrescribedIceSimulation, ::Val{:energy}) =
+    sim.integrator.p.params.ρ .* sim.integrator.p.params.c .* sim.integrator.u.T_sfc .* sim.integrator.p.params.h
+
+get_field(sim::PrescribedIceSimulation, ::Val{:water}) = nothing

experiments/AMIP/modular/coupler_driver_modular.jl

@@ -387,8 +387,7 @@ if energy_check
         mode_name == "slabplanet" && !CA.is_distributed(ClimaComms.context(boundary_space)),
         "Only non-distributed slabplanet allowable for energy_check"
     )
-    conservation_checks =
-        (; energy = EnergyConservationCheck([], [], [], [], [], []), water = WaterConservationCheck([], [], [], []))
+    conservation_checks = (; energy = EnergyConservationCheck(model_sims), water = WaterConservationCheck(model_sims))
 end
 
 ## coupler simulation
@@ -516,7 +515,7 @@ function solve_coupler!(cs)
         end
 
         ## compute global energy
-        !isnothing(cs.conservation_checks) ? check_conservation!(cs, get_slab_energy, get_land_energy) : nothing
+        !isnothing(cs.conservation_checks) ? check_conservation!(cs) : nothing
 
         ## run component models sequentially for one coupling timestep (Δt_cpl)
         ClimaComms.barrier(comms_ctx)