fixes; runs with correct (?) fluxes

Author: juliasloan25

experiments/ClimaEarth/components/ocean/clima_seaice.jl

@@ -50,7 +50,7 @@ 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(land_fraction, ocean; output_dir)
+function ClimaSeaIceSimulation(land_fraction, 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)
@@ -59,6 +59,23 @@ function ClimaSeaIceSimulation(land_fraction, ocean; output_dir)
 
     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
@@ -168,6 +185,7 @@ function FluxCalculator.update_turbulent_fluxes!(sim::ClimaSeaIceSimulation, fie
 
     (; 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!(
@@ -208,8 +226,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
-    OC.interior(si_flux_heat, :, :, 1) .=
-        OC.interior(si_flux_heat, :, :, 1) .+ (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
@@ -220,7 +238,7 @@ function Interfacer.update_field!(sim::ClimaSeaIceSimulation, ::Val{:area_fracti
 end
 
 """
-    FieldExchanger.update_sim!(sim::ClimaSeaIceSimulation, csf, area_fraction)
+    FieldExchanger.update_sim!(sim::ClimaSeaIceSimulation, csf)
 
 Update the sea ice simulation with the provided fields, which have been filled in
 by the coupler.
@@ -238,20 +256,47 @@ 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, area_fraction)
+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.F_radiative,
+        csf.SW_d,
     )
-    remapped_F_radiative = sim.remapping.scratch_arr1
+    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
-    OC.interior(si_flux_heat, :, :, 1) .= remapped_F_radiative
+    # 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
 
@@ -352,9 +397,9 @@ end
     i, j = @index(Global, NTuple)
 
     # ℑxᶠᵃᵃ: interpolate faces to centers
-    oc_flux_u +=
+    oc_flux_u[i, j, 1] +=
         ρτxio[i, j, 1] * ρₒ⁻¹ * OC.Operators.ℑxᶠᵃᵃ(i, j, 1, grid, ice_concentration)
-    oc_flux_v +=
+    oc_flux_v[i, j, 1] +=
         ρτyio[i, j, 1] * ρₒ⁻¹ * OC.Operators.ℑyᵃᶠᵃ(i, j, 1, grid, ice_concentration)
 end
 

experiments/ClimaEarth/setup_run.jl

@@ -308,7 +308,7 @@ function CoupledSimulation(config_dict::AbstractDict)
         shared_surface_space = share_surface_space ? boundary_space : nothing
         if land_model == "bucket"
 
-            # TODO only for cmip
+            # TODO move this into resolve_area_fractions and call after init, think about restarts
             polar_mask = CC.Fields.zeros(boundary_space)
             lat = CC.Fields.coordinate_field(boundary_space).lat
             polar_mask .= abs.(lat) .>= FT(80)
@@ -375,6 +375,7 @@ function CoupledSimulation(config_dict::AbstractDict)
             @assert sim_mode <: CMIPMode
 
             # TODO how should we initialize ocean fraction when using ClimaSeaIce?
+            # TODO init everything with AF=1, then resolve after getting sea ice
             sic_data = try
                 joinpath(
                     @clima_artifact("historical_sst_sic", comms_ctx),
@@ -423,6 +424,14 @@ function CoupledSimulation(config_dict::AbstractDict)
                     ice_fraction,
                     ocean_sim;
                     output_dir = dir_paths.ice_output_dir,
+                    start_date,
+                )
+                # TODO don't need to initialize ocean fraction correctly if we do this
+                # TODO can rename to `resolve_area_fractions!` and also make land_fraction cover poles here instead of in driver
+                FieldExchanger.resolve_ocean_ice_fractions!(
+                    ocean_sim,
+                    ice_sim,
+                    land_fraction,
                 )
             end
         else
@@ -600,6 +609,8 @@ 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)

src/FieldExchanger.jl

@@ -217,6 +217,7 @@ function update_sim!(sim::Interfacer.SurfaceModelSimulation, csf)
     # radiative fluxes
     Interfacer.update_field!(sim, Val(:SW_d), csf.SW_d)
     Interfacer.update_field!(sim, Val(:LW_d), csf.LW_d)
+    # TODO need to compute SWU, LWU here too
 
     # precipitation
     Interfacer.update_field!(sim, Val(:liquid_precipitation), csf.P_liq)