Merge pull request #1129 from CliMA/zs/bucket_initcond

add bucket initial condition

Author: szy21

.buildkite/pipeline.yml

@@ -195,6 +195,13 @@ steps:
         artifact_paths: "experiments/ClimaEarth/output/amip_default/artifacts/*"
         agents:
           slurm_mem: 20GB
+      
+      - label: "AMIP: bucket initial condition test"
+        command: "julia --color=yes --project=experiments/ClimaEarth/ experiments/ClimaEarth/run_amip.jl --config_file $CONFIG_PATH/amip_bucket_ic.yml --job_id amip_bucket_ic"
+        artifact_paths: "experiments/ClimaEarth/output/amip_bucket_ic/artifacts/*"
+        agents:
+          slurm_ntasks: 1
+          slurm_mem: 20GB
 
       - label: "AMIP target: albedo from function"
         key: "target_amip_albedo_function"

NEWS.md

@@ -6,6 +6,16 @@ ClimaCoupler.jl Release Notes
 
 ### ClimaEarth features
 
+### Read bucket initial conditions from NetCDF files
+
+Added functionality to allow the bucket initial conditions to be overwritten by interpolated NetCDF datasets.
+To use this feature from the YAML interface, just pass the path of the file to `land_initial_condition`. 
+We expect the file to contain the following variables:
+`W`, for subsurface water storage (2D),
+`Ws`, for surface water content (2D),
+`T`, for soil temperature (3D),
+`S`, for snow water equivalent (2D).
+
 ### Sea-surface temperature and sea ice concentration data can now be automatically downloaded
 
 Sea-surface temperature and sea ice concentration require external files. Now, a

config/ci_configs/amip_bucket_ic.yml

@@ -0,0 +1,18 @@
+apply_limiter: false
+dt: "150secs"
+dt_cpl: "150secs"
+dt_rad: "1hours"
+dt_save_to_sol: "1days"
+dz_bottom: 30
+dz_top: 3000
+h_elem: 4
+land_initial_condition: "experiments/ClimaEarth/input/bucket_ic_august.nc"
+mode_name: "amip"
+moist: "equil"
+precip_model: "0M"
+rad: "gray"
+rayleigh_sponge: true
+t_end: "300secs"
+vert_diff: "true"
+z_elem: 50
+z_stretch: false

experiments/ClimaEarth/cli_options.jl

@@ -155,6 +155,10 @@ function argparse_settings()
         help = "Access land surface albedo information from data file. [`map_static` (default), `function`, `map_temporal`]"
         arg_type = String
         default = "map_static" # to be replaced by land config file, when available
+        "--land_initial_condition"
+        help = "A file path for a NetCDF file (read documentation about requirements)"
+        arg_type = String
+        default = ""
         "--land_temperature_anomaly"
         help = "Type of temperature anomaly for bucket model. [`amip`, `aquaplanet` (default)]"
         arg_type = String

experiments/ClimaEarth/components/land/climaland_bucket.jl

@@ -9,6 +9,7 @@ import ClimaLand as CL
 import ClimaLand.Parameters as LP
 import ClimaDiagnostics as CD
 import ClimaCoupler: Checkpointer, FluxCalculator, Interfacer
+using NCDatasets
 
 ###
 ### Functions required by ClimaCoupler.jl for a SurfaceModelSimulation
@@ -52,6 +53,7 @@ function bucket_init(
     tspan::Tuple{Float64, Float64},
     config::String,
     albedo_type::String,
+    land_initial_condition::String,
     land_temperature_anomaly::String,
     output_dir::String;
     space,
@@ -138,6 +140,53 @@ function bucket_init(
     Y.bucket.Ws .= 0.0
     Y.bucket.σS .= 0.0
 
+    # Overwrite initial conditions with interpolated values from a netcdf file using
+    # the `SpaceVaryingInputs` tool. We expect the file to contain the following variables:
+    # - `W`, for subsurface water storage (2D),
+    # - `Ws`, for surface water content (2D),
+    # - `T`, for soil temperature (3D),
+    # - `S`, for snow water equivalent (2D).
+
+    if !isempty(land_initial_condition)
+        ds = NCDataset(land_initial_condition)
+        has_all_variables = all(key -> haskey(ds, key), ["W", "Ws", "T", "S"])
+        @assert has_all_variables "The land iniital condition file is expected to contain the variables W, Ws, T, and S (read documentation about requirements)."
+        close(ds)
+
+        surface_space = domain.space.surface
+        subsurface_space = domain.space.subsurface
+        regridder_type = :InterpolationsRegridder
+        extrapolation_bc = (Interpolations.Periodic(), Interpolations.Flat(), Interpolations.Flat())
+        Y.bucket.W .= SpaceVaryingInput(
+            land_initial_condition,
+            "W",
+            surface_space;
+            regridder_type,
+            regridder_kwargs = (; extrapolation_bc,),
+        )
+        Y.bucket.Ws .= SpaceVaryingInput(
+            land_initial_condition,
+            "Ws",
+            surface_space;
+            regridder_type,
+            regridder_kwargs = (; extrapolation_bc,),
+        )
+        Y.bucket.T .= SpaceVaryingInput(
+            land_initial_condition,
+            "T",
+            subsurface_space;
+            regridder_type,
+            regridder_kwargs = (; extrapolation_bc,),
+        )
+        Y.bucket.σS .= SpaceVaryingInput(
+            land_initial_condition,
+            "S",
+            surface_space;
+            regridder_type,
+            regridder_kwargs = (; extrapolation_bc,),
+        )
+    end
+
     # Set initial aux variable values
     set_initial_cache! = CL.make_set_initial_cache(model)
     set_initial_cache!(p, Y, tspan[1])

experiments/ClimaEarth/input/bucket_ic_august.nc

@@ -118,6 +118,7 @@ add_extra_diagnostics!(config_dict)
     turb_flux_partition,
     land_domain_type,
     land_albedo_type,
+    land_initial_condition,
     land_temperature_anomaly,
     energy_check,
     conservation_softfail,
@@ -253,6 +254,7 @@ if mode_name == "amip"
         tspan,
         land_domain_type,
         land_albedo_type,
+        land_initial_condition,
         land_temperature_anomaly,
         land_output_dir;
         dt = component_dt_dict["dt_land"],
@@ -351,6 +353,7 @@ elseif mode_name in ("slabplanet", "slabplanet_aqua", "slabplanet_terra")
         tspan,
         land_domain_type,
         land_albedo_type,
+        land_initial_condition,
         land_temperature_anomaly,
         land_output_dir;
         dt = component_dt_dict["dt_land"],
@@ -401,6 +404,7 @@ elseif mode_name == "slabplanet_eisenman"
         tspan,
         land_domain_type,
         land_albedo_type,
+        land_initial_condition,
         land_temperature_anomaly,
         land_output_dir;
         dt = component_dt_dict["dt_land"],

experiments/ClimaEarth/run_amip.jl

@@ -156,11 +156,6 @@ atmos_config_object.toml_dict["max_area_limiter_scale"]["value"] = 0
 
 comms_ctx = Utilities.get_comms_context(Dict("device" => "auto"))
 
-#=
-## Data File Paths
-=#
-land_mask_data = joinpath(@clima_artifact("landsea_mask_60arcseconds", comms_ctx), "landsea_mask.nc")
-
 #=
 ## Component Model Initialization
 =#
@@ -187,6 +182,10 @@ thermo_params = get_thermo_params(atmos_sim)
 ## init a 2D boundary space at the surface
 boundary_space = CC.Spaces.horizontal_space(atmos_sim.domain.face_space) # TODO: specify this in the coupler and pass it to all component models #665
 
+# Land initial condition
+# Use the default land initial condition (not reading from a file)
+land_initial_condition = ""
+
 #=
 ### Land-sea Fraction
 This is a static field that contains the area fraction of land and sea, ranging from 0 to 1. If applicable, sea ice is included in the sea fraction. at this stage.
@@ -206,6 +205,7 @@ land_sim = bucket_init(
     tspan,
     "sphere",
     "map_static",
+    land_initial_condition,
     "aquaplanet",
     land_output_dir;
     dt = Δt_cpl,

experiments/ClimaEarth/run_cloudy_slabplanet.jl

@@ -85,6 +85,7 @@ function get_coupler_args(config_dict::Dict)
     # ClimaLand-specific information
     land_domain_type = config_dict["land_domain_type"]
     land_albedo_type = config_dict["land_albedo_type"]
+    land_initial_condition = config_dict["land_initial_condition"]
     land_temperature_anomaly = config_dict["land_temperature_anomaly"]
     use_land_diagnostics = config_dict["use_land_diagnostics"]
 
@@ -116,6 +117,7 @@ function get_coupler_args(config_dict::Dict)
         plot_diagnostics,
         land_domain_type,
         land_albedo_type,
+        land_initial_condition,
         land_temperature_anomaly,
         use_land_diagnostics,
     )