Merge pull request #1519 from CliMA/js/zs/sf-af

fixes for sea ice and fluxes

Author: juliasloan25

experiments/ClimaEarth/components/ocean/prescr_seaice.jl

@@ -124,6 +124,7 @@ function PrescribedIceSimulation(
     # 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 = ifelse.(ice_fraction .> FT(0.5), FT(1), FT(0))
 
     params = IceSlabParameters{FT}()
 
@@ -263,17 +264,20 @@ function ice_rhs!(dY, Y, p, t)
 
     # Update the cached area fraction with the current SIC
     evaluate!(p.area_fraction, p.SIC_timevaryinginput, t)
+    @. p.area_fraction = ifelse(p.area_fraction > FT(0.5), FT(1), FT(0))
 
     # 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)
     @. p.area_fraction = max(min(p.area_fraction, FT(1) - p.land_fraction), FT(0))
 
     # Calculate the conductive flux, and set it to zero if the area fraction is zero
     F_conductive = @. params.k_ice / (params.h) * (params.T_base - Y.T_sfc) # fluxes are defined to be positive when upward
-    @. F_conductive = ifelse(p.area_fraction ≈ 0, zero(F_conductive), F_conductive)
-
     rhs = @. (-p.F_turb_energy - p.F_radiative + F_conductive) /
        (params.h * params.ρ * params.c)
+
+    # Zero out tendencies where there is no ice, so that ice temperature remains constant there
+    @. rhs = ifelse(p.area_fraction ≈ 0, zero(rhs), rhs)
+
     # If tendencies lead to temperature above freezing, set temperature to freezing
     @. dY.T_sfc = min(rhs, (params.T_freeze - Y.T_sfc) / float(p.dt))
 end

experiments/ClimaEarth/components/ocean/slab_ocean.jl

@@ -199,6 +199,9 @@ function slab_ocean_rhs!(dY, Y, cache, t)
     p, F_turb_energy, F_radiative = cache
     rhs = @. -(F_turb_energy + F_radiative) / (p.h * p.ρ * p.c)
 
+    # Zero out tendencies where there is no ocean, so that temperature remains constant there
+    @. rhs = ifelse(cache.area_fraction ≈ 0, zero(rhs), rhs)
+
     # Note that the area fraction has already been applied to the fluxes,
     #  so we don't need to multiply by it here.
     @. dY.T_sfc = rhs * p.evolving_switch

experiments/ClimaEarth/setup_run.jl

@@ -373,6 +373,7 @@ function CoupledSimulation(config_dict::AbstractDict)
 
         ## ocean model using prescribed data
         ice_fraction = Interfacer.get_field(ice_sim, Val(:area_fraction))
+        ice_fraction = ifelse.(ice_fraction .> FT(0.5), FT(1), FT(0))
         ocean_fraction = FT(1) .- ice_fraction .- land_fraction
 
         if sim_mode <: CMIPMode

experiments/ClimaEarth/test/component_model_tests/prescr_seaice_tests.jl

@@ -71,12 +71,12 @@ for FT in (Float32, Float64)
         # check that nothing changes with no fluxes (zero conduction when T_base == initial T_sfc)
         T_base_eq = IceSlabParameters{FT}().T_freeze - FT(5.0)
         dY, Y, p = test_sea_ice_rhs(T_base = T_base_eq)
-        @test all([i for i in extrema(dY)] .≈ [FT(0.0), FT(0.0)])
+        @test all(T -> T == FT(0), Array(parent(dY.T_sfc)))
 
         # check expected dT due to radiative flux only (again set T_base == initial T_sfc)
         dY, Y, p = test_sea_ice_rhs(F_radiative = 1.0, T_base = T_base_eq)
         dT_expected = -1.0 / (p.params.h * p.params.ρ * p.params.c)
-        @test all(extrema(dY) .≈ FT(dT_expected))
+        @test all(T -> T == FT(0) || T ≈ dT_expected, Array(parent(dY.T_sfc)))
 
         # check that tendency will not result in above freezing T
         dY, Y, p = test_sea_ice_rhs(F_radiative = 0.0, T_base = 330.0) # Float32 requires a large number here!

src/FluxCalculator.jl

@@ -356,26 +356,18 @@ function compute_surface_fluxes!(
     @. ustar = ifelse(area_fraction ≈ 0, zero(ustar), ustar)
     @. buoyancy_flux = ifelse(area_fraction ≈ 0, zero(buoyancy_flux), buoyancy_flux)
 
-    # multiply fluxes by area fraction
-    F_turb_ρτxz .*= area_fraction
-    F_turb_ρτyz .*= area_fraction
-    F_sh .*= area_fraction
-    F_lh .*= area_fraction
-    F_turb_moisture .*= area_fraction
-
     # update the fluxes, which are now area-weighted, of this surface model
     fields = (; F_turb_ρτxz, F_turb_ρτyz, F_lh, F_sh, F_turb_moisture)
     FluxCalculator.update_turbulent_fluxes!(sim, fields)
 
     # update fluxes in the coupler fields
     # add the flux contributing from this surface to the coupler field
-    # note that the fluxes area-weighted, so if a surface model is
-    #  not present at a point, the fluxes are zero
-    @. csf.F_turb_ρτxz += F_turb_ρτxz
-    @. csf.F_turb_ρτyz += F_turb_ρτyz
-    @. csf.F_lh += F_lh
-    @. csf.F_sh += F_sh
-    @. csf.F_turb_moisture += F_turb_moisture
+    # note that the fluxes are area-weighted here when we provide them to the atmosphere
+    @. csf.F_turb_ρτxz += F_turb_ρτxz * area_fraction
+    @. csf.F_turb_ρτyz += F_turb_ρτyz * area_fraction
+    @. csf.F_lh += F_lh * area_fraction
+    @. csf.F_sh += F_sh * area_fraction
+    @. csf.F_turb_moisture += F_turb_moisture * area_fraction
 
     # NOTE: This is still an area weighted contribution, which maybe doesn't make
     # too much sense for these quantities...