diff --git a/experiments/AMIP/modular/components/land/bucket_init.jl b/experiments/AMIP/modular/components/land/bucket_init.jl index 370573712..83dffc75c 100644 --- a/experiments/AMIP/modular/components/land/bucket_init.jl +++ b/experiments/AMIP/modular/components/land/bucket_init.jl @@ -163,19 +163,19 @@ Initializes the bucket model variables. """ function bucket_init( ::Type{FT}, - tspan::Tuple{FT, FT}, + tspan::Tuple{Float64, Float64}, config::String, albedo_type::String, land_temperature_anomaly::String, comms_ctx::AbstractCommsContext, regrid_dirpath::String; space, - dt::FT, - saveat::FT, + dt::Float64, + saveat::Float64, area_fraction, stepper = CTS.RK4(), date_ref::DateTime, - t_start::FT, + t_start::Float64, ) where {FT} if config != "sphere" println( @@ -210,7 +210,7 @@ function bucket_init( z_0b = FT(1e-3) κ_soil = FT(0.7) ρc_soil = FT(2e8) - t_crit = dt # This is the timescale on which snow exponentially damps to zero, in the case where all + t_crit = FT(dt) # This is the timescale on which snow exponentially damps to zero, in the case where all # the snow would melt in time t_crit. It prevents us from having to specially time step in cases where # all the snow melts in a single timestep. params = BucketModelParameters(κ_soil, ρc_soil, albedo, σS_c, W_f, z_0m, z_0b, t_crit, earth_param_set) diff --git a/experiments/AMIP/modular/components/ocean/eisenman_seaice_init.jl b/experiments/AMIP/modular/components/ocean/eisenman_seaice_init.jl index 35488bf2d..21fd472aa 100644 --- a/experiments/AMIP/modular/components/ocean/eisenman_seaice_init.jl +++ b/experiments/AMIP/modular/components/ocean/eisenman_seaice_init.jl @@ -103,13 +103,15 @@ function solve_eisenman_model!(Y, Ya, p, thermo_params, Δt) # ice thickness and mixed layer temperature changes due to atmosphereic and ocean fluxes ice_covered = parent(h_ice)[1] > 0 + + FT = eltype(T_ml) if ice_covered # ice-covered F_base = @. C0_base * (T_ml - T_base) - ΔT_ml = @. -(F_base - ocean_qflux) * Δt / (hρc_ml) - Δh_ice = @. (F_atm - F_base - ocean_qflux) * Δt / L_ice - @. e_base .+= F_base * Δt + ΔT_ml = @. -(F_base - ocean_qflux) * FT(Δt) / (hρc_ml) + Δh_ice = @. (F_atm - F_base - ocean_qflux) * FT(Δt) / L_ice + @. e_base .+= F_base * FT(Δt) else # ice-free - ΔT_ml = @. -(F_atm - ocean_qflux) * Δt / (hρc_ml) + ΔT_ml = @. -(F_atm - ocean_qflux) * FT(Δt) / (hρc_ml) Δh_ice = 0 end @@ -317,11 +319,13 @@ function get_field(sim::EisenmanIceSimulation, ::Val{:energy}) e_base = cache.Ya.e_base ocean_qflux = cache.Ya.ocean_qflux + FT = eltype(sim.integrator.u) + hρc_ml = p_o.h * p_o.ρ * p_o.c e_ml = @. p_o.h * p_o.ρ * p_o.c * sim.integrator.u.T_ml # heat e_ice = @. p_i.L_ice * sim.integrator.u.h_ice # phase - e_qflux = @. ocean_qflux * sim.integrator.t + e_qflux = @. ocean_qflux * FT(sim.integrator.t) return @. e_ml + e_ice + e_qflux + e_base diff --git a/experiments/AMIP/modular/components/ocean/prescr_seaice_init.jl b/experiments/AMIP/modular/components/ocean/prescr_seaice_init.jl index 7b9d59688..67c60c696 100644 --- a/experiments/AMIP/modular/components/ocean/prescr_seaice_init.jl +++ b/experiments/AMIP/modular/components/ocean/prescr_seaice_init.jl @@ -79,7 +79,7 @@ function ice_rhs!(du, u, p, _) # do not count tendencies that lead to temperatures above freezing, and mask out no-ice areas area_mask = Regridder.binary_mask.(area_fraction, threshold = eps(FT)) - unphysical = @. Regridder.binary_mask.(T_freeze - (Y.T_sfc + FT(rhs) * p.dt), threshold = FT(0)) .* area_mask + unphysical = @. Regridder.binary_mask.(T_freeze - (Y.T_sfc + FT(rhs) * FT(p.dt)), threshold = FT(0)) .* area_mask parent(dY.T_sfc) .= parent(rhs .* unphysical) @. p.q_sfc = TD.q_vap_saturation_generic.(p.thermo_params, Y.T_sfc, p.ρ_sfc, TD.Ice()) diff --git a/src/ConservationChecker.jl b/src/ConservationChecker.jl index 82fb47d45..55fd8399c 100644 --- a/src/ConservationChecker.jl +++ b/src/ConservationChecker.jl @@ -100,16 +100,16 @@ function check_conservation!( parent(coupler_sim.fields.F_radiative_TOA) .= parent(Interfacer.get_field(sim, Val(:F_radiative_TOA))) if isempty(ccs.toa_net_source) - radiation_sources_accum = sum(coupler_sim.fields.F_radiative_TOA .* coupler_sim.Δt_cpl) # ∫ J / m^2 dA + radiation_sources_accum = sum(coupler_sim.fields.F_radiative_TOA .* FT(coupler_sim.Δt_cpl)) # ∫ J / m^2 dA else radiation_sources_accum = - sum(coupler_sim.fields.F_radiative_TOA .* coupler_sim.Δt_cpl) .+ ccs.toa_net_source[end] # ∫ J / m^2 dA + sum(coupler_sim.fields.F_radiative_TOA .* FT(coupler_sim.Δt_cpl)) .+ ccs.toa_net_source[end] # ∫ J / m^2 dA end push!(ccs.toa_net_source, radiation_sources_accum) # save atmos previous = getproperty(ccs, sim_name) - current = sum(FT.(Interfacer.get_field(sim, Val(:energy)))) # # ∫ J / m^3 dV + current = sum(Interfacer.get_field(sim, Val(:energy))) # # ∫ J / m^3 dV push!(previous, current) total += current + radiation_sources_accum @@ -122,7 +122,7 @@ function check_conservation!( current = FT(0) else previous = getproperty(ccs, sim_name) - current = sum(FT.(Interfacer.get_field(sim, Val(:energy)) .* area_fraction)) # # ∫ J / m^3 dV + current = sum(Interfacer.get_field(sim, Val(:energy)) .* area_fraction) # # ∫ J / m^3 dV end push!(previous, current) total += current @@ -211,7 +211,8 @@ function surface_water_gain_from_rates(cs::Interfacer.CoupledSimulation) evaporation = cs.fields.F_turb_moisture # kg / m^2 / s / layer depth precipitation_l = cs.fields.P_liq precipitation_s = cs.fields.P_snow - @. -(evaporation + precipitation_l + precipitation_s) * cs.Δt_cpl # kg / m^2 / layer depth + FT = eltype(evaporation) + @. -(evaporation + precipitation_l + precipitation_s) * FT(cs.Δt_cpl) # kg / m^2 / layer depth end # setup the GKS socket application environment as nul for better performance and to avoid GKS connection errors while plotting diff --git a/src/Diagnostics.jl b/src/Diagnostics.jl index 269e31795..6104ba599 100644 --- a/src/Diagnostics.jl +++ b/src/Diagnostics.jl @@ -6,7 +6,7 @@ This module contains functions for defining, gathering and outputting online mod module Diagnostics using ClimaCore: Spaces, Fields, InputOutput -using ClimaCoupler.Interfacer: CoupledSimulation +using ClimaCoupler.Interfacer: CoupledSimulation, float_type using Dates using ClimaCoupler.TimeManager: AbstractFrequency, Monthly, EveryTimestep, trigger_callback using ClimaComms @@ -78,12 +78,12 @@ end Collects diagnostics in diags names. """ function collect_diags(cs::CoupledSimulation, dg::DiagnosticsGroup) - + FT = float_type(cs) diags = (;) diag_names = propertynames(dg.field_vector) for name in diag_names - diags = (; diags..., zip((name,), (get_var(cs, Val(name)),))...) + diags = (; diags..., zip((name,), (FT.(get_var(cs, Val(name))),))...) end return Fields.FieldVector(; diags...) diff --git a/test/component_model_tests/eisenman_seaice_tests.jl b/test/component_model_tests/eisenman_seaice_tests.jl index 5bc16cdb2..dd2c9e423 100644 --- a/test/component_model_tests/eisenman_seaice_tests.jl +++ b/test/component_model_tests/eisenman_seaice_tests.jl @@ -80,7 +80,7 @@ for FT in (Float32, Float64) @test all(parent(Ya.e_base) .≈ 0) @test all(parent(Y.T_ml) .≈ params_ice.T_base) # ocean temperature below ice remains unchanged - h_ice_new = @. h_ice_0 + F_atm * Δt / params_ice.L_ice + h_ice_new = @. h_ice_0 + F_atm * FT(Δt) / params_ice.L_ice @test all(parent(Y.h_ice) .≈ parent(h_ice_new)) # ice growth T_sfc_new = @. params_ice.T_base .+ (-F_atm) / (params_ice.k_ice / (h_ice_new) + ∂F_atm∂T_sfc) @test all(parent(Y.T_sfc) .≈ parent(T_sfc_new)) # surface temperature decreases @@ -112,7 +112,7 @@ for FT in (Float32, Float64) @test all(parent(Ya.e_base) .≈ 0) # no contribution from basal fluxes @test all(parent(Y.T_ml) .≈ params_ice.T_base) # ocean temperature stays at freezing point - h_ice_new = @. h_ice_0 + F_atm * Δt / params_ice.L_ice + h_ice_new = @. h_ice_0 + F_atm * FT(Δt) / params_ice.L_ice @test all(parent(Y.h_ice) .≈ parent(h_ice_new)) # no ice growth T_sfc_new = @. params_ice.T_base .+ (-F_atm) / (params_ice.k_ice / (h_ice_new) + ∂F_atm∂T_sfc) @test all(parent(Y.T_sfc) .≈ parent(T_sfc_new)) # surface temperature decreases @@ -150,7 +150,7 @@ for FT in (Float32, Float64) @test all(parent(Ya.e_base) .≈ 0) # no contribution from basal fluxes @test all(parent(Y.T_ml) .≈ params_ice.T_base) # ocean temperature stays at freezing point - h_ice_new = @. h_ice_0 + F_atm * Δt / params_ice.L_ice + h_ice_new = @. h_ice_0 + F_atm * FT(Δt) / params_ice.L_ice @test all(parent(Y.h_ice) .≈ parent(h_ice_new)) # ice growth @test all(parent(Y.T_sfc) .≈ params_ice.T_base) # ice surface temperature doesn't exceed freezing @test all(parent(Y.T_sfc) .< parent(T_sfc_0 .+ ΔT_sfc)) @@ -180,7 +180,7 @@ for FT in (Float32, Float64) F_atm = @. Ya.F_rad + Ya.F_turb @test all(parent(Y.h_ice) .≈ 0) # ice melts - T_ml_new = @. T_ml_0 - (F_atm) * Δt / (params_ocean.h * params_ocean.ρ * params_ocean.c) - + T_ml_new = @. T_ml_0 - (F_atm) * FT(Δt) / (params_ocean.h * params_ocean.ρ * params_ocean.c) - h_ice_0 * params_ice.L_ice / (params_ocean.h * params_ocean.ρ * params_ocean.c) @test all(parent(Y.T_ml) .≈ parent(T_ml_new)) # ocean temperature increases due to F_atm (reduced by the latent heat of melting) @@ -210,7 +210,7 @@ for FT in (Float32, Float64) F_atm = @. Ya.F_rad + Ya.F_turb @test all(parent(Ya.e_base) .≈ 0) # no contribution from basal fluxes - T_ml_new = @. T_ml_0 - (F_atm) * Δt / (params_ocean.h * params_ocean.ρ * params_ocean.c) + T_ml_new = @. T_ml_0 - (F_atm) * FT(Δt) / (params_ocean.h * params_ocean.ρ * params_ocean.c) @test all(parent(Y.T_ml) .≈ parent(T_ml_new)) # ocean temperature increases @test all(parent(Y.h_ice) .≈ 0) # no ice @test all(parent(Y.T_sfc) .≈ parent(T_ml_new)) # surface temperature = ocean temperature @@ -240,10 +240,11 @@ for FT in (Float32, Float64) F_atm = @. Ya.F_rad + Ya.F_turb ∂F_atm∂T_sfc = get_∂F_rad_energy∂T_sfc(T_sfc_0, params_ice) .+ Ya.∂F_turb_energy∂T_sfc + FT = eltype(F_atm) @test all(parent(Ya.e_base) .≈ 0) # no contribution from basal fluxes @test all(parent(Y.T_ml) .≈ params_ice.T_freeze) # ocean temperature remains at the freezing point h_ice_new = - @. (F_atm - (T_ml_0 - params_ice.T_freeze) * params_ocean.h * params_ocean.ρ * params_ocean.c) * Δt / + @. (F_atm - (T_ml_0 - params_ice.T_freeze) * params_ocean.h * params_ocean.ρ * params_ocean.c) * FT(Δt) / params_ice.L_ice @test all(parent(Y.h_ice) .≈ parent(h_ice_new)) # ice growth @test all(parent(Y.h_ice) .> 0) @@ -277,11 +278,11 @@ for FT in (Float32, Float64) solve_eisenman_model!(Y[colidx], Ya[colidx], params, thermo_params, Δt) end - @test all(parent(Ya.e_base) .≈ params_ice.C0_base * ΔT_ml * Δt) # non-zero contribution from basal flux - T_ml_new = T_ml_0 .- params_ice.C0_base .* ΔT_ml * Δt / (params_ocean.h * params_ocean.ρ * params_ocean.c) + @test all(parent(Ya.e_base) .≈ params_ice.C0_base * ΔT_ml * FT(Δt)) # non-zero contribution from basal flux + T_ml_new = T_ml_0 .- params_ice.C0_base .* ΔT_ml * FT(Δt) / (params_ocean.h * params_ocean.ρ * params_ocean.c) @test all(parent(Y.T_ml) .≈ parent(T_ml_new)) # ocean temperature decreases @test all(parent(Y.T_ml) .< parent(T_ml_0)) - h_ice_new = @. h_ice_0 - params_ice.C0_base * 10 * Δt / params_ice.L_ice + h_ice_new = @. h_ice_0 - params_ice.C0_base * 10 * FT(Δt) / params_ice.L_ice @test all(parent(Y.h_ice) .≈ parent(h_ice_new)) @test all(parent(Y.T_sfc) .≈ params_ice.T_base) # surface temperature unchanged (T_ml can affect in only if ice free) end @@ -312,11 +313,11 @@ for FT in (Float32, Float64) end @test all(parent(Ya.e_base) .≈ 0) # no contribution from basal flux - T_ml_new = @. T_ml_0 + Ya.ocean_qflux * Δt / (params_ocean.h * params_ocean.ρ * params_ocean.c) + T_ml_new = @. T_ml_0 + Ya.ocean_qflux * FT(Δt) / (params_ocean.h * params_ocean.ρ * params_ocean.c) @test all(parent(Y.T_ml) .≈ parent(T_ml_new)) # qflux increases ocean temperature @test all(parent(Y.T_ml) .> parent(T_ml_0)) @test all(parent(Y.T_sfc) .≈ params_ice.T_base) # surface temperature unchanged (T_ml can affect in only if ice free) - h_ice_new = @. h_ice_0 - Ya.ocean_qflux * Δt / params_ice.L_ice + h_ice_new = @. h_ice_0 - Ya.ocean_qflux * FT(Δt) / params_ice.L_ice @test all(parent(Y.h_ice) .≈ parent(h_ice_new)) # # qflux decreases sea-ice thickness @test all(parent(Y.h_ice) .< parent(h_ice_0)) end @@ -347,10 +348,10 @@ for FT in (Float32, Float64) @test all(parent(h_ice) .≈ 0.001) step!(sim, 2 * Δt) h_ice = sim.integrator.u.h_ice - @test all(parent(h_ice) .≈ 0.002) + @test all(abs.(parent(h_ice) .- 0.002) .< 10eps(FT)) total_energy_calc = (get_field(sim, Val(:energy)) .- total_energy_0) - total_energy_expeted = 300 .* ones(boundary_space) .* 2 .* Δt + total_energy_expeted = 300 .* ones(boundary_space) .* 2 .* FT(Δt) @test all(parent(total_energy_calc) .≈ parent(total_energy_expeted)) end diff --git a/test/conservation_checker_tests.jl b/test/conservation_checker_tests.jl index 4a8f3f51d..cf072ace8 100644 --- a/test/conservation_checker_tests.jl +++ b/test/conservation_checker_tests.jl @@ -25,51 +25,36 @@ struct TestAtmos{I} <: Interfacer.AtmosModelSimulation i::I end name(s::TestAtmos) = "TestAtmos" -function get_field(s::TestAtmos, ::Val{:F_radiative_TOA}) - FT = Domains.float_type(Meshes.domain(s.i.space.grid.topology.mesh)) - ones(s.i.space) .* FT.(200) -end +get_field(s::TestAtmos, ::Val{:F_radiative_TOA}) = ones(s.i.space) .* 200 +get_field(s::TestAtmos, ::Val{:water}) = ones(s.i.space) .* 1 function get_field(s::TestAtmos, ::Val{:energy}) FT = Domains.float_type(Meshes.domain(s.i.space.grid.topology.mesh)) - ones(s.i.space) .* FT.(1e6) -end -function get_field(s::TestAtmos, ::Val{:water}) - FT = Domains.float_type(Meshes.domain(s.i.space.grid.topology.mesh)) - ones(s.i.space) .* FT.(1) + ones(s.i.space) .* FT(1e6) end struct TestOcean{I} <: Interfacer.SurfaceModelSimulation i::I end name(s::TestOcean) = "TestOcean" +get_field(s::TestOcean, ::Val{:water}) = ones(s.i.space) .* 0 function get_field(s::TestOcean, ::Val{:energy}) FT = Domains.float_type(Meshes.domain(s.i.space.grid.topology.mesh)) - ones(s.i.space) .* FT.(1e6) -end -function get_field(s::TestOcean, ::Val{:water}) - FT = Domains.float_type(Meshes.domain(s.i.space.grid.topology.mesh)) - ones(s.i.space) .* FT.(0) + ones(s.i.space) .* FT(1e6) end function get_field(s::TestOcean, ::Val{:area_fraction}) FT = Domains.float_type(Meshes.domain(s.i.space.grid.topology.mesh)) - ones(s.i.space) .* FT.(0.25) + ones(s.i.space) .* FT(0.25) end struct TestLand{I} <: Interfacer.SurfaceModelSimulation i::I end name(s::TestLand) = "TestLand" -function get_field(s::TestLand, ::Val{:energy}) - FT = Domains.float_type(Meshes.domain(s.i.space.grid.topology.mesh)) - ones(s.i.space) .* FT.(0) -end -function get_field(s::TestLand, ::Val{:water}) - FT = Domains.float_type(Meshes.domain(s.i.space.grid.topology.mesh)) - ones(s.i.space) .* FT.(0) -end +get_field(s::TestLand, ::Val{:energy}) = ones(s.i.space) .* 0 +get_field(s::TestLand, ::Val{:water}) = ones(s.i.space) .* 0 function get_field(s::TestLand, ::Val{:area_fraction}) FT = Domains.float_type(Meshes.domain(s.i.space.grid.topology.mesh)) - ones(s.i.space) .* FT.(0.25) + ones(s.i.space) .* FT(0.25) end diff --git a/test/diagnostics_tests.jl b/test/diagnostics_tests.jl index 2ba04fc9c..0435e77be 100644 --- a/test/diagnostics_tests.jl +++ b/test/diagnostics_tests.jl @@ -19,7 +19,7 @@ import ClimaCoupler.Diagnostics: post_save, save_time_format -get_var(cs::Interfacer.CoupledSimulation, ::Val{:x}) = float_type(cs)(1) +get_var(cs::Interfacer.CoupledSimulation, ::Val{:x}) = 1 for FT in (Float32, Float64) @testset "init_diagnostics for FT=$FT" begin diff --git a/test/field_exchanger_tests.jl b/test/field_exchanger_tests.jl index 8fcad85cf..5eb57521c 100644 --- a/test/field_exchanger_tests.jl +++ b/test/field_exchanger_tests.jl @@ -74,7 +74,10 @@ end struct TestSurfaceSimulationLand{C} <: SurfaceModelSimulation cache::C end -get_field(::TestSurfaceSimulationLand, ::Val{:area_fraction}) = 0.5 +function get_field(sim::TestSurfaceSimulationLand, ::Val{:area_fraction}) + FT = eltype(sim.cache.turbulent_energy_flux) + return FT(0.5) +end function update_field!(sim::TestSurfaceSimulationLand, ::Val{:turbulent_energy_flux}, field) parent(sim.cache.turbulent_energy_flux) .= parent(field) end @@ -110,7 +113,6 @@ for FT in (Float32, Float64) end end - @testset "import_combined_surface_fields! for FT=$FT" begin # coupler cache setup boundary_space = TestHelper.create_space(FT) diff --git a/test/interfacer_tests.jl b/test/interfacer_tests.jl index e1f503d0c..d8f9dc622 100644 --- a/test/interfacer_tests.jl +++ b/test/interfacer_tests.jl @@ -28,10 +28,10 @@ struct DummySimulation3{S} <: LandModelSimulation end get_field(sim::SurfaceModelSimulation, ::Val{:var}) = ones(sim.space) -get_field(sim::SurfaceModelSimulation, ::Val{:surface_temperature}) = ones(sim.space) .* FT(300) +get_field(sim::SurfaceModelSimulation, ::Val{:surface_temperature}) = ones(sim.space) .* 300 function get_field(sim::SurfaceModelSimulation, ::Val{:var_float}) FT = Domains.float_type(Meshes.domain(sim.space.grid.topology.mesh)) - FT(2) + return FT(2) end for FT in (Float32, Float64)