use zero(FT) instead of zero(field) in broadcasts

experiments/ClimaEarth/components/land/climaland_bucket.jl

@@ -152,9 +152,9 @@ function BucketSimulation(
         Y.bucket.T .= orog_adjusted_T
     end
 
-    Y.bucket.W .= 0.15
-    Y.bucket.Ws .= 0.0
-    Y.bucket.σS .= 0.0
+    Y.bucket.W .= FT(0.15)
+    Y.bucket.Ws .= zero(FT)
+    Y.bucket.σS .= zero(FT)
 
     # Overwrite initial conditions with interpolated values from a netcdf file using
     # the `SpaceVaryingInputs` tool. We expect the file to contain the following variables:

experiments/ClimaEarth/components/land/climaland_integrated.jl

@@ -607,10 +607,8 @@ function FluxCalculator.compute_surface_fluxes!(
     )
 
     # Zero out the fluxes where the area fraction is zero
-    @. csf.scalar_temp1 =
-        ifelse(area_fraction == 0, zero(csf.scalar_temp1), csf.scalar_temp1)
-    @. csf.scalar_temp2 =
-        ifelse(area_fraction == 0, zero(csf.scalar_temp2), csf.scalar_temp2)
+    @. csf.scalar_temp1 = ifelse(area_fraction == 0, zero(FT), csf.scalar_temp1)
+    @. csf.scalar_temp2 = ifelse(area_fraction == 0, zero(FT), csf.scalar_temp2)
 
     # Update the coupler field in-place
     @. csf.F_lh += csf.scalar_temp1 * area_fraction
@@ -628,8 +626,7 @@ function FluxCalculator.compute_surface_fluxes!(
             p.snow.snow_cover_fraction .* snow_dest.vapor_flux
         ) .* ρ_liq,
     )
-    @. csf.scalar_temp1 =
-        ifelse(area_fraction == 0, zero(csf.scalar_temp1), csf.scalar_temp1)
+    @. csf.scalar_temp1 = ifelse(area_fraction == 0, zero(FT), csf.scalar_temp1)
     @. csf.F_turb_moisture += csf.scalar_temp1 * area_fraction
 
     # Combine turbulent momentum fluxes from each component of the land model
@@ -640,17 +637,15 @@ function FluxCalculator.compute_surface_fluxes!(
         soil_dest.ρτxz .* (1 .- p.snow.snow_cover_fraction) .+
         p.snow.snow_cover_fraction .* snow_dest.ρτxz,
     )
-    @. csf.scalar_temp1 =
-        ifelse(area_fraction == 0, zero(csf.scalar_temp1), csf.scalar_temp1)
+    @. csf.scalar_temp1 = ifelse(area_fraction == 0, zero(FT), csf.scalar_temp1)
     @. csf.F_turb_ρτxz += csf.scalar_temp1 * area_fraction
 
     Interfacer.remap!(
         csf.scalar_temp1,
         soil_dest.ρτyz .* (1 .- p.snow.snow_cover_fraction) .+
         p.snow.snow_cover_fraction .* snow_dest.ρτyz,
     )
-    @. csf.scalar_temp1 =
-        ifelse(area_fraction == 0, zero(csf.scalar_temp1), csf.scalar_temp1)
+    @. csf.scalar_temp1 = ifelse(area_fraction == 0, zero(FT), csf.scalar_temp1)
     @. csf.F_turb_ρτyz += csf.scalar_temp1 * area_fraction
 
     # Combine the buoyancy flux from each component of the land model
@@ -661,15 +656,13 @@ function FluxCalculator.compute_surface_fluxes!(
         soil_dest.buoy_flux .* (1 .- p.snow.snow_cover_fraction) .+
         p.snow.snow_cover_fraction .* snow_dest.buoy_flux,
     )
-    @. csf.scalar_temp1 =
-        ifelse(area_fraction == 0, zero(csf.scalar_temp1), csf.scalar_temp1)
+    @. csf.scalar_temp1 = ifelse(area_fraction == 0, zero(FT), csf.scalar_temp1)
     @. csf.buoyancy_flux += csf.scalar_temp1 * area_fraction
 
     # Compute ustar from the momentum fluxes and surface air density
     #  ustar = sqrt(ρτ / ρ)
     @. csf.scalar_temp1 = sqrt(sqrt(csf.F_turb_ρτxz^2 + csf.F_turb_ρτyz^2) / csf.ρ_atmos)
-    @. csf.scalar_temp1 =
-        ifelse(area_fraction == 0, zero(csf.scalar_temp1), csf.scalar_temp1)
+    @. csf.scalar_temp1 = ifelse(area_fraction == 0, zero(FT), csf.scalar_temp1)
     # If ustar is zero, set it to eps to avoid division by zero in the atmosphere
     @. csf.ustar += max(csf.scalar_temp1 * area_fraction, eps(FT))
 
@@ -683,8 +676,7 @@ function FluxCalculator.compute_surface_fluxes!(
     surface_params = LP.surface_fluxes_parameters(sim.model.soil.parameters.earth_param_set)
     @. csf.scalar_temp1 =
         -csf.ustar^3 / SFP.von_karman_const(surface_params) / non_zero(csf.buoyancy_flux)
-    @. csf.scalar_temp1 =
-        ifelse(area_fraction == 0, zero(csf.scalar_temp1), csf.scalar_temp1)
+    @. csf.scalar_temp1 = ifelse(area_fraction == 0, zero(FT), csf.scalar_temp1)
     # When L_MO is infinite, avoid multiplication by zero to prevent NaN
     @. csf.L_MO +=
         ifelse(isinf(csf.scalar_temp1), csf.scalar_temp1, csf.scalar_temp1 * area_fraction)

experiments/ClimaEarth/components/ocean/oceananigans.jl

@@ -238,9 +238,9 @@ function FieldExchanger.resolve_area_fractions!(
         polar_mask .= abs.(lat) .>= FT(80)
 
         # 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)
+        @. land_fraction = ifelse(polar_mask == FT(1), one(FT), land_fraction)
+        @. ice_fraction = ifelse(polar_mask == FT(1), zero(FT), ice_fraction)
+        @. ocean_fraction = ifelse(polar_mask == FT(1), zero(FT), ocean_fraction)
     end
 
     # Update the ice concentration field in the ocean simulation

experiments/ClimaEarth/components/ocean/prescr_seaice.jl

@@ -177,7 +177,7 @@ function PrescribedIceSimulation(
     evaluate!(ice_fraction, SIC_timevaryinginput, tspan[1])
 
     # Make ice fraction binary rather than fractional
-    ice_fraction = ifelse.(ice_fraction .> FT(0.5), FT(1), FT(0))
+    @. ice_fraction = ifelse(ice_fraction > FT(0.5), one(FT), zero(FT))
 
     params = IceSlabParameters{FT}(coupled_param_dict)
 
@@ -322,7 +322,7 @@ 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))
+    @. p.area_fraction = ifelse(p.area_fraction > FT(0.5), one(FT), zero(FT))
 
     # 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)
@@ -337,7 +337,7 @@ function ice_rhs!(dY, Y, p, t)
         F_conductive
     ) / (h * ρ * 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)
+    @. rhs = ifelse(p.area_fraction ≈ 0, zero(FT), rhs)
 
     # If tendencies lead to temperature above freezing, set temperature to freezing
     @. dY.T_bulk = min(rhs, (T_freeze - Y.T_bulk) / float(p.dt))

experiments/ClimaEarth/components/ocean/slab_ocean.jl

@@ -259,10 +259,11 @@ end
 function slab_ocean_rhs!(dY, Y, cache, t)
     params, F_turb_energy, SW_d, LW_d = cache
     (; α, ϵ, σ, h, ρ, c) = params
+    FT = eltype(Y)
     rhs = @. (-F_turb_energy + (1 - α) * SW_d + ϵ * (LW_d - σ * Y.T_sfc^4)) / (h * ρ * c)
 
     # Zero out tendencies where there is no ocean, so that temperature remains constant there
-    @. rhs = ifelse(cache.area_fraction ≈ 0, zero(rhs), rhs)
+    @. rhs = ifelse(cache.area_fraction ≈ 0, zero(FT), rhs)
 
     # Note that the area fraction has already been applied to the fluxes,
     #  so we don't need to multiply by it here.

experiments/ClimaEarth/setup_run.jl

@@ -269,7 +269,7 @@ function CoupledSimulation(config_dict::AbstractDict)
     land_mask_data =
         joinpath(@clima_artifact("landsea_mask_60arcseconds", comms_ctx), "landsea_mask.nc")
     land_fraction = SpaceVaryingInput(land_mask_data, "landsea", boundary_space)
-    land_fraction = ifelse.(land_fraction .> eps(FT), FT(1), FT(0))
+    @. land_fraction = ifelse(land_fraction > eps(FT), one(FT), zero(FT))
 
     #=
     ### Surface Models: AMIP and SlabPlanet Modes

experiments/ClimaEarth/test/compare.jl

@@ -26,7 +26,7 @@ function _error(arr1::AbstractArray, arr2::AbstractArray; ABS_TOL = 100eps(eltyp
     arr2 = Array(arr2) .* isfinite.(Array(arr2))
     diff = abs.(arr1 .- arr2)
     denominator = abs.(arr1)
-    error = ifelse.(denominator .> ABS_TOL, diff ./ denominator, diff)
+    error = @. ifelse(denominator > ABS_TOL, diff ./ denominator, diff)
     return error
 end
 

experiments/ClimaEarth/test/fluxes_test.jl

@@ -29,6 +29,7 @@ include(joinpath("..", "setup_run.jl"))
     cs = CoupledSimulation(config_dict)
     step!(cs)
     boundary_space = Interfacer.boundary_space(cs)
+    FT = CC.Spaces.undertype(boundary_space)
 
     # Unpack component models
     (; atmos_sim, land_sim, ocean_sim, ice_sim) = cs.model_sims
@@ -46,10 +47,10 @@ include(joinpath("..", "setup_run.jl"))
     p = land_sim.integrator.p
     land_fraction = Interfacer.get_field(land_sim, Val(:area_fraction))
     land_flux = Interfacer.remap(land_sim.integrator.p.bucket.R_n, boundary_space)
-    @. land_flux = ifelse(land_fraction ≈ 0, zero(land_flux), land_flux)
+    @. land_flux = ifelse(land_fraction ≈ 0, zero(FT), land_flux)
 
     err_land = @. atmos_flux - land_flux
-    @. err_land = ifelse(land_fraction ≈ 0, zero(err_land), err_land)
+    @. err_land = ifelse(land_fraction ≈ 0, zero(FT), err_land)
     @show "Bucket flux error: $(maximum(abs.(err_land)))"
     @test maximum(abs.(err_land)) < 5
 
@@ -65,7 +66,7 @@ include(joinpath("..", "setup_run.jl"))
     ice_rad_flux =
         (1 .- α) .* p.SW_d .+
         ϵ .* (p.LW_d .- σ .* Interfacer.get_field(ice_sim, Val(:surface_temperature)) .^ 4)
-    @. ice_rad_flux = ifelse(p.area_fraction ≈ 0, zero(ice_rad_flux), ice_rad_flux)
+    @. ice_rad_flux = ifelse(p.area_fraction ≈ 0, zero(FT), ice_rad_flux)
     ice_fraction = Interfacer.get_field(ice_sim, Val(:area_fraction))
 
     # Prescribed ocean: SST is prescribed, but for this test we can still compute
@@ -79,7 +80,7 @@ include(joinpath("..", "setup_run.jl"))
             σ .* Interfacer.get_field(ocean_sim, Val(:surface_temperature)) .^ 4
         )
     ocean_fraction = Interfacer.get_field(ocean_sim, Val(:area_fraction))
-    @. ocean_rad_flux = ifelse(ocean_fraction ≈ 0, zero(ocean_rad_flux), ocean_rad_flux)
+    @. ocean_rad_flux = ifelse(ocean_fraction ≈ 0, zero(FT), ocean_rad_flux)
 
     # Combine component fluxes by area-weighted sum (incl. bucket sign convention):
     combined_fluxes =

src/FieldExchanger.jl

@@ -38,7 +38,7 @@ function update_surface_fractions!(cs::Interfacer.CoupledSimulation)
     if haskey(cs.model_sims, :land_sim)
         land_fraction = Interfacer.get_field(cs.model_sims.land_sim, Val(:area_fraction))
     else
-        cs.fields.scalar_temp1 .= 0
+        cs.fields.scalar_temp1 .= zero(FT)
         land_fraction = cs.fields.scalar_temp1
     end
 
@@ -56,7 +56,7 @@ function update_surface_fractions!(cs::Interfacer.CoupledSimulation)
         )
         ice_fraction = Interfacer.get_field(ice_sim, Val(:area_fraction))
     else
-        cs.fields.scalar_temp1 .= 0
+        cs.fields.scalar_temp1 .= zero(FT)
         ice_fraction = cs.fields.scalar_temp1
     end
 
@@ -74,7 +74,7 @@ function update_surface_fractions!(cs::Interfacer.CoupledSimulation)
             resolve_area_fractions!(ocean_sim, cs.model_sims.ice_sim, land_fraction)
         end
     else
-        cs.fields.scalar_temp1 .= 0
+        cs.fields.scalar_temp1 .= zero(FT)
         ocean_fraction = cs.fields.scalar_temp1
     end
 
@@ -331,7 +331,8 @@ is computed from the combined upward longwave radiation.
 """
 function combine_surfaces!(combined_field, sims, field_name, scalar_temp)
     boundary_space = axes(combined_field)
-    combined_field .= 0
+    FT = CC.Spaces.undertype(boundary_space)
+    combined_field .= zero(FT)
     for sim in sims
         if sim isa Interfacer.SurfaceModelSimulation
             # Store the area fraction of this simulation in `scalar_temp`
@@ -342,7 +343,7 @@ function combine_surfaces!(combined_field, sims, field_name, scalar_temp)
                 scalar_temp .*
                 ifelse.(
                     scalar_temp .≈ 0,
-                    zero(combined_field),
+                    zero(FT),
                     Interfacer.get_field(sim, field_name, boundary_space),
                 )
         end
@@ -357,7 +358,7 @@ function combine_surfaces!(csf, sims, field_name::Val{:surface_temperature})
     boundary_space = axes(T_sfc)
     FT = CC.Spaces.undertype(boundary_space)
 
-    T_sfc .= FT(0)
+    T_sfc .= zero(FT)
     for sim in sims
         if sim isa Interfacer.SurfaceModelSimulation
             # Store the area fraction and emissivity of this simulation in temp fields
@@ -373,7 +374,7 @@ function combine_surfaces!(csf, sims, field_name::Val{:surface_temperature})
                 area_fraction .*
                 ifelse.(
                     area_fraction .≈ 0,
-                    zero(T_sfc),
+                    zero(FT),
                     emissivity_sim .*
                     Interfacer.get_field(sim, field_name, boundary_space) .^ FT(4),
                 )

src/FluxCalculator.jl

@@ -332,14 +332,14 @@ function compute_surface_fluxes!(
     # Zero out fluxes where the area fraction is zero
     # Multiplying by `area_fraction` is not sufficient because the fluxes may
     # be NaN where the area fraction is zero.
-    @. F_turb_ρτxz = ifelse(area_fraction ≈ 0, zero(F_turb_ρτxz), F_turb_ρτxz)
-    @. F_turb_ρτyz = ifelse(area_fraction ≈ 0, zero(F_turb_ρτyz), F_turb_ρτyz)
-    @. F_sh = ifelse(area_fraction ≈ 0, zero(F_sh), F_sh)
-    @. F_lh = ifelse(area_fraction ≈ 0, zero(F_lh), F_lh)
-    @. F_turb_moisture = ifelse(area_fraction ≈ 0, zero(F_turb_moisture), F_turb_moisture)
-    @. L_MO = ifelse(area_fraction ≈ 0, zero(L_MO), L_MO)
-    @. ustar = ifelse(area_fraction ≈ 0, zero(ustar), ustar)
-    @. buoyancy_flux = ifelse(area_fraction ≈ 0, zero(buoyancy_flux), buoyancy_flux)
+    @. F_turb_ρτxz = ifelse(area_fraction ≈ 0, zero(FT), F_turb_ρτxz)
+    @. F_turb_ρτyz = ifelse(area_fraction ≈ 0, zero(FT), F_turb_ρτyz)
+    @. F_sh = ifelse(area_fraction ≈ 0, zero(FT), F_sh)
+    @. F_lh = ifelse(area_fraction ≈ 0, zero(FT), F_lh)
+    @. F_turb_moisture = ifelse(area_fraction ≈ 0, zero(FT), F_turb_moisture)
+    @. L_MO = ifelse(area_fraction ≈ 0, zero(FT), L_MO)
+    @. ustar = ifelse(area_fraction ≈ 0, zero(FT), ustar)
+    @. buoyancy_flux = ifelse(area_fraction ≈ 0, zero(FT), buoyancy_flux)
 
     # 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)