Adding a new ocean setup

experiments/ClimaEarth/components/atmosphere/climaatmos.jl

@@ -267,11 +267,10 @@ Interfacer.get_field(sim::ClimaAtmosSimulation, ::Val{:air_temperature}) =
         sim.integrator.p.params.thermodynamics_params,
         CC.Fields.level(sim.integrator.p.precomputed.ᶜts, 1),
     )
-Interfacer.get_field(sim::ClimaAtmosSimulation, ::Val{:air_density}) =
-    TD.air_density.(
-        sim.integrator.p.params.thermodynamics_params,
-        CC.Fields.level(sim.integrator.p.precomputed.ᶜts, 1),
-    )
+Interfacer.get_field(sim::ClimaAtmosSimulation, ::Val{:air_density}) = TD.air_density.(
+    sim.integrator.p.params.thermodynamics_params,
+    CC.Fields.level(sim.integrator.p.precomputed.ᶜts, 1),
+)
 # When CO2 is stored as a 1D Array in the tracers cache, we access
 # it from there. Otherwise, we pull a fixed value from ClimaParams.
 Interfacer.get_field(sim::ClimaAtmosSimulation, ::Val{:co2}) =
@@ -292,14 +291,14 @@ function Interfacer.get_field(sim::ClimaAtmosSimulation, ::Val{:diffuse_fraction
     else
         direct_flux_dn = radiation_model.face_sw_direct_flux_dn[1, :]
         FT = eltype(total_flux_dn)
-        diffuse_fraction =
-            clamp.(
-                (
-                    (x, y) -> y > zero(y) ? x / y : zero(y)
-                ).(total_flux_dn .- direct_flux_dn, total_flux_dn),
-                zero(FT),
-                one(FT),
-            )
+        diffuse_fraction = clamp.(
+            ((x, y) -> y > zero(y) ? x / y : zero(y)).(
+                total_flux_dn .- direct_flux_dn,
+                total_flux_dn,
+            ),
+            zero(FT),
+            one(FT),
+        )
     end
     return CC.Fields.array2field(diffuse_fraction, lowest_face_space)
 end
@@ -382,25 +381,23 @@ function Interfacer.update_field!(
     q_sfc_atmos = Interfacer.remap(atmos_surface_space, csf.scalar_temp4)
 
     # Store `ρ_sfc_atmos` in an atmosphere scratch field on the surface space
-    temp_field_surface =
-        FluxCalculator.extrapolate_ρ_to_sfc.(
-            thermo_params,
-            sim.integrator.p.precomputed.sfc_conditions.ts,
-            T_sfc_atmos,
-        )
+    temp_field_surface = FluxCalculator.extrapolate_ρ_to_sfc.(
+        thermo_params,
+        sim.integrator.p.precomputed.sfc_conditions.ts,
+        T_sfc_atmos,
+    )
     ρ_sfc_atmos = temp_field_surface
 
     if sim.integrator.p.atmos.moisture_model isa CA.DryModel
         sim.integrator.p.precomputed.sfc_conditions.ts .=
             TD.PhaseDry_ρT.(thermo_params, ρ_sfc_atmos, T_sfc_atmos)
     else
-        sim.integrator.p.precomputed.sfc_conditions.ts .=
-            TD.PhaseNonEquil_ρTq.(
-                thermo_params,
-                ρ_sfc_atmos,
-                T_sfc_atmos,
-                TD.PhasePartition.(q_sfc_atmos),
-            )
+        sim.integrator.p.precomputed.sfc_conditions.ts .= TD.PhaseNonEquil_ρTq.(
+            thermo_params,
+            ρ_sfc_atmos,
+            T_sfc_atmos,
+            TD.PhasePartition.(q_sfc_atmos),
+        )
     end
 end
 Interfacer.get_field(sim::ClimaAtmosSimulation, ::Val{:height_int}) =
@@ -491,8 +488,7 @@ function FluxCalculator.update_turbulent_fluxes!(sim::ClimaAtmosSimulation, fiel
     Interfacer.remap!(temp_field_surface, F_turb_ρτxz) # F_turb_ρτxz_atmos
     F_turb_ρτyz_atmos = Interfacer.remap(atmos_surface_space, F_turb_ρτyz) # F_turb_ρτyz_atmos
     sim.integrator.p.precomputed.sfc_conditions.ρ_flux_uₕ .= (
-        surface_normal .⊗
-        CA.C12.(
+        surface_normal .⊗ CA.C12.(
             temp_field_surface .* vec_ct12_ct1 .+ F_turb_ρτyz_atmos .* vec_ct12_ct2,
             surface_local_geometry,
         )

experiments/ClimaEarth/components/land/climaland_integrated.jl

@@ -271,20 +271,18 @@ function ClimaLandSimulation(
         # Set initial conditions for the state
         @. Y.soil.ϑ_l = θ_r + (ν - θ_r) / 2
         Y.soil.θ_i .= FT(0.0)
-        ρc_s =
-            CL.Soil.volumetric_heat_capacity.(
-                Y.soil.ϑ_l,
-                Y.soil.θ_i,
-                ρc_ds,
-                earth_param_set,
-            )
-        Y.soil.ρe_int .=
-            CL.Soil.volumetric_internal_energy.(
-                Y.soil.θ_i,
-                ρc_s,
-                orog_adjusted_T,
-                earth_param_set,
-            )
+        ρc_s = CL.Soil.volumetric_heat_capacity.(
+            Y.soil.ϑ_l,
+            Y.soil.θ_i,
+            ρc_ds,
+            earth_param_set,
+        )
+        Y.soil.ρe_int .= CL.Soil.volumetric_internal_energy.(
+            Y.soil.θ_i,
+            ρc_s,
+            orog_adjusted_T,
+            earth_param_set,
+        )
 
         Y.snow.S .= FT(0)
         Y.snow.S_l .= FT(0)
@@ -439,13 +437,12 @@ function Interfacer.update_field!(
     Interfacer.remap!(sim.integrator.p.scratch1, ρ_atmos)
     Interfacer.remap!(sim.integrator.p.scratch2, T_atmos)
     Interfacer.remap!(sim.integrator.p.scratch3, q_atmos)
-    sim.integrator.p.drivers.thermal_state .=
-        TD.PhaseEquil_ρTq.(
-            thermo_params,
-            sim.integrator.p.scratch1,
-            sim.integrator.p.scratch2,
-            sim.integrator.p.scratch3,
-        )
+    sim.integrator.p.drivers.thermal_state .= TD.PhaseEquil_ρTq.(
+        thermo_params,
+        sim.integrator.p.scratch1,
+        sim.integrator.p.scratch2,
+        sim.integrator.p.scratch3,
+    )
 end
 
 function Interfacer.step!(sim::ClimaLandSimulation, t)

experiments/ClimaEarth/components/ocean/clima_seaice.jl

@@ -15,6 +15,26 @@ include("climaocean_helpers.jl")
 # Rename ECCO password env variable to match ClimaOcean.jl
 haskey(ENV, "ECCO_PASSWORD") && (ENV["ECCO_WEBDAV_PASSWORD"] = ENV["ECCO_PASSWORD"])
 
+function OC.set!(model::CSI.SeaIceModel, ::OMIPEvolvedInitialConditions)
+
+    # Dowload the initialization file if not preset
+    if !isfile("omip_initialization.jld2")
+        url = "https://www.dropbox.com/scl/fi/114xafnmoekgc2da3cco7/omip_initialization.jld2?rlkey=e3wnjgwbr2c2zaszysjbptgnh&st=9zagmgux&dl=0"
+        path = "omip_initialization.jld2"
+        download(url, path)
+    end
+
+    file = jldopen("omip_initialization.jld2")
+
+    OC.set!(
+        model,
+        h = file["hi"],
+        ℵ = file["ℵi"]
+    )
+
+    return nothing
+end
+
 """
     ClimaSeaIceSimulation{SIM, A, REMAP, NT, IP}
 
@@ -67,6 +87,7 @@ function ClimaSeaIceSimulation(
     ocean;
     output_dir,
     start_date = nothing,
+    initial_conditions = OMIPEvolvedInitialConditions(),
     coupled_param_dict = CP.create_toml_dict(eltype(ocean.area_fraction)),
     Δt = 5 * 60.0, # 5 minutes
 )
@@ -170,7 +191,7 @@ function sea_ice_simulation(
     top_surface_temperature = OC.Field{OC.Center, OC.Center, Nothing}(grid)
     top_heat_boundary_condition = MeltingConstrainedFluxBalance()
     kᴺ = size(grid, 3)
-    surface_ocean_salinity = OC.interior(ocean.model.tracers.S, :, :, kᴺ:kᴺ)
+    surface_ocean_salinity = OC.interior(ocean.model.tracers.S,:,:,(kᴺ:kᴺ))
     bottom_heat_boundary_condition = IceWaterThermalEquilibrium(surface_ocean_salinity)
 
     ice_thermodynamics = CSI.SlabSeaIceThermodynamics(
@@ -298,8 +319,8 @@ function FluxCalculator.update_turbulent_fluxes!(sim::ClimaSeaIceSimulation, fie
 
     # Update the sea ice only where the concentration is greater than zero.
     si_flux_heat = sim.ice.model.external_heat_fluxes.top[1]
-    OC.interior(si_flux_heat, :, :, 1) .+=
-        (OC.interior(ice_concentration, :, :, 1) .> 0) .* (remapped_F_lh .+ remapped_F_sh)
+    OC.interior(si_flux_heat,:,:,1) .+=
+        (OC.interior(ice_concentration,:,:,1) .> 0) .* (remapped_F_lh .+ remapped_F_sh)
 
     return nothing
 end
@@ -353,8 +374,8 @@ function FieldExchanger.update_sim!(sim::ClimaSeaIceSimulation, csf)
     ϵ = Interfacer.get_field(sim, Val(:emissivity)) # scalar
 
     # Update only where ice concentration is greater than zero.
-    OC.interior(si_flux_heat, :, :, 1) .=
-        (OC.interior(ice_concentration, :, :, 1) .> 0) .*
+    OC.interior(si_flux_heat,:,:,1) .=
+        (OC.interior(ice_concentration,:,:,1) .> 0) .*
         (-(1 .- α) .* remapped_SW_d .- ϵ .* remapped_LW_d)
     return nothing
 end
@@ -433,13 +454,13 @@ function FluxCalculator.ocean_seaice_fluxes!(
     )
 
     oc_flux_T = surface_flux(ocean_sim.ocean.model.tracers.T)
-    OC.interior(oc_flux_T, :, :, 1) .+=
-        OC.interior(ice_concentration, :, :, 1) .* OC.interior(Qi, :, :, 1) .* ρₒ⁻¹ ./ cₒ
+    OC.interior(oc_flux_T,:,:,1) .+=
+        OC.interior(ice_concentration,:,:,1) .* OC.interior(Qi,:,:,1) .* ρₒ⁻¹ ./ cₒ
 
     oc_flux_S = surface_flux(ocean_sim.ocean.model.tracers.S)
-    OC.interior(oc_flux_S, :, :, 1) .+=
-        OC.interior(ice_concentration, :, :, 1) .*
-        OC.interior(ice_sim.ocean_ice_fluxes.salt, :, :, 1)
+    OC.interior(oc_flux_S,:,:,1) .+=
+        OC.interior(ice_concentration,:,:,1) .*
+        OC.interior(ice_sim.ocean_ice_fluxes.salt,:,:,1)
 
     return nothing
 end