diff --git a/experiments/ClimaEarth/components/land/climaland_bucket.jl b/experiments/ClimaEarth/components/land/climaland_bucket.jl index 107c045dbc..90e5e20fe7 100644 --- a/experiments/ClimaEarth/components/land/climaland_bucket.jl +++ b/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: diff --git a/experiments/ClimaEarth/components/land/climaland_integrated.jl b/experiments/ClimaEarth/components/land/climaland_integrated.jl index 479afaf233..05b7f86aed 100644 --- a/experiments/ClimaEarth/components/land/climaland_integrated.jl +++ b/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,8 +637,7 @@ 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!( @@ -649,8 +645,7 @@ function FluxCalculator.compute_surface_fluxes!( 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) diff --git a/experiments/ClimaEarth/components/ocean/oceananigans.jl b/experiments/ClimaEarth/components/ocean/oceananigans.jl index c90d477a4d..a5e438440f 100644 --- a/experiments/ClimaEarth/components/ocean/oceananigans.jl +++ b/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 diff --git a/experiments/ClimaEarth/components/ocean/prescr_seaice.jl b/experiments/ClimaEarth/components/ocean/prescr_seaice.jl index fecd8a224b..ec4772b65f 100644 --- a/experiments/ClimaEarth/components/ocean/prescr_seaice.jl +++ b/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)) diff --git a/experiments/ClimaEarth/components/ocean/slab_ocean.jl b/experiments/ClimaEarth/components/ocean/slab_ocean.jl index a2efe94705..94eb2e1334 100644 --- a/experiments/ClimaEarth/components/ocean/slab_ocean.jl +++ b/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. diff --git a/experiments/ClimaEarth/setup_run.jl b/experiments/ClimaEarth/setup_run.jl index 5daf5611d6..f067774679 100644 --- a/experiments/ClimaEarth/setup_run.jl +++ b/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 diff --git a/experiments/ClimaEarth/test/compare.jl b/experiments/ClimaEarth/test/compare.jl index 3b63b0d12b..4c34577adc 100644 --- a/experiments/ClimaEarth/test/compare.jl +++ b/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 diff --git a/experiments/ClimaEarth/test/fluxes_test.jl b/experiments/ClimaEarth/test/fluxes_test.jl index a4e062c57b..6943745a71 100644 --- a/experiments/ClimaEarth/test/fluxes_test.jl +++ b/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 = diff --git a/src/FieldExchanger.jl b/src/FieldExchanger.jl index d3ec982352..07533cf905 100644 --- a/src/FieldExchanger.jl +++ b/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), ) diff --git a/src/FluxCalculator.jl b/src/FluxCalculator.jl index f8c8458a80..57494ddd6c 100644 --- a/src/FluxCalculator.jl +++ b/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)