Checking conservation in the full land model

.buildkite/longruns_gpu/pipeline.yml

@@ -53,36 +53,6 @@ steps:
         env:
           CLIMACOMMS_DEVICE: "CUDA"
 
-      - label: ":sunglasses: California regional simulation"
-        command:
-          - julia --color=yes --project=.buildkite experiments/long_runs/land_region.jl
-        artifact_paths: "california_longrun_gpu/*pdf"
-        agents:
-          slurm_gpus: 1
-          slurm_time: 01:00:00
-        env:
-          CLIMACOMMS_DEVICE: "CUDA"
-
-      - label: "Soil"
-        command:
-          - julia --color=yes --project=.buildkite experiments/long_runs/soil.jl
-        artifact_paths: "soil_longrun_gpu/*pdf"
-        agents:
-          slurm_gpus: 1
-          slurm_time: 3:00:00
-        env:
-          CLIMACOMMS_DEVICE: "CUDA"
-
-      - label: "Global bucket simulation"
-        command:
-          - julia --color=yes --project=.buildkite experiments/long_runs/bucket.jl
-        artifact_paths: "bucket_longrun_gpu/*pdf"
-        agents:
-          slurm_gpus: 1
-          slurm_time: 00:30:00
-        env:
-          CLIMACOMMS_DEVICE: "CUDA"
-
   - group: "Longer runs of Global Land Models"
     if: build.env("LONGER_RUN") != null
     steps:

experiments/long_runs/snowy_land.jl

@@ -343,7 +343,7 @@ function setup_simulation(; greet = false)
     # we simulate from and until the beginning of
     # March so that a full season is included in seasonal metrics.
     start_date = LONGER_RUN ? DateTime("2000-03-01") : DateTime("2008-03-01")
-    stop_date = LONGER_RUN ? DateTime("2019-03-01") : DateTime("2010-03-01")
+    stop_date = LONGER_RUN ? DateTime("2019-03-01") : DateTime("2008-04-01")
     Δt = 450.0
     nelements = (101, 15)
     if greet
@@ -358,8 +358,26 @@ function setup_simulation(; greet = false)
         ClimaLand.ModelSetup.global_domain(FT; comms_ctx = context, nelements)
     params = LP.LandParameters(FT)
     model = setup_model(FT, start_date, stop_date, Δt, domain, params)
-
-    simulation = LandSimulation(FT, start_date, stop_date, Δt, model; outdir)
+    diagnostics = ClimaLand.default_diagnostics(
+        model,
+        start_date;
+        output_writer = ClimaDiagnostics.Writers.NetCDFWriter(
+            domain.space.subsurface,
+            outdir;
+            start_date,
+        ),
+        conservation = true,
+        conservation_period = Day(1),
+    )
+    simulation = LandSimulation(
+        FT,
+        start_date,
+        stop_date,
+        Δt,
+        model;
+        outdir,
+        diagnostics,
+    )
     return simulation
 end
 
@@ -393,20 +411,5 @@ include(
         "figures_function.jl",
     ),
 )
-make_figures(root_path, outdir, short_names)
-
-include("leaderboard/leaderboard.jl")
-diagnostics_folder_path = outdir
-leaderboard_base_path = root_path
-for data_source in ("ERA5", "ILAMB")
-    compute_monthly_leaderboard(
-        leaderboard_base_path,
-        diagnostics_folder_path,
-        data_source,
-    )
-    compute_seasonal_leaderboard(
-        leaderboard_base_path,
-        diagnostics_folder_path,
-        data_source,
-    )
-end
+## Conservation
+check_conservation(root_path, outdir)

src/diagnostics/default_diagnostics.jl

@@ -269,7 +269,9 @@ function default_diagnostics(
     start_date;
     output_writer,
     output_vars = :short,
-    average_period = :monthly,
+    average_period = :daily,
+    conservation = false,
+    conservation_period = Day(10),
 ) where {FT}
 
     define_diagnostics!(land_model)
@@ -388,4 +390,25 @@ function default_diagnostics(
             start_date,
         )
     end
+
+    if conservation
+        additional_diags = ["epa", "epac", "wvpa", "wvpac"]
+        additional_outputs = vcat(
+            map(additional_diags) do short_name
+                output_schedule_func =
+                    conservation_period isa Period ?
+                    EveryCalendarDtSchedule(conservation_period; start_date) : EveryDtSchedule(conservation_period)
+                return ScheduledDiagnostic(
+                    variable = get_diagnostic_variable(short_name),
+                    compute_schedule_func = EveryStepSchedule(),
+                    output_schedule_func = output_schedule_func,
+                    output_writer = output_writer,
+                )
+            end...,
+        )
+    else
+        additional_outputs = []
+    end
+
+    return [default_outputs..., additional_outputs...]
 end

src/diagnostics/land_compute_methods.jl

@@ -55,6 +55,11 @@ end
 @diagnostic_compute "water_volume_per_area_change" EnergyHydrology Y.soil.∫F_vol_liq_water_dt
 @diagnostic_compute "energy_per_area_change" EnergyHydrology Y.soil.∫F_e_dt
 
+@diagnostic_compute "water_volume_per_area" LandModel p.total_water
+@diagnostic_compute "energy_per_area" LandModel p.total_energy
+@diagnostic_compute "water_volume_per_area_change" LandModel Y.∫F_vol_liq_water_dt
+@diagnostic_compute "energy_per_area_change" LandModel Y.∫F_e_dt
+
 
 ### BucketModel ###
 

src/integrated/land.jl

@@ -222,6 +222,8 @@ lsm_aux_vars(m::LandModel) = (
     :ϵ_sfc,
     :α_sfc,
     :α_ground,
+    :total_energy,
+    :total_water,
 )
 
 """
@@ -249,6 +251,8 @@ lsm_aux_types(m::LandModel{FT}) where {FT} = (
     FT,
     FT,
     NamedTuple{(:PAR, :NIR), Tuple{FT, FT}},
+    FT,
+    FT,
 )
 
 """
@@ -276,8 +280,157 @@ lsm_aux_domain_names(m::LandModel) = (
     :surface,
     :surface,
     :surface,
+    :surface,
+    :surface,
 )
 
+function initialize_prognostic(
+    model::LandModel{FT},
+    coords::NamedTuple,
+) where {FT}
+    components = land_components(model)
+    Y_state_list = map(components) do (component)
+        submodel = getproperty(model, component)
+        getproperty(initialize_prognostic(submodel, coords), component)
+    end
+    conservation_vars = (:∫F_vol_liq_water_dt, :∫F_e_dt)
+    conservation_types = (FT, FT)
+    conservation_domain_names = (:surface, :surface)
+    conservation_state_list = initialize_vars(
+        conservation_vars,
+        conservation_types,
+        conservation_domain_names,
+        coords,
+        :conservation_check,
+    )
+    Y = ClimaCore.Fields.FieldVector(;
+        NamedTuple{(components..., conservation_vars...)}((
+            Y_state_list...,
+            conservation_state_list.conservation_check...,
+        ))...,
+    )
+    return Y
+end
+
+function make_update_aux(land::LandModel)
+    components = land_components(land)
+    update_aux_function_list =
+        map(x -> make_update_aux(getproperty(land, x)), components)
+    function update_aux!(p, Y, t)
+        for f! in update_aux_function_list
+            f!(p, Y, t)
+        end
+        sfc_cache = p.scratch1 .* 0
+        ClimaLand.total_liq_water_vol_per_area!(
+            p.total_water,
+            land,
+            Y,
+            p,
+            t,
+            sfc_cache,
+        )
+        ClimaLand.total_energy_per_area!(
+            p.total_energy,
+            land,
+            Y,
+            p,
+            t,
+            sfc_cache,
+        )
+    end
+    return update_aux!
+end
+
+function make_imp_tendency(land::LandModel)
+    components = land_components(land)
+    compute_imp_tendency_list =
+        map(x -> make_compute_imp_tendency(getproperty(land, x)), components)
+    update_aux! = make_update_aux(land)
+    update_boundary_fluxes! = make_update_boundary_fluxes(land)
+    function imp_tendency!(dY, Y, p, t)
+        update_aux!(p, Y, t)
+        update_boundary_fluxes!(p, Y, t)
+        for f! in compute_imp_tendency_list
+            f!(dY, Y, p, t)
+        end
+        atmos = land.snow.boundary_conditions.atmos
+        e_flux_falling_snow =
+            Snow.energy_flux_falling_snow(atmos, p, land.snow.parameters)
+        e_flux_falling_rain =
+            Snow.energy_flux_falling_rain(atmos, p, land.snow.parameters)
+        top_heat_flux = @. lazy(
+            e_flux_falling_snow +
+            e_flux_falling_rain +
+            (
+                p.snow.turbulent_fluxes.lhf +
+                p.snow.turbulent_fluxes.shf +
+                p.snow.R_n
+            ) * p.snow.snow_cover_fraction +
+            (1 - p.snow.snow_cover_fraction) * (
+                p.soil.R_n +
+                p.soil.turbulent_fluxes.lhf +
+                p.soil.turbulent_fluxes.shf
+            ) +
+            p.canopy.turbulent_fluxes.shf +
+            p.canopy.turbulent_fluxes.lhf -
+            p.canopy.radiative_transfer.SW_n -
+            p.canopy.radiative_transfer.LW_n,
+        )
+        bottom_heat_flux = p.soil.bottom_bc.heat
+        #runoff_energy_flux = ;
+        @. dY.∫F_e_dt = -(top_heat_flux - bottom_heat_flux)
+
+        top_water_flux = @. lazy(
+            p.drivers.P_liq +
+            p.drivers.P_snow +
+            p.snow.turbulent_fluxes.vapor_flux * p.snow.snow_cover_fraction +
+            (1 - p.snow.snow_cover_fraction) * (
+                p.soil.turbulent_fluxes.vapor_flux_liq +
+                p.soil.turbulent_fluxes.vapor_flux_ice
+            ) +
+            p.canopy.turbulent_fluxes.transpiration,
+        )
+        runoff_water_flux = @. lazy(p.soil.R_s + p.soil.R_ss)
+        bottom_water_flux = p.soil.bottom_bc.water
+        @. dY.∫F_vol_liq_water_dt =
+            -(top_water_flux + runoff_water_flux - bottom_water_flux)
+
+    end
+    return imp_tendency!
+end
+
+function make_exp_tendency(land::LandModel{FT}) where {FT}
+    components = land_components(land)
+    compute_exp_tendency_list =
+        map(x -> make_compute_exp_tendency(getproperty(land, x)), components)
+    update_aux! = make_update_aux(land)
+    update_boundary_fluxes! = make_update_boundary_fluxes(land)
+    function exp_tendency!(dY, Y, p, t)
+        update_aux!(p, Y, t)
+        # Treat the following terms explicitly!
+
+        # update root extraction
+        update_root_extraction!(p, Y, t, land) # defined in src/integrated/soil_canopy_root_interactions.jl
+
+        # Compute the ground heat flux in place:
+        update_soil_snow_ground_heat_flux!(
+            p,
+            Y,
+            land.soil.parameters,
+            land.snow.parameters,
+            land.soil.domain,
+            FT,
+        )
+        # Treat the following terms explicitly!
+        update_boundary_fluxes!(p, Y, t)
+        for f! in compute_exp_tendency_list
+            f!(dY, Y, p, t)
+        end
+    end
+    return exp_tendency!
+end
+
+
 """
     make_update_boundary_fluxes(
         land::LandModel{FT, MM, SM, RM, SnM},
@@ -340,6 +493,7 @@ function make_update_boundary_fluxes(
             land.soil.domain,
             FT,
         )
+
         #Now update snow boundary conditions, which rely on the ground heat flux
         update_snow_bf!(p, Y, t)
 
@@ -527,7 +681,6 @@ function soil_boundary_fluxes!(
         p.excess_heat_flux +
         p.snow.snow_cover_fraction * p.ground_heat_flux +
         infiltration_energy_flux
-
     return nothing
 end
 
@@ -615,7 +768,7 @@ function snow_boundary_fluxes!(
     e_flux_falling_rain =
         Snow.energy_flux_falling_rain(bc.atmos, p, model.parameters)
 
-    # positive fluxes are TOWARDS atmos, but R_n positive if snow absorbs energy
+    # positive fluxes are TOWARDS atmos
     @. p.snow.total_energy_flux =
         e_flux_falling_snow +
         (

src/shared_utilities/implicit_timestepping.jl

@@ -102,6 +102,8 @@ function FieldMatrixWithSolver(Y::ClimaCore.Fields.FieldVector)
     explicit_vars = (
         @name(soil.∫F_vol_liq_water_dt),
         @name(soil.∫F_e_dt),
+        @name(∫F_vol_liq_water_dt),
+        @name(∫F_e_dt),
         @name(soilco2.C),
         @name(soil.θ_i),
         @name(canopy.hydraulics.ϑ_l),