diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml
index 32ac7be35e..821f87dd64 100644
--- a/.buildkite/pipeline.yml
+++ b/.buildkite/pipeline.yml
@@ -38,19 +38,19 @@ steps:
       - "julia --project=experiments/ClimaEarth/ -e 'using Pkg; Pkg.precompile()'"
       - "julia --project=experiments/ClimaEarth/ -e 'using Pkg; Pkg.status()'"
 
-      - echo "--- Instantiate ClimaCore experiments env"
-      - "julia --project=experiments/ClimaCore/ -e 'using Pkg; Pkg.develop(path=\".\")'"
-      - "julia --project=experiments/ClimaCore/ -e 'using Pkg; Pkg.instantiate(;verbose=true)'"
-      - "julia --project=experiments/ClimaCore/ -e 'using Pkg; Pkg.add(\"MPI\")'"
-      - "julia --project=experiments/ClimaCore/ -e 'using Pkg; Pkg.precompile()'"
-      - "julia --project=experiments/ClimaCore/ -e 'using Pkg; Pkg.status()'"
-
-      - echo "--- Instantiate test env"
-      - "julia --project=test/ -e 'using Pkg; Pkg.develop(path=\".\")'"
-      - "julia --project=test/ -e 'using Pkg; Pkg.instantiate(;verbose=true)'"
-      - "julia --project=test/ -e 'using Pkg; Pkg.add(\"MPI\")'"
-      - "julia --project=test/ -e 'using Pkg; Pkg.precompile()'"
-      - "julia --project=test/ -e 'using Pkg; Pkg.status()'"
+      # - echo "--- Instantiate ClimaCore experiments env"
+      # - "julia --project=experiments/ClimaCore/ -e 'using Pkg; Pkg.develop(path=\".\")'"
+      # - "julia --project=experiments/ClimaCore/ -e 'using Pkg; Pkg.instantiate(;verbose=true)'"
+      # - "julia --project=experiments/ClimaCore/ -e 'using Pkg; Pkg.add(\"MPI\")'"
+      # - "julia --project=experiments/ClimaCore/ -e 'using Pkg; Pkg.precompile()'"
+      # - "julia --project=experiments/ClimaCore/ -e 'using Pkg; Pkg.status()'"
+
+      # - echo "--- Instantiate test env"
+      # - "julia --project=test/ -e 'using Pkg; Pkg.develop(path=\".\")'"
+      # - "julia --project=test/ -e 'using Pkg; Pkg.instantiate(;verbose=true)'"
+      # - "julia --project=test/ -e 'using Pkg; Pkg.add(\"MPI\")'"
+      # - "julia --project=test/ -e 'using Pkg; Pkg.precompile()'"
+      # - "julia --project=test/ -e 'using Pkg; Pkg.status()'"
 
     concurrency: 1
     concurrency_group: 'depot/climacoupler-ci'
@@ -63,311 +63,329 @@ steps:
 
   - wait
 
-  - group: "Unit Tests"
-    steps:
-
-      - label: "MPI Utilities unit tests"
-        key: "utilities_mpi_tests"
-        command: "srun julia --color=yes --project=test/ test/utilities_tests.jl"
-        timeout_in_minutes: 5
-        env:
-          CLIMACOMMS_CONTEXT: "MPI"
-        agents:
-          slurm_ntasks: 2
-          slurm_mem: 16GB
-
-      - label: "MPI Interfacer unit tests"
-        key: "interfacer_mpi_tests"
-        command: "srun julia --color=yes --project=test/ test/interfacer_tests.jl"
-        timeout_in_minutes: 5
-        env:
-          CLIMACOMMS_CONTEXT: "MPI"
-        agents:
-          slurm_ntasks: 2
-          slurm_mem: 16GB
-
-  - group: "GPU: unit tests"
-    steps:
-      - label: "GPU runtests"
-        command: "julia --color=yes --project=test/ test/runtests.jl"
-        timeout_in_minutes: 10
-        env:
-          CLIMACOMMS_DEVICE: "CUDA"
-        agents:
-          slurm_ntasks: 1
-          slurm_gres: "gpu:1"
-          slurm_mem: 24GB
-
-  - group: "ClimaEarth tests"
-    steps:
-      - label: "ClimaEarth runtests"
-        command: "julia --color=yes --project=experiments/ClimaEarth/ experiments/ClimaEarth/test/runtests.jl"
-        agents:
-          slurm_mem: 16GB
-
-      - label: "MPI restarts"
-        key: "mpi_restarts"
-        command: "srun julia --color=yes --project=experiments/ClimaEarth/ experiments/ClimaEarth/test/restart.jl"
-        env:
-          CLIMACOMMS_CONTEXT: "MPI"
-        timeout_in_minutes: 50
-        soft_fail:
-          - exit_status: -1
-          - exit_status: 255
-        agents:
-          slurm_ntasks: 2
-          slurm_mem: 32GB
-
-      - label: "GPU restarts"
-        command: "julia --color=yes --project=experiments/ClimaEarth/ experiments/ClimaEarth/test/restart.jl"
-        env:
-          CLIMACOMMS_DEVICE: "CUDA"
-        agents:
-          slurm_ntasks: 1
-          slurm_gres: "gpu:1"
-          slurm_mem: 32GB
-
-  - group: "Integration Tests"
-    steps:
-      # SLABPLANET EXPERIMENTS
-
-      # Slabplanet default:
-      # - this is the most lightweight example with conservation and visual checks, with CLI specification as follows
-      #   - numerics: dt = dt_cpl = 200s, nelem = 4
-      #   - physics: bulk aerodynamic surface fluxes, gray radiation, idealized insolation, equil moisture model, 0-moment microphysics
-      #   - input data: monotonous remapping (land mask, SST, SIC)
-      #   - slurm: unthreaded, 1 ntask
-      #   - diagnostics: check and plot energy conservation, output plots after 9 days
-      - label: "Slabplanet: default"
-        key: "slabplanet_default"
-        command: "julia --color=yes --project=experiments/ClimaEarth/ experiments/ClimaEarth/run_amip.jl --config_file $CONFIG_PATH/slabplanet_default.yml --job_id slabplanet_default"
-        artifact_paths: "experiments/ClimaEarth/output/slabplanet_default/artifacts/*"
-        agents:
-          slurm_mem: 20GB
-
-      - label: "Slabplanet: dry, no radiation, fixed ocean T"
-        key: "slabplanet_dry_norad"
-        command: "julia --color=yes --project=experiments/ClimaEarth/ experiments/ClimaEarth/run_amip.jl --config_file $CONFIG_PATH/slabplanet_dry_norad.yml --job_id slabplanet_dry_norad"
-        artifact_paths: "experiments/ClimaEarth/output/slabplanet_dry_norad/artifacts/*"
-        agents:
-          slurm_mem: 20GB
-
-      - label: "Slabplanet: extra atmos diagnostics"
-        key: "slabplanet_atmos_diags"
-        command: "julia --color=yes --project=experiments/ClimaEarth/ experiments/ClimaEarth/run_amip.jl --config_file $CONFIG_PATH/slabplanet_atmos_diags.yml --job_id slabplanet_atmos_diags"
-        artifact_paths: "experiments/ClimaEarth/output/slabplanet_atmos_diags/artifacts/*"
-        agents:
-          slurm_mem: 20GB
-
-      - label: "Slabplanet terra: atmos and bucket"
-        key: "slabplanet_terra"
-        command: "julia --color=yes --project=experiments/ClimaEarth/ experiments/ClimaEarth/run_amip.jl --config_file $CONFIG_PATH/slabplanet_terra.yml --job_id slabplanet_terra"
-        artifact_paths: "experiments/ClimaEarth/output/slabplanet_terra/artifacts/*"
-        agents:
-          slurm_mem: 20GB
-
-      - label: "Slabplanet aqua: atmos and slab ocean"
-        key: "slabplanet_aqua"
-        command: "julia --color=yes --project=experiments/ClimaEarth/ experiments/ClimaEarth/run_amip.jl --config_file $CONFIG_PATH/slabplanet_aqua.yml --job_id slabplanet_aqua"
-        artifact_paths: "experiments/ClimaEarth/output/slabplanet_aqua/artifacts/*"
-        agents:
-          slurm_mem: 20GB
-
-      # AMIP EXPERIMENTS
-
-      # Test default behavior with no config file or job ID provided
-      - label: "AMIP: default"
-        key: "amip_default"
-        command: "julia --color=yes --project=experiments/ClimaEarth/ experiments/ClimaEarth/run_amip.jl"
-        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: integrated land non-spun up initial condition test"
-        command: "julia --color=yes --project=experiments/ClimaEarth/ experiments/ClimaEarth/run_amip.jl --config_file $CONFIG_PATH/amip_land_ic.yml --job_id amip_land_ic"
-        artifact_paths: "experiments/ClimaEarth/output/amip_land_ic/artifacts/*"
-        agents:
-          slurm_ntasks: 1
-          slurm_mem: 20GB
-
-      - label: "AMIP - Float64 + hourly checkpoint"
-        key: "amip"
-        command: "julia --color=yes --project=experiments/ClimaEarth/ experiments/ClimaEarth/run_amip.jl --config_file $CONFIG_PATH/amip_coarse_ft64_hourly_checkpoints.yml --job_id amip_coarse_ft64_hourly_checkpoints"
-        artifact_paths: "experiments/ClimaEarth/output/amip_coarse_ft64_hourly_checkpoints/artifacts/*"
-        env:
-          FLAME_PLOT: ""
-          BUILD_HISTORY_HANDLE: ""
-        agents:
-          slurm_ntasks: 1
-          slurm_mem: 20GB
-
-      - label: "AMIP - Component dts test"
-        command: "julia --color=yes --project=experiments/ClimaEarth/ experiments/ClimaEarth/run_amip.jl --config_file $CONFIG_PATH/amip_component_dts.yml --job_id target_amip_component_dts"
-        artifact_paths: "experiments/ClimaEarth/output/target_amip_component_dts/artifacts/*"
-        agents:
-          slurm_ntasks: 1
-          slurm_mem: 20GB
-
-      - label: "MPI AMIP"
-        command: "srun julia --color=yes --project=experiments/ClimaEarth/ experiments/ClimaEarth/run_amip.jl --config_file $CONFIG_PATH/amip_coarse_mpi.yml --job_id amip_coarse_mpi"
-        artifact_paths: "experiments/ClimaEarth/output/amip_coarse_mpi/artifacts/*"
-        timeout_in_minutes: 30
-        env:
-          CLIMACOMMS_CONTEXT: "MPI"
-        agents:
-          slurm_ntasks: 4
-          slurm_mem_per_cpu: 12GB
-
-      # short high-res performance test
-      - label: "Unthreaded AMIP FINE" # also reported by longruns with a flame graph
-        key: "unthreaded_amip_fine"
-        command: "julia --color=yes --project=experiments/ClimaEarth/ experiments/ClimaEarth/run_amip.jl --config_file $CONFIG_PATH/amip_n1_shortrun.yml --job_id target_amip_n1_shortrun"
-        artifact_paths: "experiments/ClimaEarth/output/target_amip_n1_shortrun/artifacts/*"
-        env:
-          BUILD_HISTORY_HANDLE: ""
-        agents:
-          slurm_mem: 20GB
-
-      # CLIMACORE EXPERIMENTS
-
-      - label: "sea_breeze"
-        command: "julia --color=yes --project=experiments/ClimaCore experiments/ClimaCore/sea_breeze/run.jl"
-        artifact_paths: "experiments/ClimaCore/sea_breeze/output/*"
-        agents:
-          slurm_mem: 20GB
+  # - group: "Unit Tests"
+  #   steps:
+
+  #     - label: "MPI Utilities unit tests"
+  #       key: "utilities_mpi_tests"
+  #       command: "srun julia --color=yes --project=test/ test/utilities_tests.jl"
+  #       timeout_in_minutes: 5
+  #       env:
+  #         CLIMACOMMS_CONTEXT: "MPI"
+  #       agents:
+  #         slurm_ntasks: 2
+  #         slurm_mem: 16GB
+
+  #     - label: "MPI Interfacer unit tests"
+  #       key: "interfacer_mpi_tests"
+  #       command: "srun julia --color=yes --project=test/ test/interfacer_tests.jl"
+  #       timeout_in_minutes: 5
+  #       env:
+  #         CLIMACOMMS_CONTEXT: "MPI"
+  #       agents:
+  #         slurm_ntasks: 2
+  #         slurm_mem: 16GB
+
+  # - group: "GPU: unit tests"
+  #   steps:
+  #     - label: "GPU runtests"
+  #       command: "julia --color=yes --project=test/ test/runtests.jl"
+  #       timeout_in_minutes: 10
+  #       env:
+  #         CLIMACOMMS_DEVICE: "CUDA"
+  #       agents:
+  #         slurm_ntasks: 1
+  #         slurm_gres: "gpu:1"
+  #         slurm_mem: 24GB
+
+  # - group: "ClimaEarth tests"
+  #   steps:
+  #     - label: "ClimaEarth runtests"
+  #       command: "julia --color=yes --project=experiments/ClimaEarth/ experiments/ClimaEarth/test/runtests.jl"
+  #       agents:
+  #         slurm_mem: 16GB
+
+  #     - label: "MPI restarts"
+  #       key: "mpi_restarts"
+  #       command: "srun julia --color=yes --project=experiments/ClimaEarth/ experiments/ClimaEarth/test/restart.jl"
+  #       env:
+  #         CLIMACOMMS_CONTEXT: "MPI"
+  #       timeout_in_minutes: 50
+  #       soft_fail:
+  #         - exit_status: -1
+  #         - exit_status: 255
+  #       agents:
+  #         slurm_ntasks: 2
+  #         slurm_mem: 32GB
+
+  #     - label: "GPU restarts"
+  #       command: "julia --color=yes --project=experiments/ClimaEarth/ experiments/ClimaEarth/test/restart.jl"
+  #       env:
+  #         CLIMACOMMS_DEVICE: "CUDA"
+  #       agents:
+  #         slurm_ntasks: 1
+  #         slurm_gres: "gpu:1"
+  #         slurm_mem: 32GB
+
+  # - group: "Integration Tests"
+  #   steps:
+  #     # SLABPLANET EXPERIMENTS
+
+  #     # Slabplanet default:
+  #     # - this is the most lightweight example with conservation and visual checks, with CLI specification as follows
+  #     #   - numerics: dt = dt_cpl = 200s, nelem = 4
+  #     #   - physics: bulk aerodynamic surface fluxes, gray radiation, idealized insolation, equil moisture model, 0-moment microphysics
+  #     #   - input data: monotonous remapping (land mask, SST, SIC)
+  #     #   - slurm: unthreaded, 1 ntask
+  #     #   - diagnostics: check and plot energy conservation, output plots after 9 days
+  #     - label: "Slabplanet: default"
+  #       key: "slabplanet_default"
+  #       command: "julia --color=yes --project=experiments/ClimaEarth/ experiments/ClimaEarth/run_amip.jl --config_file $CONFIG_PATH/slabplanet_default.yml --job_id slabplanet_default"
+  #       artifact_paths: "experiments/ClimaEarth/output/slabplanet_default/artifacts/*"
+  #       agents:
+  #         slurm_mem: 20GB
+
+  #     - label: "Slabplanet: dry, no radiation, fixed ocean T"
+  #       key: "slabplanet_dry_norad"
+  #       command: "julia --color=yes --project=experiments/ClimaEarth/ experiments/ClimaEarth/run_amip.jl --config_file $CONFIG_PATH/slabplanet_dry_norad.yml --job_id slabplanet_dry_norad"
+  #       artifact_paths: "experiments/ClimaEarth/output/slabplanet_dry_norad/artifacts/*"
+  #       agents:
+  #         slurm_mem: 20GB
+
+  #     - label: "Slabplanet: extra atmos diagnostics"
+  #       key: "slabplanet_atmos_diags"
+  #       command: "julia --color=yes --project=experiments/ClimaEarth/ experiments/ClimaEarth/run_amip.jl --config_file $CONFIG_PATH/slabplanet_atmos_diags.yml --job_id slabplanet_atmos_diags"
+  #       artifact_paths: "experiments/ClimaEarth/output/slabplanet_atmos_diags/artifacts/*"
+  #       agents:
+  #         slurm_mem: 20GB
+
+  #     - label: "Slabplanet terra: atmos and bucket"
+  #       key: "slabplanet_terra"
+  #       command: "julia --color=yes --project=experiments/ClimaEarth/ experiments/ClimaEarth/run_amip.jl --config_file $CONFIG_PATH/slabplanet_terra.yml --job_id slabplanet_terra"
+  #       artifact_paths: "experiments/ClimaEarth/output/slabplanet_terra/artifacts/*"
+  #       agents:
+  #         slurm_mem: 20GB
+
+  #     - label: "Slabplanet aqua: atmos and slab ocean"
+  #       key: "slabplanet_aqua"
+  #       command: "julia --color=yes --project=experiments/ClimaEarth/ experiments/ClimaEarth/run_amip.jl --config_file $CONFIG_PATH/slabplanet_aqua.yml --job_id slabplanet_aqua"
+  #       artifact_paths: "experiments/ClimaEarth/output/slabplanet_aqua/artifacts/*"
+  #       agents:
+  #         slurm_mem: 20GB
+
+  #     # AMIP EXPERIMENTS
+
+  #     # Test default behavior with no config file or job ID provided
+  #     - label: "AMIP: default"
+  #       key: "amip_default"
+  #       command: "julia --color=yes --project=experiments/ClimaEarth/ experiments/ClimaEarth/run_amip.jl"
+  #       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: integrated land non-spun up initial condition test"
+  #       command: "julia --color=yes --project=experiments/ClimaEarth/ experiments/ClimaEarth/run_amip.jl --config_file $CONFIG_PATH/amip_land_ic.yml --job_id amip_land_ic"
+  #       artifact_paths: "experiments/ClimaEarth/output/amip_land_ic/artifacts/*"
+  #       agents:
+  #         slurm_ntasks: 1
+  #         slurm_mem: 20GB
+
+  #     - label: "AMIP - Float64 + hourly checkpoint"
+  #       key: "amip"
+  #       command: "julia --color=yes --project=experiments/ClimaEarth/ experiments/ClimaEarth/run_amip.jl --config_file $CONFIG_PATH/amip_coarse_ft64_hourly_checkpoints.yml --job_id amip_coarse_ft64_hourly_checkpoints"
+  #       artifact_paths: "experiments/ClimaEarth/output/amip_coarse_ft64_hourly_checkpoints/artifacts/*"
+  #       env:
+  #         FLAME_PLOT: ""
+  #         BUILD_HISTORY_HANDLE: ""
+  #       agents:
+  #         slurm_ntasks: 1
+  #         slurm_mem: 20GB
+
+
+  #     - label: "AMIP - Component dts test"
+  #       command: "julia --color=yes --project=experiments/ClimaEarth/ experiments/ClimaEarth/run_amip.jl --config_file $CONFIG_PATH/amip_component_dts.yml --job_id target_amip_component_dts"
+  #       artifact_paths: "experiments/ClimaEarth/output/target_amip_component_dts/artifacts/*"
+  #       agents:
+  #         slurm_ntasks: 1
+  #         slurm_mem: 20GB
+
+  #     - label: "MPI AMIP"
+  #       command: "srun julia --color=yes --project=experiments/ClimaEarth/ experiments/ClimaEarth/run_amip.jl --config_file $CONFIG_PATH/amip_coarse_mpi.yml --job_id amip_coarse_mpi"
+  #       artifact_paths: "experiments/ClimaEarth/output/amip_coarse_mpi/artifacts/*"
+  #       timeout_in_minutes: 30
+  #       env:
+  #         CLIMACOMMS_CONTEXT: "MPI"
+  #       agents:
+  #         slurm_ntasks: 4
+  #         slurm_mem_per_cpu: 12GB
+
+  #     # short high-res performance test
+  #     - label: "Unthreaded AMIP FINE" # also reported by longruns with a flame graph
+  #       key: "unthreaded_amip_fine"
+  #       command: "julia --color=yes --project=experiments/ClimaEarth/ experiments/ClimaEarth/run_amip.jl --config_file $CONFIG_PATH/amip_n1_shortrun.yml --job_id target_amip_n1_shortrun"
+  #       artifact_paths: "experiments/ClimaEarth/output/target_amip_n1_shortrun/artifacts/*"
+  #       env:
+  #         BUILD_HISTORY_HANDLE: ""
+  #       agents:
+  #         slurm_mem: 20GB
+
+  #     # CLIMACORE EXPERIMENTS
+
+  #     - label: "sea_breeze"
+  #       command: "julia --color=yes --project=experiments/ClimaCore experiments/ClimaCore/sea_breeze/run.jl"
+  #       artifact_paths: "experiments/ClimaCore/sea_breeze/output/*"
+  #       agents:
+  #         slurm_mem: 20GB
+
+  #     - label: "heat-diffusion"
+  #       command: "julia --color=yes --project=experiments/ClimaCore/ experiments/ClimaCore/heat-diffusion/run.jl"
+  #       artifact_paths: "experiments/ClimaCore/output/heat-diffusion/artifacts/*"
+  #       agents:
+  #         slurm_mem: 20GB
+
+  # - group: "GPU integration tests"
+  #   steps:
+  #     # GPU RUNS: slabplanet
+  #     - label: "GPU Slabplanet: albedo from function"
+  #       key: "gpu_slabplanet_albedo_function"
+  #       command: "julia --color=yes --project=experiments/ClimaEarth/ experiments/ClimaEarth/run_amip.jl --config_file $CONFIG_PATH/slabplanet_albedo_function.yml --job_id gpu_slabplanet_albedo_function"
+  #       artifact_paths: "experiments/ClimaEarth/output/gpu_slabplanet_albedo_function/artifacts/*"
+  #       env:
+  #         CLIMACOMMS_DEVICE: "CUDA"
+  #       agents:
+  #         slurm_mem: 20GB
+  #         slurm_gpus: 1
+
+  #     # GPU RUNS: AMIP
+  #     - label: "GPU AMIP: ED only + integrated land"
+  #       key: "gpu_amip_edonly_integrated_land"
+  #       command: "julia --color=yes --project=experiments/ClimaEarth/ experiments/ClimaEarth/run_amip.jl --config_file $CONFIG_PATH/amip_edonly_integrated_land.yml --job_id gpu_amip_edonly_integrated_land"
+  #       artifact_paths: "experiments/ClimaEarth/output/gpu_amip_edonly_integrated_land/artifacts/*"
+  #       env:
+  #         CLIMACOMMS_DEVICE: "CUDA"
+  #       agents:
+  #         slurm_mem: 20GB
+  #         slurm_gpus: 1
+
+  #     - label: "GPU AMIP: ED only + bucket"
+  #       key: "gpu_amip_edonly_bucket"
+  #       command: "julia --color=yes --project=experiments/ClimaEarth/ experiments/ClimaEarth/run_amip.jl --config_file $CONFIG_PATH/amip_edonly_bucket.yml --job_id gpu_amip_edonly_bucket"
+  #       artifact_paths: "experiments/ClimaEarth/output/gpu_amip_edonly_bucket/artifacts/*"
+  #       env:
+  #         CLIMACOMMS_DEVICE: "CUDA"
+  #       agents:
+  #         slurm_mem: 20GB
+  #         slurm_gpus: 1
+
+  #     - label: "GPU AMIP: diag. EDMF + integrated land"
+  #       key: "gpu_amip_diagedmf_integrated_land"
+  #       command: "julia --color=yes --project=experiments/ClimaEarth/ experiments/ClimaEarth/run_amip.jl --config_file $CONFIG_PATH/amip_diagedmf_integrated_land.yml --job_id gpu_amip_diagedmf_integrated_land"
+  #       artifact_paths: "experiments/ClimaEarth/output/gpu_amip_diagedmf_integrated_land/artifacts/*"
+  #       env:
+  #         CLIMACOMMS_DEVICE: "CUDA"
+  #       agents:
+  #         slurm_mem: 20GB
+  #         slurm_gpus: 1
+
+  #     - label: "GPU AMIP: diag. EDMF + bucket"
+  #       key: "gpu_amip_diagedmf_bucket"
+  #       command: "julia --color=yes --project=experiments/ClimaEarth/ experiments/ClimaEarth/run_amip.jl --config_file $CONFIG_PATH/amip_diagedmf_bucket.yml --job_id gpu_amip_diagedmf_bucket"
+  #       artifact_paths: "experiments/ClimaEarth/output/gpu_amip_diagedmf_bucket/artifacts/*"
+  #       env:
+  #         CLIMACOMMS_DEVICE: "CUDA"
+  #       agents:
+  #         slurm_mem: 20GB
+  #         slurm_gpus: 1
+
+  #     - label: "GPU AMIP test: albedo from function"
+  #       key: "gpu_amip_albedo_function"
+  #       command: "julia --color=yes --project=experiments/ClimaEarth/ experiments/ClimaEarth/run_amip.jl --config_file $CONFIG_PATH/amip_albedo_function.yml --job_id gpu_amip_albedo_function"
+  #       artifact_paths: "experiments/ClimaEarth/output/gpu_amip_albedo_function/artifacts/*"
+  #       env:
+  #         CLIMACOMMS_DEVICE: "CUDA"
+  #       agents:
+  #         slurm_mem: 20GB
+  #         slurm_gpus: 1
+
+  #     - label: "GPU AMIP: albedo from temporal map + 0M"
+  #       key: "gpu_amip_albedo_temporal_map"
+  #       command: "julia --color=yes --project=experiments/ClimaEarth/ experiments/ClimaEarth/run_amip.jl --config_file $CONFIG_PATH/amip_albedo_temporal_map.yml --job_id gpu_amip_albedo_temporal_map"
+  #       artifact_paths: "experiments/ClimaEarth/output/gpu_amip_albedo_temporal_map/artifacts/*"
+  #       env:
+  #         CLIMACOMMS_DEVICE: "CUDA"
+  #       agents:
+  #         slurm_mem: 20GB
+  #         slurm_gpus: 1
+
+  #     - label: "GPU AMIP: albedo from temporal map + 1M"
+  #       key: "gpu_amip_albedo_temporal_map_1M"
+  #       command: "julia --color=yes --project=experiments/ClimaEarth/ experiments/ClimaEarth/run_amip.jl --config_file $CONFIG_PATH/amip_albedo_temporal_map_1M.yml --job_id gpu_amip_albedo_temporal_map_1M"
+  #       artifact_paths: "experiments/ClimaEarth/output/gpu_amip_albedo_temporal_map_1M/artifacts/*"
+  #       env:
+  #         CLIMACOMMS_DEVICE: "CUDA"
+  #       agents:
+  #         slurm_mem: 20GB
+  #         slurm_gpus: 1
 
-      - label: "heat-diffusion"
-        command: "julia --color=yes --project=experiments/ClimaCore/ experiments/ClimaCore/heat-diffusion/run.jl"
-        artifact_paths: "experiments/ClimaCore/output/heat-diffusion/artifacts/*"
-        agents:
-          slurm_mem: 20GB
-
-  - group: "GPU integration tests"
+  - group: "CMIP"
     steps:
-      # GPU RUNS: slabplanet
-      - label: "GPU Slabplanet: albedo from function"
-        key: "gpu_slabplanet_albedo_function"
-        command: "julia --color=yes --project=experiments/ClimaEarth/ experiments/ClimaEarth/run_amip.jl --config_file $CONFIG_PATH/slabplanet_albedo_function.yml --job_id gpu_slabplanet_albedo_function"
-        artifact_paths: "experiments/ClimaEarth/output/gpu_slabplanet_albedo_function/artifacts/*"
-        env:
-          CLIMACOMMS_DEVICE: "CUDA"
-        agents:
-          slurm_mem: 20GB
-          slurm_gpus: 1
 
-      # GPU RUNS: AMIP
-      - label: "GPU AMIP: ED only + integrated land"
-        key: "gpu_amip_edonly_integrated_land"
-        command: "julia --color=yes --project=experiments/ClimaEarth/ experiments/ClimaEarth/run_amip.jl --config_file $CONFIG_PATH/amip_edonly_integrated_land.yml --job_id gpu_amip_edonly_integrated_land"
-        artifact_paths: "experiments/ClimaEarth/output/gpu_amip_edonly_integrated_land/artifacts/*"
+      - label: "GPU CMIP: ClimaAtmos + bucket land + Oceananigans + PrescribedSeaIce"
+        key: "gpu_cmip_oceananigans_prescrseaice"
+        command: "julia --color=yes --project=experiments/ClimaEarth/ experiments/ClimaEarth/run_amip.jl --config_file $CONFIG_PATH/cmip_oceananigans_prescrseaice.yml --job_id gpu_cmip_oceananigans_prescrseaice"
+        artifact_paths: "experiments/ClimaEarth/output/gpu_cmip_oceananigans_prescrseaice/artifacts/*"
         env:
           CLIMACOMMS_DEVICE: "CUDA"
         agents:
           slurm_mem: 20GB
           slurm_gpus: 1
 
-      - label: "GPU AMIP: ED only + bucket"
-        key: "gpu_amip_edonly_bucket"
-        command: "julia --color=yes --project=experiments/ClimaEarth/ experiments/ClimaEarth/run_amip.jl --config_file $CONFIG_PATH/amip_edonly_bucket.yml --job_id gpu_amip_edonly_bucket"
-        artifact_paths: "experiments/ClimaEarth/output/gpu_amip_edonly_bucket/artifacts/*"
-        env:
-          CLIMACOMMS_DEVICE: "CUDA"
-        agents:
-          slurm_mem: 20GB
-          slurm_gpus: 1
-
-      - label: "GPU AMIP: diag. EDMF + integrated land"
-        key: "gpu_amip_diagedmf_integrated_land"
-        command: "julia --color=yes --project=experiments/ClimaEarth/ experiments/ClimaEarth/run_amip.jl --config_file $CONFIG_PATH/amip_diagedmf_integrated_land.yml --job_id gpu_amip_diagedmf_integrated_land"
-        artifact_paths: "experiments/ClimaEarth/output/gpu_amip_diagedmf_integrated_land/artifacts/*"
-        env:
-          CLIMACOMMS_DEVICE: "CUDA"
+      - label: "CPU CMIP: ClimaAtmos + bucket land + Oceananigans + ClimaSeaIce"
+        key: "cmip_oceananigans_climaseaice"
+        command: "julia --color=yes --project=experiments/ClimaEarth/ experiments/ClimaEarth/run_amip.jl --config_file $CONFIG_PATH/cmip_oceananigans_climaseaice.yml --job_id cmip_oceananigans_climaseaice"
+        artifact_paths: "experiments/ClimaEarth/output/cmip_oceananigans_climaseaice/artifacts/*"
         agents:
           slurm_mem: 20GB
-          slurm_gpus: 1
-
-      - label: "GPU AMIP: diag. EDMF + bucket"
-        key: "gpu_amip_diagedmf_bucket"
-        command: "julia --color=yes --project=experiments/ClimaEarth/ experiments/ClimaEarth/run_amip.jl --config_file $CONFIG_PATH/amip_diagedmf_bucket.yml --job_id gpu_amip_diagedmf_bucket"
-        artifact_paths: "experiments/ClimaEarth/output/gpu_amip_diagedmf_bucket/artifacts/*"
-        env:
-          CLIMACOMMS_DEVICE: "CUDA"
-        agents:
-          slurm_mem: 20GB
-          slurm_gpus: 1
 
-      - label: "GPU AMIP test: albedo from function"
-        key: "gpu_amip_albedo_function"
-        command: "julia --color=yes --project=experiments/ClimaEarth/ experiments/ClimaEarth/run_amip.jl --config_file $CONFIG_PATH/amip_albedo_function.yml --job_id gpu_amip_albedo_function"
-        artifact_paths: "experiments/ClimaEarth/output/gpu_amip_albedo_function/artifacts/*"
+      - label: "GPU CMIP: ClimaAtmos + bucket land + Oceananigans + ClimaSeaIce"
+        key: "gpu_cmip_oceananigans_climaseaice"
+        command: "julia --color=yes --project=experiments/ClimaEarth/ experiments/ClimaEarth/run_amip.jl --config_file $CONFIG_PATH/cmip_oceananigans_climaseaice.yml --job_id gpu_cmip_oceananigans_climaseaice"
+        artifact_paths: "experiments/ClimaEarth/output/gpu_cmip_oceananigans_climaseaice/artifacts/*"
         env:
           CLIMACOMMS_DEVICE: "CUDA"
         agents:
           slurm_mem: 20GB
           slurm_gpus: 1
 
-      - label: "GPU AMIP: albedo from temporal map + 0M"
-        key: "gpu_amip_albedo_temporal_map"
-        command: "julia --color=yes --project=experiments/ClimaEarth/ experiments/ClimaEarth/run_amip.jl --config_file $CONFIG_PATH/amip_albedo_temporal_map.yml --job_id gpu_amip_albedo_temporal_map"
-        artifact_paths: "experiments/ClimaEarth/output/gpu_amip_albedo_temporal_map/artifacts/*"
-        env:
-          CLIMACOMMS_DEVICE: "CUDA"
-        agents:
-          slurm_mem: 20GB
-          slurm_gpus: 1
-
-      - label: "GPU AMIP: albedo from temporal map + 1M"
-        key: "gpu_amip_albedo_temporal_map_1M"
-        command: "julia --color=yes --project=experiments/ClimaEarth/ experiments/ClimaEarth/run_amip.jl --config_file $CONFIG_PATH/amip_albedo_temporal_map_1M.yml --job_id gpu_amip_albedo_temporal_map_1M"
-        artifact_paths: "experiments/ClimaEarth/output/gpu_amip_albedo_temporal_map_1M/artifacts/*"
-        env:
-          CLIMACOMMS_DEVICE: "CUDA"
-        agents:
-          slurm_mem: 20GB
-          slurm_gpus: 1
-
-  - group: "CMIP"
-    steps:
-
-      - label: "GPU CMIP"
-        key: "gpu_my_first_cmip"
-        command: "julia --color=yes --project=experiments/ClimaEarth/ experiments/ClimaEarth/run_amip.jl --config_file $CONFIG_PATH/my_first_cmip.yml --job_id my_first_cmip"
-        artifact_paths: "experiments/ClimaEarth/output/my_first_cmip/artifacts/*"
-        env:
-          CLIMACOMMS_DEVICE: "CUDA"
-        agents:
-          slurm_mem: 20GB
-          slurm_gpus: 1
-
-  - group: "Calibration experiments"
-    steps:
-      - label: "Perfect model calibration test"
-        key: "amip_pm_calibration"
-        command:
-          - "julia --color=yes --project=experiments/ClimaEarth experiments/calibration/run_calibration.jl"
-        artifact_paths: "experiments/calibration/output/*"
-        env:
-          CLIMACOMMS_DEVICE: "CUDA"
-          CLIMACOMMS_CONTEXT: "SINGLETON"
-          SHORT_RUN: ""
-        agents:
-          slurm_mem: 64GB
-          slurm_ntasks: 3
-          slurm_gpus_per_task: 1
-          slurm_cpus_per_task: 4
-
-  - wait
-
-  # plot job performance history
-  - label: ":chart_with_downwards_trend: build history"
-    command:
-      - build_history staging # name of branch to plot
-    artifact_paths:
-      - "build_history.html"
+  # - group: "Calibration experiments"
+  #   steps:
+  #     - label: "Perfect model calibration test"
+  #       key: "amip_pm_calibration"
+  #       command:
+  #         - "julia --color=yes --project=experiments/ClimaEarth experiments/calibration/run_calibration.jl"
+  #       artifact_paths: "experiments/calibration/output/*"
+  #       env:
+  #         CLIMACOMMS_DEVICE: "CUDA"
+  #         CLIMACOMMS_CONTEXT: "SINGLETON"
+  #         SHORT_RUN: ""
+  #       agents:
+  #         slurm_mem: 64GB
+  #         slurm_ntasks: 3
+  #         slurm_gpus_per_task: 1
+  #         slurm_cpus_per_task: 4
+
+  # - wait
+
+  # # plot job performance history
+  # - label: ":chart_with_downwards_trend: build history"
+  #   command:
+  #     - build_history staging # name of branch to plot
+  #   artifact_paths:
+  #     - "build_history.html"
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 5b4aef7ade..94f1d45f49 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -17,6 +17,9 @@ jobs:
   test:
     name: ci ${{ matrix.version }} - ${{ matrix.os }}
     runs-on: ${{ matrix.os }}
+    env:
+      ECCO_USERNAME: ${{ secrets.ECCO_USERNAME }}
+      ECCO_WEBDAV_PASSWORD: ${{ secrets.ECCO_WEBDAV_PASSWORD }}
     strategy:
       fail-fast: false
       matrix:
diff --git a/config/ci_configs/cmip_oceananigans_climaseaice.yml b/config/ci_configs/cmip_oceananigans_climaseaice.yml
new file mode 100644
index 0000000000..f3568a0eb1
--- /dev/null
+++ b/config/ci_configs/cmip_oceananigans_climaseaice.yml
@@ -0,0 +1,24 @@
+FLOAT_TYPE: "Float32"
+albedo_model: "CouplerAlbedo"
+atmos_config_file: "config/atmos_configs/climaatmos_edonly.yml"
+bucket_albedo_type: "map_temporal"
+coupler_toml: ["toml/amip.toml"]
+dt: "120secs"
+dt_cpl: "120secs"
+dz_bottom: 100.0
+energy_check: false
+h_elem: 8
+ice_model: "clima_seaice"
+mode_name: "cmip"
+netcdf_output_at_levels: true
+output_default_diagnostics: true
+radiation_reset_rng_seed: true
+rayleigh_sponge: true
+start_date: "20100101"
+surface_setup: "PrescribedSurface"
+t_end: "1days"
+topo_smoothing: true
+topography: "Earth"
+viscous_sponge: true
+z_elem: 39
+z_max: 60000.0
diff --git a/config/ci_configs/my_first_cmip.yml b/config/ci_configs/cmip_oceananigans_prescrseaice.yml
similarity index 100%
rename from config/ci_configs/my_first_cmip.yml
rename to config/ci_configs/cmip_oceananigans_prescrseaice.yml
diff --git a/docs/src/fieldexchanger.md b/docs/src/fieldexchanger.md
index 6be11ffa0e..eb66d06c12 100644
--- a/docs/src/fieldexchanger.md
+++ b/docs/src/fieldexchanger.md
@@ -39,6 +39,6 @@ the atmosphere and each surface model.
 
 ```@docs
     ClimaCoupler.FieldExchanger.combine_surfaces!
-    ClimaCoupler.FieldExchanger.resolve_ocean_ice_fractions!
+    ClimaCoupler.FieldExchanger.resolve_area_fractions!
     ClimaCoupler.FieldExchanger.import_atmos_fields!
 ```
diff --git a/experiments/ClimaEarth/Project.toml b/experiments/ClimaEarth/Project.toml
index ff292ea7e6..8adaa52bd1 100644
--- a/experiments/ClimaEarth/Project.toml
+++ b/experiments/ClimaEarth/Project.toml
@@ -41,7 +41,7 @@ ClimaAtmos = "0.27, 0.28, 0.29, 0.30, 0.31"
 ClimaCalibrate = "0.1"
 ClimaDiagnostics = "0.2.6"
 ClimaLand = "1.0"
-ClimaOcean = "0.8"
+ClimaOcean = "0.8.6"
 ClimaParams = "1.0"
 ClimaSeaIce = "0.3"
 ClimaTimeSteppers = "0.7, 0.8"
diff --git a/experiments/ClimaEarth/cli_options.jl b/experiments/ClimaEarth/cli_options.jl
index b7d35f5b46..f2feed1fce 100644
--- a/experiments/ClimaEarth/cli_options.jl
+++ b/experiments/ClimaEarth/cli_options.jl
@@ -175,6 +175,11 @@ function argparse_settings()
         help = "Directory containing ERA5 initial condition files (subseasonal mode). Filenames inferred from start_date [none (default)]. Generated with `https://github.com/CliMA/WeatherQuest`"
         arg_type = String
         default = nothing
+        # Ice model specific
+        "--ice_model"
+        help = "Sea ice model to use. [`prescribed` (default), `clima_seaice`]"
+        arg_type = String
+        default = "prescribed"
     end
     return s
 end
diff --git a/experiments/ClimaEarth/components/ocean/clima_seaice.jl b/experiments/ClimaEarth/components/ocean/clima_seaice.jl
new file mode 100644
index 0000000000..1c2661c6b0
--- /dev/null
+++ b/experiments/ClimaEarth/components/ocean/clima_seaice.jl
@@ -0,0 +1,416 @@
+import Oceananigans as OC
+import ClimaSeaIce as CSI
+import ClimaOcean as CO
+import ClimaCoupler: Checkpointer, FieldExchanger, FluxCalculator, Interfacer, Utilities
+import ClimaComms
+import ClimaCore as CC
+import Thermodynamics as TD
+import ClimaOcean.EN4: download_dataset
+using KernelAbstractions: @kernel, @index, @inbounds
+
+include("climaocean_helpers.jl")
+
+# Rename ECCO password env variable to match ClimaOcean.jl
+haskey(ENV, "ECCO_PASSWORD") && (ENV["ECCO_WEBDAV_PASSWORD"] = ENV["ECCO_PASSWORD"])
+
+"""
+    ClimaSeaIceSimulation{SIM, A, OPROP, REMAP}
+
+The ClimaCoupler simulation object used to run with ClimaSeaIce.
+This type is used by the coupler to indicate that this simulation
+is an surface/ocean simulation for dispatch.
+
+It contains the following objects:
+- `ice::SIM`: The ClimaSeaIce simulation object.
+- `area_fraction::A`: A ClimaCore Field representing the surface area fraction of this component model on the exchange grid.
+- `melting_speed::MS`: An constant characteristic speed for melting/freezing.
+- `remapping::REMAP`: Objects needed to remap from the exchange (spectral) grid to Oceananigans spaces.
+- `ocean_ice_fluxes::NT`: A NamedTuple of fluxes between the ocean and sea ice, computed at each coupling step.
+"""
+struct ClimaSeaIceSimulation{SIM, A, MS, REMAP, NT} <: Interfacer.SeaIceModelSimulation
+    ice::SIM
+    area_fraction::A
+    melting_speed::MS
+    remapping::REMAP
+    ocean_ice_fluxes::NT
+end
+
+"""
+    ClimaSeaIceSimulation()
+
+Creates an ClimaSeaIceSimulation object containing a model, an integrator, and
+a surface area fraction field.
+This type is used to indicate that this simulation is an ocean simulation for
+dispatch in coupling.
+
+Initially, no sea ice is present.
+
+Since this model does not solve for prognostic temperature, we use a
+prescribed heat flux boundary condition at the top, which is used to solve for
+temperature at the surface of the sea ice. The surface temperature from the
+previous step is provided to the coupler to be used in computing fluxes.
+
+Specific details about the default model configuration
+can be found in the documentation for `ClimaOcean.ocean_simulation`.
+"""
+function ClimaSeaIceSimulation(ocean; output_dir, start_date = nothing)
+    # Initialize the sea ice with the same grid as the ocean
+    grid = ocean.ocean.model.grid
+    arch = OC.Architectures.architecture(grid)
+    advection = ocean.ocean.model.advection.T
+    top_heat_boundary_condition = CSI.MeltingConstrainedFluxBalance()
+
+    ice = CO.sea_ice_simulation(grid, ocean.ocean; advection, top_heat_boundary_condition)
+
+    # Initialize nonzero sea ice if start date provided
+    if !isnothing(start_date)
+        sic_metadata = CO.DataWrangling.Metadatum(
+            :sea_ice_concentration,
+            dataset = CO.DataWrangling.ECCO.ECCO4Monthly(),
+            date = start_date,
+        )
+        h_metadata = CO.DataWrangling.Metadatum(
+            :sea_ice_thickness,
+            dataset = CO.DataWrangling.ECCO.ECCO4Monthly(),
+            date = start_date,
+        )
+
+        OC.set!(ice.model.ice_concentration, sic_metadata)
+        OC.set!(ice.model.ice_thickness, h_metadata)
+    end
+
+    melting_speed = 1e-4
+
+    # Since ocean and sea ice share the same grid, we can also share the remapping objects
+    remapping = ocean.remapping
+
+    # Before version 0.96.22, the NetCDFWriter was broken on GPU
+    if arch isa OC.CPU || pkgversion(OC) >= v"0.96.22"
+        # Save all tracers and velocities to a NetCDF file at daily frequency
+        outputs = OC.prognostic_fields(ice.model)
+        jld_writer = OC.JLD2Writer(
+            ice.model,
+            outputs;
+            schedule = OC.TimeInterval(86400), # Daily output
+            filename = joinpath(output_dir, "seaice_diagnostics.jld2"),
+            overwrite_existing = true,
+            array_type = Array{Float32},
+        )
+        ice.output_writers[:diagnostics] = jld_writer
+    end
+
+    # Allocate space for the sea ice-ocean (io) fluxes
+    io_bottom_heat_flux = OC.Field{OC.Center, OC.Center, Nothing}(grid)
+    io_frazil_heat_flux = OC.Field{OC.Center, OC.Center, Nothing}(grid)
+    io_salt_flux = OC.Field{OC.Center, OC.Center, Nothing}(grid)
+    x_momentum = OC.Field{OC.Face, OC.Center, Nothing}(grid)
+    y_momentum = OC.Field{OC.Center, OC.Face, Nothing}(grid)
+
+    ocean_ice_fluxes = (
+        interface_heat = io_bottom_heat_flux,
+        frazil_heat = io_frazil_heat_flux,
+        salt = io_salt_flux,
+        x_momentum = x_momentum,
+        y_momentum = y_momentum,
+    )
+
+    # Get the initial area fraction from the fractional ice concentration
+    boundary_space = axes(ocean.area_fraction)
+    FT = CC.Spaces.undertype(boundary_space)
+    area_fraction = Interfacer.remap(ice.model.ice_concentration, boundary_space)
+
+    sim = ClimaSeaIceSimulation(
+        ice,
+        area_fraction,
+        melting_speed,
+        remapping,
+        ocean_ice_fluxes,
+    )
+    return sim
+end
+
+###############################################################################
+### Functions required by ClimaCoupler.jl for a SurfaceModelSimulation
+###############################################################################
+
+# Timestep the simulation forward to time `t`
+Interfacer.step!(sim::ClimaSeaIceSimulation, t) =
+    OC.time_step!(sim.ice, float(t) - sim.ice.model.clock.time)
+
+Interfacer.get_field(sim::ClimaSeaIceSimulation, ::Val{:area_fraction}) = sim.area_fraction
+Interfacer.get_field(sim::ClimaSeaIceSimulation, ::Val{:ice_concentration}) =
+    sim.ice.model.ice_concentration
+
+# At the moment, we return always Float32. This is because we always want to run
+# Oceananingans with Float64, so we have no way to know the float type here. Sticking with
+# Float32 ensures that nothing is accidentally promoted to Float64. We will need to change
+# this anyway.
+Interfacer.get_field(sim::ClimaSeaIceSimulation, ::Val{:roughness_buoyancy}) =
+    Float32(5.8e-5)
+Interfacer.get_field(sim::ClimaSeaIceSimulation, ::Val{:roughness_momentum}) =
+    Float32(5.8e-5)
+Interfacer.get_field(sim::ClimaSeaIceSimulation, ::Val{:beta}) = Float32(1)
+Interfacer.get_field(sim::ClimaSeaIceSimulation, ::Val{:emissivity}) = Float32(1)
+Interfacer.get_field(sim::ClimaSeaIceSimulation, ::Val{:surface_direct_albedo}) =
+    Float32(0.7)
+Interfacer.get_field(sim::ClimaSeaIceSimulation, ::Val{:surface_diffuse_albedo}) =
+    Float32(0.7)
+
+# Approximate the sea ice surface temperature as the temperature computed from the
+#  fluxes at the previous timestep.
+Interfacer.get_field(sim::ClimaSeaIceSimulation, ::Val{:surface_temperature}) =
+    273.15 + sim.ice.model.ice_thermodynamics.top_surface_temperature
+
+"""
+    FluxCalculator.update_turbulent_fluxes!(sim::ClimaSeaIceSimulation, fields)
+
+Update the turbulent fluxes in the simulation using the values stored in the coupler fields.
+These include latent heat flux, sensible heat flux, momentum fluxes, and moisture flux.
+
+The input `fields` are already area-weighted, so there's no need to weight them again.
+
+Note that currently the moisture flux has no effect on the sea ice model, which has
+constant salinity.
+
+A note on sign conventions:
+SurfaceFluxes and ClimaSeaIce both use the convention that a positive flux is an upward flux.
+No sign change is needed during the exchange, except for moisture/salinity fluxes:
+SurfaceFluxes provides moisture moving from atmosphere to ocean as a negative flux at the surface,
+and ClimaSeaIce represents moisture moving from atmosphere to ocean as a positive salinity flux,
+so a sign change is needed when we convert from moisture to salinity flux.
+"""
+function FluxCalculator.update_turbulent_fluxes!(sim::ClimaSeaIceSimulation, fields)
+    # Only LatitudeLongitudeGrid are supported because otherwise we have to rotate the vectors
+
+    (; F_lh, F_sh, F_turb_ρτxz, F_turb_ρτyz, F_turb_moisture) = fields
+    grid = sim.ice.model.grid
+    ice_concentration = sim.ice.model.ice_concentration
+
+    # Remap momentum fluxes onto reduced 2D Center, Center fields using scratch arrays and fields
+    CC.Remapping.interpolate!(
+        sim.remapping.scratch_arr1,
+        sim.remapping.remapper_cc,
+        F_turb_ρτxz,
+    )
+    OC.set!(sim.remapping.scratch_cc1, sim.remapping.scratch_arr1) # zonal momentum flux
+    CC.Remapping.interpolate!(
+        sim.remapping.scratch_arr2,
+        sim.remapping.remapper_cc,
+        F_turb_ρτyz,
+    )
+    OC.set!(sim.remapping.scratch_cc2, sim.remapping.scratch_arr2) # meridional momentum flux
+
+    # Rename for clarity; these are now Center, Center Oceananigans fields
+    F_turb_ρτxz_cc = sim.remapping.scratch_cc1
+    F_turb_ρτyz_cc = sim.remapping.scratch_cc2
+
+    # Set the momentum flux BCs at the correct locations using the remapped scratch fields
+    # Note that this requires the sea ice model to always be run with dynamics turned on
+    si_flux_u = sim.ice.model.dynamics.external_momentum_stresses.top.u
+    si_flux_v = sim.ice.model.dynamics.external_momentum_stresses.top.v
+    set_from_extrinsic_vector!(
+        (; u = si_flux_u, v = si_flux_v),
+        grid,
+        F_turb_ρτxz_cc,
+        F_turb_ρτyz_cc,
+    )
+
+    # Remap the latent and sensible heat fluxes using scratch arrays
+    CC.Remapping.interpolate!(sim.remapping.scratch_arr1, sim.remapping.remapper_cc, F_lh) # latent heat flux
+    CC.Remapping.interpolate!(sim.remapping.scratch_arr2, sim.remapping.remapper_cc, F_sh) # sensible heat flux
+
+    # Rename for clarity; recall F_turb_energy = F_lh + F_sh
+    remapped_F_lh = sim.remapping.scratch_arr1
+    remapped_F_sh = sim.remapping.scratch_arr2
+
+    # Update the sea ice only where the concentration is greater than zero.
+    si_flux_heat = sim.ice.model.external_heat_fluxes.top
+    OC.interior(si_flux_heat, :, :, 1) .+=
+        (OC.interior(ice_concentration, :, :, 1) .> 0) .* (remapped_F_lh .+ remapped_F_sh)
+
+    return nothing
+end
+
+function Interfacer.update_field!(sim::ClimaSeaIceSimulation, ::Val{:area_fraction}, field)
+    sim.area_fraction .= field
+    return nothing
+end
+
+"""
+    FieldExchanger.update_sim!(sim::ClimaSeaIceSimulation, csf)
+
+Update the sea ice simulation with the provided fields, which have been filled in
+by the coupler.
+
+Update the portion of the surface_fluxes for T and S that is due to radiation and
+precipitation. The rest will be updated in `update_turbulent_fluxes!`.
+
+Note that currently precipitation has no effect on the sea ice model, which has
+constant salinity.
+
+A note on sign conventions:
+ClimaAtmos and ClimaSeaIce both use the convention that a positive flux is an upward flux.
+No sign change is needed during the exchange, except for precipitation/salinity fluxes.
+ClimaAtmos provides precipitation as a negative flux at the surface, and
+ClimaSeaIce represents precipitation as a positive salinity flux,
+so a sign change is needed when we convert from precipitation to salinity flux.
+"""
+function FieldExchanger.update_sim!(sim::ClimaSeaIceSimulation, csf)
+    ice_concentration = sim.ice.model.ice_concentration
+
+    # Remap radiative flux onto scratch array; rename for clarity
+    CC.Remapping.interpolate!(
+        sim.remapping.scratch_arr1,
+        sim.remapping.remapper_cc,
+        csf.SW_d,
+    )
+    remapped_SW_d = sim.remapping.scratch_arr1
+
+    CC.Remapping.interpolate!(
+        sim.remapping.scratch_arr2,
+        sim.remapping.remapper_cc,
+        csf.LW_d,
+    )
+    remapped_LW_d = sim.remapping.scratch_arr2
+
+    # Update only the part due to radiative fluxes. For the full update, the component due
+    # to latent and sensible heat is missing and will be updated in update_turbulent_fluxes.
+    si_flux_heat = sim.ice.model.external_heat_fluxes.top
+    # TODO: get sigma from parameters
+    σ = 5.67e-8
+    α = Interfacer.get_field(sim, Val(:surface_direct_albedo)) # scalar
+    ϵ = 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) .* .-(1 .- α) .* remapped_SW_d .-
+        ϵ .* (
+            remapped_LW_d .-
+            σ .*
+            (
+                273.15 .+ OC.interior(
+                    sim.ice.model.ice_thermodynamics.top_surface_temperature,
+                    :,
+                    :,
+                    1,
+                )
+            ) .^ 4
+        )
+    return nothing
+end
+
+"""
+    ocean_seaice_fluxes!(ocean_sim::OceananigansSimulation, ice_sim::ClimaSeaIceSimulation)
+
+Compute the fluxes between the ocean and sea ice, storing them in the `ocean_ice_fluxes`
+fields of the ocean and sea ice simulations.
+
+This function assumes both simulations share the same grid, so no remapping is done.
+
+Both simulations have had their atmospheric fluxes updated already in this timestep
+(see `update_sim!` and `update_turbulent_fluxes!`), so we add the contributions from the
+ocean-sea ice interactions to the existing fluxes, rather than overwriting all fluxes.
+
+!!! note
+    This function must be called after the turbulent fluxes have been updated in both
+    simulations. Here only the contributions from the sea ice/ocean interactions
+    are added to the fluxes.
+"""
+function FluxCalculator.ocean_seaice_fluxes!(
+    ocean_sim::OceananigansSimulation,
+    ice_sim::ClimaSeaIceSimulation,
+)
+    melting_speed = ice_sim.melting_speed
+    ocean_properties = ocean_sim.ocean_properties
+    ice_concentration = Interfacer.get_field(ice_sim, Val(:ice_concentration))
+
+    # Update the sea ice concentration in the ocean simulation
+    ocean_sim.ice_concentration .= ice_concentration
+
+    # Compute the fluxes and store them in the both simulations
+    CO.OceanSeaIceModels.InterfaceComputations.compute_sea_ice_ocean_fluxes!(
+        ice_sim.ocean_ice_fluxes,
+        ocean_sim.ocean,
+        ice_sim.ice,
+        melting_speed,
+        ocean_properties,
+    )
+
+    ## Update the internals of the sea ice model
+    # Set the bottom heat flux to the sum of the frazil and interface heat fluxes
+    bottom_heat_flux = ice_sim.ice.model.external_heat_fluxes.bottom
+
+    Qf = ice_sim.ocean_ice_fluxes.frazil_heat        # frazil heat flux
+    Qi = ice_sim.ocean_ice_fluxes.interface_heat     # interfacial heat flux
+    bottom_heat_flux .= Qf .+ Qi
+
+    ## Update the internals of the ocean model
+    ρₒ⁻¹ = 1 / ocean_sim.ocean_properties.reference_density
+    cₒ = ocean_sim.ocean_properties.heat_capacity
+
+    # Compute fluxes for u, v, T, and S from momentum, heat, and freshwater fluxes
+    oc_flux_u = surface_flux(ocean_sim.ocean.model.velocities.u)
+    oc_flux_v = surface_flux(ocean_sim.ocean.model.velocities.v)
+
+    ρτxio = ice_sim.ocean_ice_fluxes.x_momentum # sea_ice - ocean zonal momentum flux
+    ρτyio = ice_sim.ocean_ice_fluxes.y_momentum # sea_ice - ocean meridional momentum flux
+
+    # Update the momentum flux contributions from ocean/sea ice fluxes
+    grid = ocean_sim.ocean.model.grid
+    arch = OC.Architectures.architecture(grid)
+    OC.Utils.launch!(
+        arch,
+        grid,
+        :xy,
+        _add_ocean_ice_stress!,
+        oc_flux_u,
+        oc_flux_v,
+        grid,
+        ρτxio,
+        ρτyio,
+        ρₒ⁻¹,
+        ice_concentration,
+    )
+
+    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_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)
+
+    return nothing
+end
+
+@kernel function _add_ocean_ice_stress!(
+    oc_flux_u,
+    oc_flux_v,
+    grid,
+    ρτxio,
+    ρτyio,
+    ρₒ⁻¹,
+    ice_concentration,
+)
+    i, j = @index(Global, NTuple)
+
+    # ℑxᶠᵃᵃ: interpolate faces to centers
+    oc_flux_u[i, j, 1] +=
+        ρτxio[i, j, 1] * ρₒ⁻¹ * OC.Operators.ℑxᶠᵃᵃ(i, j, 1, grid, ice_concentration)
+    oc_flux_v[i, j, 1] +=
+        ρτyio[i, j, 1] * ρₒ⁻¹ * OC.Operators.ℑyᵃᶠᵃ(i, j, 1, grid, ice_concentration)
+end
+
+"""
+    get_model_prog_state(sim::ClimaSeaIceSimulation)
+
+Returns the model state of a simulation as a `ClimaCore.FieldVector`.
+It's okay to leave this unimplemented for now, but we won't be able to use the
+restart system.
+
+TODO extend this for non-ClimaCore states.
+"""
+function Checkpointer.get_model_prog_state(sim::ClimaSeaIceSimulation)
+    @warn "get_model_prog_state not implemented for ClimaSeaIceSimulation"
+end
diff --git a/experiments/ClimaEarth/components/ocean/oceananigans.jl b/experiments/ClimaEarth/components/ocean/oceananigans.jl
index b665fb981b..d8cef89fe7 100644
--- a/experiments/ClimaEarth/components/ocean/oceananigans.jl
+++ b/experiments/ClimaEarth/components/ocean/oceananigans.jl
@@ -10,7 +10,7 @@ using KernelAbstractions: @kernel, @index, @inbounds
 include("climaocean_helpers.jl")
 
 """
-    OceananigansSimulation{SIM, A, OPROP, REMAP}
+    OceananigansSimulation{SIM, A, OPROP, REMAP, SIC}
 
 The ClimaCoupler simulation object used to run with Oceananigans.
 This type is used by the coupler to indicate that this simulation
@@ -21,12 +21,14 @@ It contains the following objects:
 - `area_fraction::A`: A ClimaCore Field representing the surface area fraction of this component model on the exchange grid.
 - `ocean_properties::OPROP`: A NamedTuple of ocean properties and parameters
 - `remapping::REMAP`: Objects needed to remap from the exchange (spectral) grid to Oceananigans spaces.
+- `ice_concentration::SIC`: An Oceananigans Field representing the sea ice concentration on the ocean/sea ice grid.
 """
-struct OceananigansSimulation{SIM, A, OPROP, REMAP} <: Interfacer.OceanModelSimulation
+struct OceananigansSimulation{SIM, A, OPROP, REMAP, SIC} <: Interfacer.OceanModelSimulation
     ocean::SIM
     area_fraction::A
     ocean_properties::OPROP
     remapping::REMAP
+    ice_concentration::SIC
 end
 
 """
@@ -41,7 +43,7 @@ Specific details about the default model configuration
 can be found in the documentation for `ClimaOcean.ocean_simulation`.
 """
 function OceananigansSimulation(
-    area_fraction,
+    boundary_space,
     start_date,
     stop_date;
     output_dir,
@@ -76,7 +78,7 @@ function OceananigansSimulation(
     bottom_height = CO.regrid_bathymetry(
         underlying_grid;
         minimum_depth = 30,
-        interpolation_passes = 20,
+        interpolation_passes = 1, # TODO revert
         major_basins = 1,
     )
 
@@ -86,6 +88,7 @@ function OceananigansSimulation(
         active_cells_map = true,
     )
 
+    # TODO use_restoring only if not using ClimaSeaIce
     use_restoring = start_date + Dates.Month(1) < stop_date
 
     if use_restoring
@@ -138,7 +141,7 @@ function OceananigansSimulation(
     lat_cc = reshape(lat_cc, 1, length(lat_cc))
     target_points_cc = @. CC.Geometry.LatLongPoint(lat_cc, long_cc)
     # TODO: We can remove the `nothing` after CC > 0.14.33
-    remapper_cc = CC.Remapping.Remapper(axes(area_fraction), target_points_cc, nothing)
+    remapper_cc = CC.Remapping.Remapper(boundary_space, target_points_cc, nothing)
 
     # Construct two 2D Center/Center fields to use as scratch space while remapping
     scratch_cc1 = OC.Field{OC.Center, OC.Center, Nothing}(grid)
@@ -151,14 +154,13 @@ function OceananigansSimulation(
     interpolated_values_dim..., _buffer_length = size(remapper_cc._interpolated_values)
     scratch_arr1 = ArrayType(zeros(FT, interpolated_values_dim...))
     scratch_arr2 = ArrayType(zeros(FT, interpolated_values_dim...))
+    scratch_arr3 = ArrayType(zeros(FT, interpolated_values_dim...))
 
-    remapping = (; remapper_cc, scratch_cc1, scratch_cc2, scratch_arr1, scratch_arr2)
+    remapping =
+        (; remapper_cc, scratch_cc1, scratch_cc2, scratch_arr1, scratch_arr2, scratch_arr3)
 
-    ocean_properties = (;
-        ocean_reference_density = 1020,
-        ocean_heat_capacity = 3991,
-        ocean_fresh_water_density = 999.8,
-    )
+    ocean_properties =
+        (; reference_density = 1020, heat_capacity = 3991, fresh_water_density = 999.8)
 
     # Before version 0.96.22, the NetCDFWriter was broken on GPU
     if arch isa OC.CPU || pkgversion(OC) >= v"0.96.22"
@@ -176,20 +178,38 @@ function OceananigansSimulation(
         ocean.output_writers[:diagnostics] = netcdf_writer
     end
 
-    sim = OceananigansSimulation(ocean, area_fraction, ocean_properties, remapping)
-    return sim
+    # Initialize with 0 ice concentration; this will be updated in `resolve_area_fractions!`
+    # if the ocean is coupled to a non-prescribed sea ice model.
+    ice_concentration = OC.Field{OC.Center, OC.Center, Nothing}(grid)
+
+    # Create a dummy area fraction that will get overwritten in `update_surface_fractions!`
+    area_fraction = ones(boundary_space)
+
+    return OceananigansSimulation(
+        ocean,
+        area_fraction,
+        ocean_properties,
+        remapping,
+        ice_concentration,
+    )
 end
 
 """
-    FieldExchanger.resolve_ocean_ice_fractions!(ocean_sim, ice_sim, land_fraction)
+    FieldExchanger.resolve_area_fractions!(ocean_sim, ice_sim, land_fraction)
 
 Ensure the ocean and ice area fractions are consistent with each other.
 This matters in the case of a LatitudeLongitudeGrid, which is only
 defined between -80 and 80 degrees latitude. In this case, we want to
 set the ice fraction to `1 - land_fraction` on [-90, -80] and [80, 90]
 degrees latitude, and make sure the ocean fraction is 0 there.
+
+The land fraction is expected to be set to 1 at the poles before calling this function,
+and doesn't need to be set again since its fraction is static.
+
+This function also updates the ice concentration field in the ocean simulation
+so that it can be used for weighting flux updates.
 """
-function FieldExchanger.resolve_ocean_ice_fractions!(
+function FieldExchanger.resolve_area_fractions!(
     ocean_sim::OceananigansSimulation,
     ice_sim,
     land_fraction,
@@ -205,10 +225,17 @@ function FieldExchanger.resolve_ocean_ice_fractions!(
         polar_mask = CC.Fields.zeros(boundary_space)
         polar_mask .= abs.(lat) .>= FT(80)
 
-        # Set ice fraction to 1 - land_fraction and ocean fraction to 0 where polar_mask is 1
-        @. ice_fraction = ifelse.(polar_mask == FT(1), FT(1) - land_fraction, ice_fraction)
+        # Set land fraction to 1 and ice/ocean fraction to 0 where polar_mask is 1
+        @. land_fraction = ifelse.(polar_mask == FT(1), FT(1), land_fraction)
+        @. ice_fraction = ifelse.(polar_mask == FT(1), FT(0), ice_fraction)
         @. ocean_fraction = ifelse.(polar_mask == FT(1), FT(0), ocean_fraction)
     end
+
+    # Update the ice concentration field in the ocean simulation
+    ice_sim isa ClimaSeaIceSimulation && (
+        ocean_sim.ice_concentration .=
+            Interfacer.get_field(ice_sim, Val(:ice_concentration))
+    )
     return nothing
 end
 
@@ -223,7 +250,6 @@ Interfacer.step!(sim::OceananigansSimulation, t) =
 Interfacer.get_field(sim::OceananigansSimulation, ::Val{:area_fraction}) = sim.area_fraction
 
 # TODO: Better values for this
-
 # At the moment, we return always Float32. This is because we always want to run
 # Oceananingans with Float64, so we have no way to know the float type here. Sticking with
 # Float32 ensures that nothing is accidentally promoted to Float64. We will need to change
@@ -246,9 +272,15 @@ Interfacer.get_field(sim::OceananigansSimulation, ::Val{:surface_temperature}) =
 """
     FluxCalculator.update_turbulent_fluxes!(sim::OceananigansSimulation, fields)
 
-Update the turbulent fluxes in the simulation using the values stored in the coupler fields.
+Update the turbulent fluxes in the simulation using the values computed at this time step.
 These include latent heat flux, sensible heat flux, momentum fluxes, and moisture flux.
 
+Rather than setting the surface fluxes and overwriting previous values, this function adds only
+the contributions from the turbulent fluxes. `update_sim!` sets the surface fluxes due to
+radiation and precipitation. Additional contributions may be made in `ocean_seaice_fluxes!`.
+An exception is the momentum fluxes, which are set directly here since they are not updated
+in `update_sim!`.
+
 A note on sign conventions:
 SurfaceFluxes and Oceananigans both use the convention that a positive flux is an upward flux.
 No sign change is needed during the exchange, except for moisture/salinity fluxes:
@@ -259,6 +291,7 @@ so a sign change is needed when we convert from moisture to salinity flux.
 function FluxCalculator.update_turbulent_fluxes!(sim::OceananigansSimulation, fields)
     (; F_lh, F_sh, F_turb_ρτxz, F_turb_ρτyz, F_turb_moisture) = fields
     grid = sim.ocean.model.grid
+    ice_concentration = sim.ice_concentration
 
     # Remap momentum fluxes onto reduced 2D Center, Center fields using scratch arrays and fields
     CC.Remapping.interpolate!(
@@ -278,6 +311,13 @@ function FluxCalculator.update_turbulent_fluxes!(sim::OceananigansSimulation, fi
     F_turb_ρτxz_cc = sim.remapping.scratch_cc1
     F_turb_ρτyz_cc = sim.remapping.scratch_cc2
 
+    # Weight by (1 - sea ice concentration)
+    # TODO does this work with OC fields?
+    OC.interior(F_turb_ρτxz_cc, :, :, 1) .=
+        OC.interior(F_turb_ρτxz_cc, :, :, 1) .* (1.0 .- ice_concentration)
+    OC.interior(F_turb_ρτyz_cc, :, :, 1) .=
+        OC.interior(F_turb_ρτyz_cc, :, :, 1) .* (1.0 .- ice_concentration)
+
     # Set the momentum flux BCs at the correct locations using the remapped scratch fields
     oc_flux_u = surface_flux(sim.ocean.model.velocities.u)
     oc_flux_v = surface_flux(sim.ocean.model.velocities.v)
@@ -288,8 +328,7 @@ function FluxCalculator.update_turbulent_fluxes!(sim::OceananigansSimulation, fi
         F_turb_ρτyz_cc,
     )
 
-    (; ocean_reference_density, ocean_heat_capacity, ocean_fresh_water_density) =
-        sim.ocean_properties
+    (; reference_density, heat_capacity, fresh_water_density) = sim.ocean_properties
 
     # Remap the latent and sensible heat fluxes using scratch arrays
     CC.Remapping.interpolate!(sim.remapping.scratch_arr1, sim.remapping.remapper_cc, F_lh) # latent heat flux
@@ -305,7 +344,8 @@ function FluxCalculator.update_turbulent_fluxes!(sim::OceananigansSimulation, fi
     oc_flux_T = surface_flux(sim.ocean.model.tracers.T)
     OC.interior(oc_flux_T, :, :, 1) .=
         OC.interior(oc_flux_T, :, :, 1) .+
-        (remapped_F_lh .+ remapped_F_sh) ./ (ocean_reference_density * ocean_heat_capacity)
+        (1.0 .- ice_concentration) .* (remapped_F_lh .+ remapped_F_sh) ./
+        (reference_density * heat_capacity)
 
     # Add the part of the salinity flux that comes from the moisture flux, we also need to
     # add the component due to precipitation (that was done with the radiative fluxes)
@@ -314,11 +354,12 @@ function FluxCalculator.update_turbulent_fluxes!(sim::OceananigansSimulation, fi
         sim.remapping.remapper_cc,
         F_turb_moisture,
     )
-    moisture_fresh_water_flux = sim.remapping.scratch_arr1 ./ ocean_fresh_water_density
+    moisture_fresh_water_flux = sim.remapping.scratch_arr1 ./ fresh_water_density
     oc_flux_S = surface_flux(sim.ocean.model.tracers.S)
     surface_salinity = OC.interior(sim.ocean.model.tracers.S, :, :, 1)
     OC.interior(oc_flux_S, :, :, 1) .=
-        OC.interior(oc_flux_S, :, :, 1) .- surface_salinity .* moisture_fresh_water_flux
+        OC.interior(oc_flux_S, :, :, 1) .-
+        (1.0 .- ice_concentration) .* surface_salinity .* moisture_fresh_water_flux
     return nothing
 end
 
@@ -336,6 +377,9 @@ by the coupler.
 Update the portion of the surface_fluxes for T and S that is due to radiation and
 precipitation. The rest will be updated in `update_turbulent_fluxes!`.
 
+This function sets the surface fluxes directly, overwriting any previous values.
+Additional contributions will be made in `update_turbulent_fluxes!` and `ocean_seaice_fluxes!`.
+
 A note on sign conventions:
 ClimaAtmos and Oceananigans both use the convention that a positive flux is an upward flux.
 No sign change is needed during the exchange, except for precipitation/salinity fluxes.
@@ -344,8 +388,8 @@ Oceananigans represents precipitation as a positive salinity flux,
 so a sign change is needed when we convert from precipitation to salinity flux.
 """
 function FieldExchanger.update_sim!(sim::OceananigansSimulation, csf)
-    (; ocean_reference_density, ocean_heat_capacity, ocean_fresh_water_density) =
-        sim.ocean_properties
+    (; reference_density, heat_capacity, fresh_water_density) = sim.ocean_properties
+    ice_concentration = sim.ice_concentration
 
     # Remap radiative flux onto scratch array; rename for clarity
     CC.Remapping.interpolate!(
@@ -370,13 +414,13 @@ function FieldExchanger.update_sim!(sim::OceananigansSimulation, csf)
     α = Interfacer.get_field(sim, Val(:surface_direct_albedo)) # scalar
     ϵ = Interfacer.get_field(sim, Val(:emissivity)) # scalar
     OC.interior(oc_flux_T, :, :, 1) .=
-        (
+        (1.0 .- ice_concentration) .* (
             -(1 - α) .* remapped_SW_d .-
             ϵ * (
                 remapped_LW_d .-
                 σ .* (273.15 .+ OC.interior(sim.ocean.model.tracers.T, :, :, 1)) .^ 4
             )
-        ) ./ (ocean_reference_density * ocean_heat_capacity)
+        ) ./ (reference_density * heat_capacity)
 
     # Remap precipitation fields onto scratch arrays; rename for clarity
     CC.Remapping.interpolate!(
@@ -394,11 +438,10 @@ function FieldExchanger.update_sim!(sim::OceananigansSimulation, csf)
 
     # Virtual salt flux
     oc_flux_S = surface_flux(sim.ocean.model.tracers.S)
-    precipitating_fresh_water_flux =
-        (remapped_P_liq .+ remapped_P_snow) ./ ocean_fresh_water_density
-    surface_salinity_flux =
-        OC.interior(sim.ocean.model.tracers.S, :, :, 1) .* precipitating_fresh_water_flux
-    OC.interior(oc_flux_S, :, :, 1) .= .-surface_salinity_flux
+    OC.interior(oc_flux_S, :, :, 1) .=
+        OC.interior(oc_flux_S, :, :, 1) .-
+        OC.interior(sim.ocean.model.tracers.S, :, :, 1) .* (1.0 .- ice_concentration) .*
+        (remapped_P_liq .+ remapped_P_snow) ./ fresh_water_density
     return nothing
 end
 
diff --git a/experiments/ClimaEarth/components/ocean/prescr_ocean.jl b/experiments/ClimaEarth/components/ocean/prescr_ocean.jl
index f207fbfada..755fbdb525 100644
--- a/experiments/ClimaEarth/components/ocean/prescr_ocean.jl
+++ b/experiments/ClimaEarth/components/ocean/prescr_ocean.jl
@@ -65,7 +65,6 @@ function PrescribedOceanSimulation(
     space,
     start_date,
     t_start,
-    area_fraction,
     thermo_params,
     comms_ctx;
     z0m = FT(5.8e-5),
@@ -113,7 +112,7 @@ function PrescribedOceanSimulation(
         α_diffuse = ones(space) .* α_diffuse_val,
         u_int = zeros(space),
         v_int = zeros(space),
-        area_fraction = area_fraction,
+        area_fraction = ones(space),
         phase = TD.Liquid(),
         thermo_params = thermo_params,
         SST_timevaryinginput = SST_timevaryinginput,
diff --git a/experiments/ClimaEarth/components/ocean/prescr_seaice.jl b/experiments/ClimaEarth/components/ocean/prescr_seaice.jl
index ca90e9f851..af180f3bf9 100644
--- a/experiments/ClimaEarth/components/ocean/prescr_seaice.jl
+++ b/experiments/ClimaEarth/components/ocean/prescr_seaice.jl
@@ -119,12 +119,8 @@ function PrescribedIceSimulation(
     )
 
     # Get initial SIC values and use them to calculate ice fraction
-    SIC_init = CC.Fields.zeros(space)
-    evaluate!(SIC_init, SIC_timevaryinginput, tspan[1])
-
-    # Overwrite ice fraction with the static land area fraction anywhere we have nonzero land area
-    #  max needed to avoid Float32 errors (see issue #271; Heisenbug on HPC)
-    ice_fraction = @. max(min(SIC_init, FT(1) - land_fraction), FT(0))
+    ice_fraction = CC.Fields.zeros(space)
+    evaluate!(ice_fraction, SIC_timevaryinginput, tspan[1])
 
     params = IceSlabParameters{FT}()
 
diff --git a/experiments/ClimaEarth/components/ocean/slab_ocean.jl b/experiments/ClimaEarth/components/ocean/slab_ocean.jl
index 9db88a0899..ad14fa11ff 100644
--- a/experiments/ClimaEarth/components/ocean/slab_ocean.jl
+++ b/experiments/ClimaEarth/components/ocean/slab_ocean.jl
@@ -55,7 +55,7 @@ function slab_ocean_space_init(space, params)
 end
 
 """
-    SlabOceanSimulation(::Type{FT}; tspan, dt, saveat, space, area_fraction, stepper = CTS.RK4()) where {FT}
+    SlabOceanSimulation(::Type{FT}; tspan, dt, saveat, space, stepper = CTS.RK4()) where {FT}
 
 Initializes the `DiffEq` problem, and creates a Simulation-type object containing the necessary information for `step!` in the coupling loop.
 """
@@ -65,7 +65,6 @@ function SlabOceanSimulation(
     dt,
     saveat,
     space,
-    area_fraction,
     thermo_params,
     stepper = CTS.RK4(),
     evolving = true,
@@ -80,7 +79,7 @@ function SlabOceanSimulation(
         F_turb_energy = CC.Fields.zeros(space),
         SW_d = CC.Fields.zeros(space),
         LW_d = CC.Fields.zeros(space),
-        area_fraction = area_fraction,
+        area_fraction = CC.Fields.ones(space),
         thermo_params = thermo_params,
         α_direct = CC.Fields.ones(space) .* params.α,
         α_diffuse = CC.Fields.ones(space) .* params.α,
diff --git a/experiments/ClimaEarth/setup_run.jl b/experiments/ClimaEarth/setup_run.jl
index abc82c8aa6..a60f39acc4 100644
--- a/experiments/ClimaEarth/setup_run.jl
+++ b/experiments/ClimaEarth/setup_run.jl
@@ -76,6 +76,7 @@ include("components/ocean/slab_ocean.jl")
 include("components/ocean/prescr_ocean.jl")
 include("components/ocean/prescr_seaice.jl")
 include("components/ocean/oceananigans.jl")
+include("components/ocean/clima_seaice.jl")
 
 #=
 ### Configuration Dictionaries
@@ -140,6 +141,7 @@ function CoupledSimulation(config_dict::AbstractDict)
         output_dir_root,
         parameter_files,
         era5_initial_condition_dir,
+        ice_model,
     ) = get_coupler_args(config_dict)
 
     #=
@@ -343,28 +345,11 @@ function CoupledSimulation(config_dict::AbstractDict)
             error("Invalid land model specified: $(land_model)")
         end
 
-        ## sea ice model
-        ice_sim = PrescribedIceSimulation(
-            FT;
-            tspan = tspan,
-            dt = component_dt_dict["dt_seaice"],
-            saveat = saveat,
-            space = boundary_space,
-            thermo_params = thermo_params,
-            comms_ctx,
-            start_date,
-            land_fraction,
-            sic_path = subseasonal_sic,
-        )
-
-        ## ocean model using prescribed data
-        ice_fraction = Interfacer.get_field(ice_sim, Val(:area_fraction))
-        ocean_fraction = FT(1) .- ice_fraction .- land_fraction
-
+        ## ocean model
         if sim_mode <: CMIPMode
             stop_date = date(tspan[end] - tspan[begin])
             ocean_sim = OceananigansSimulation(
-                ocean_fraction,
+                boundary_space,
                 start_date,
                 stop_date;
                 output_dir = dir_paths.ocean_output_dir,
@@ -376,12 +361,34 @@ function CoupledSimulation(config_dict::AbstractDict)
                 boundary_space,
                 start_date,
                 t_start,
-                ocean_fraction,
                 thermo_params,
                 comms_ctx;
                 sst_path = subseasonal_sst,
             )
         end
+        ## sea ice model
+        if ice_model == "clima_seaice"
+            ice_sim = ClimaSeaIceSimulation(
+                ocean_sim;
+                output_dir = dir_paths.ice_output_dir,
+                start_date,
+            )
+        elseif ice_model == "prescribed"
+            ice_sim = PrescribedIceSimulation(
+                FT;
+                tspan = tspan,
+                dt = component_dt_dict["dt_seaice"],
+                saveat = saveat,
+                space = boundary_space,
+                thermo_params = thermo_params,
+                comms_ctx,
+                start_date,
+                land_fraction,
+                sic_path = subseasonal_sic,
+            )
+        else
+            error("Invalid ice model specified: $(ice_model)")
+        end
 
     elseif (sim_mode <: AbstractSlabplanetSimulationMode)
 
@@ -413,7 +420,6 @@ function CoupledSimulation(config_dict::AbstractDict)
             dt = component_dt_dict["dt_ocean"],
             space = boundary_space,
             saveat = saveat,
-            area_fraction = (FT(1) .- land_fraction), ## NB: this ocean fraction includes areas covered by sea ice (unlike the one contained in the cs)
             thermo_params = thermo_params,
             evolving = evolving_ocean,
         )
@@ -534,8 +540,8 @@ function CoupledSimulation(config_dict::AbstractDict)
 
         The concrete steps for proper initialization are:
         =#
-
         # 1. Make sure surface model area fractions sum to 1 everywhere.
+        #  Note that ocean and ice fractions are not accurate until after this call.
         FieldExchanger.update_surface_fractions!(cs)
 
         # 2. Import atmospheric and surface fields into the coupler fields,
@@ -545,9 +551,14 @@ function CoupledSimulation(config_dict::AbstractDict)
         # 3. Update any fields in the model caches that can only be filled after the initial exchange.
         FieldExchanger.set_caches!(cs)
 
+        # TODO think about calling exchange again here to provide correct rad fluxes to surfaces
+
         # 4. Calculate and update turbulent fluxes for each surface model,
         #  and save the weighted average in coupler fields
         FluxCalculator.turbulent_fluxes!(cs)
+
+        # 5. Compute any ocean-sea ice fluxes
+        FluxCalculator.ocean_seaice_fluxes!(cs)
     end
     Utilities.show_memory_usage()
     return cs
@@ -682,12 +693,15 @@ function step!(cs::CoupledSimulation)
     ## update the surface fractions for surface models
     FieldExchanger.update_surface_fractions!(cs)
 
-    ## exchange all non-turbulent flux fields between models
+    ## exchange all non-turbulent flux fields between models, including radiative and precipitation fluxes
     FieldExchanger.exchange!(cs)
 
     ## calculate turbulent fluxes in the coupler and update the model simulations with them
     FluxCalculator.turbulent_fluxes!(cs)
 
+    ## compute any ocean-sea ice fluxes
+    FluxCalculator.ocean_seaice_fluxes!(cs)
+
     ## Maybe call the callbacks
     TimeManager.callbacks!(cs)
 
diff --git a/experiments/ClimaEarth/user_io/arg_parsing.jl b/experiments/ClimaEarth/user_io/arg_parsing.jl
index 0a404b93f7..92b5e51eee 100644
--- a/experiments/ClimaEarth/user_io/arg_parsing.jl
+++ b/experiments/ClimaEarth/user_io/arg_parsing.jl
@@ -170,8 +170,12 @@ function get_coupler_args(config_dict::Dict)
     bucket_albedo_type = config_dict["bucket_albedo_type"]
     bucket_initial_condition = config_dict["bucket_initial_condition"]
 
+    # Initial condition setting
     era5_initial_condition_dir = config_dict["era5_initial_condition_dir"]
 
+    # Ice model-specific information
+    ice_model = config_dict["ice_model"]
+
     return (;
         job_id,
         sim_mode,
@@ -203,6 +207,7 @@ function get_coupler_args(config_dict::Dict)
         bucket_initial_condition,
         parameter_files,
         era5_initial_condition_dir,
+        ice_model,
     )
 end
 
diff --git a/src/FieldExchanger.jl b/src/FieldExchanger.jl
index 3e16043030..98b723db8a 100644
--- a/src/FieldExchanger.jl
+++ b/src/FieldExchanger.jl
@@ -16,7 +16,7 @@ export update_sim!,
     exchange!,
     set_caches!,
     update_surface_fractions!,
-    resolve_ocean_ice_fractions!
+    resolve_area_fractions!
 
 """
     update_surface_fractions!(cs::Interfacer.CoupledSimulation)
@@ -69,9 +69,9 @@ function update_surface_fractions!(cs::Interfacer.CoupledSimulation)
         )
         ocean_fraction = Interfacer.get_field(ocean_sim, Val(:area_fraction))
 
-        # ensure that ocean and ice fractions are consistent
+        # Apply any additional constraints on the ocean and ice fractions if necessary
         if haskey(cs.model_sims, :ice_sim)
-            resolve_ocean_ice_fractions!(ocean_sim, cs.model_sims.ice_sim, land_fraction)
+            FieldExchanger.resolve_area_fractions!(ocean_sim, cs.model_sims.ice_sim, land_fraction)
         end
     else
         cs.fields.scalar_temp1 .= 0
@@ -84,7 +84,7 @@ function update_surface_fractions!(cs::Interfacer.CoupledSimulation)
 end
 
 """
-    resolve_ocean_ice_fractions!(ocean_sim, ice_sim, land_fraction)
+    resolve_area_fractions!(ocean_sim, ice_sim, land_fraction)
 
 Ensure that the ocean and ice fractions are consistent with each other.
 For most ocean and ice models, this does nothing since the ocean fraction is
@@ -92,7 +92,7 @@ defined as `1 - ice_fraction - land_fraction`. However, some models may have
 additional constraints on the ice and ocean fractions that need to be enforced.
 This function can be extended for such models.
 """
-function resolve_ocean_ice_fractions!(ocean_sim, ice_sim, land_fraction)
+function resolve_area_fractions!(ocean_sim, ice_sim, land_fraction)
     return nothing
 end
 
@@ -209,6 +209,11 @@ end
 Updates the surface component model cache with the current coupler fields
 *besides turbulent fluxes*, which are updated in `update_turbulent_fluxes`.
 
+Note that upwelling longwave and shortwave radiation are not computed here,
+and are expected to be computed internally by the surface model.
+Some component models extend this function and compute the upwelling longwave
+and shortwave radiation in their methods of `update_sim!`.
+
 # Arguments
 - `sim`: [Interfacer.SurfaceModelSimulation] containing a surface model simulation object.
 - `csf`: [NamedTuple] containing coupler fields.
diff --git a/src/FluxCalculator.jl b/src/FluxCalculator.jl
index 866dbcdc6a..c3595eb512 100644
--- a/src/FluxCalculator.jl
+++ b/src/FluxCalculator.jl
@@ -14,7 +14,11 @@ import ClimaCore as CC
 import ..Interfacer, ..Utilities
 
 export extrapolate_ρ_to_sfc,
-    turbulent_fluxes!, get_surface_params, update_turbulent_fluxes!, compute_surface_fluxes!
+    turbulent_fluxes!,
+    get_surface_params,
+    update_turbulent_fluxes!,
+    compute_surface_fluxes!,
+    ocean_seaice_fluxes!
 
 function turbulent_fluxes!(cs::Interfacer.CoupledSimulation)
     return turbulent_fluxes!(cs.fields, cs.model_sims, cs.thermo_params)
@@ -380,4 +384,25 @@ function compute_surface_fluxes!(
     return nothing
 end
 
+"""
+    ocean_seaice_fluxes!(cs::CoupledSimulation)
+    ocean_seaice_fluxes!(ocean_sim, ice_sim)
+
+Compute the fluxes between the ocean and sea ice simulations.
+This function does nothing by default - it should be extended
+for any ocean and sea ice models that support flux calculations.
+"""
+function ocean_seaice_fluxes!(cs::Interfacer.CoupledSimulation)
+    haskey(cs.model_sims, :ocean_sim) &&
+        haskey(cs.model_sims, :ice_sim) &&
+        ocean_seaice_fluxes!(cs.model_sims.ocean_sim, cs.model_sims.ice_sim)
+    return nothing
+end
+function ocean_seaice_fluxes!(
+    ocean_sim::Union{Interfacer.OceanModelSimulation, Interfacer.AbstractSurfaceStub},
+    ice_sim::Union{Interfacer.SeaIceModelSimulation, Interfacer.AbstractSurfaceStub},
+)
+    return nothing
+end
+
 end # module
diff --git a/src/Utilities.jl b/src/Utilities.jl
index b2f8494f09..ae4739fc13 100644
--- a/src/Utilities.jl
+++ b/src/Utilities.jl
@@ -150,12 +150,14 @@ function setup_output_dirs(;
     atmos_output_dir = joinpath(output_dir_root, "clima_atmos")
     land_output_dir = joinpath(output_dir_root, "clima_land")
     ocean_output_dir = joinpath(output_dir_root, "clima_ocean")
+    ice_output_dir = joinpath(output_dir_root, "clima_seaice")
     coupler_output_dir = joinpath(output_dir_root, "clima_coupler")
 
     if ClimaComms.iamroot(comms_ctx)
         mkpath(atmos_output_dir)
         mkpath(land_output_dir)
         mkpath(ocean_output_dir)
+        mkpath(ice_output_dir)
         mkpath(coupler_output_dir)
 
         mkpath(artifacts_dir)
@@ -175,6 +177,7 @@ function setup_output_dirs(;
         atmos_output_dir,
         land_output_dir,
         ocean_output_dir,
+        ice_output_dir,
         coupler_output_dir,
         artifacts_dir,
         regrid_dir,