Merge pull request #1497 from CliMA/js/t-radfluxes

compute T_sfc from longwave fluxes for radiation

Author: juliasloan25

experiments/ClimaEarth/components/atmosphere/climaatmos.jl

@@ -399,12 +399,9 @@ function Interfacer.update_field!(sim::ClimaAtmosSimulation, ::Val{:emissivity},
     # Remap field onto the atmosphere surface space in scratch field
     temp_field_surface = sim.integrator.p.scratch.ᶠtemp_field_level
     Interfacer.remap!(temp_field_surface, field)
-    # sets all 16 bands (rows) to the same values
-    sim.integrator.p.radiation.rrtmgp_model.surface_emissivity .= reshape(
-        CC.Fields.field2array(temp_field_surface),
-        1,
-        length(parent(temp_field_surface)),
-    )
+    # Set each row (band) of the emissivity matrix by transposing the vector returned from `field2array`
+    sim.integrator.p.radiation.rrtmgp_model.surface_emissivity .=
+        CC.Fields.field2array(temp_field_surface)'
 end
 function Interfacer.update_field!(
     sim::ClimaAtmosSimulation,
@@ -415,11 +412,9 @@ function Interfacer.update_field!(
     # Remap field onto the atmosphere surface space in scratch field
     temp_field_surface = sim.integrator.p.scratch.ᶠtemp_field_level
     Interfacer.remap!(temp_field_surface, field)
-    sim.integrator.p.radiation.rrtmgp_model.direct_sw_surface_albedo .= reshape(
-        CC.Fields.field2array(temp_field_surface),
-        1,
-        length(parent(temp_field_surface)),
-    )
+    # Set each row (band) of the albedo matrix by transposing the vector returned from `field2array`
+    sim.integrator.p.radiation.rrtmgp_model.direct_sw_surface_albedo .=
+        CC.Fields.field2array(temp_field_surface)'
 end
 
 function Interfacer.update_field!(
@@ -431,11 +426,9 @@ function Interfacer.update_field!(
     # Remap field onto the atmosphere surface space in scratch field
     temp_field_surface = sim.integrator.p.scratch.ᶠtemp_field_level
     Interfacer.remap!(temp_field_surface, field)
-    sim.integrator.p.radiation.rrtmgp_model.diffuse_sw_surface_albedo .= reshape(
-        CC.Fields.field2array(temp_field_surface),
-        1,
-        length(parent(temp_field_surface)),
-    )
+    # Set each row (band) of the albedo matrix by transposing the vector returned from `field2array`
+    sim.integrator.p.radiation.rrtmgp_model.diffuse_sw_surface_albedo .=
+        CC.Fields.field2array(temp_field_surface)'
 end
 
 function Interfacer.update_field!(

experiments/ClimaEarth/components/land/climaland_integrated.jl

@@ -295,6 +295,9 @@ function ClimaLandSimulation(
         # Initialize the surface temperature so the atmosphere can compute radiation.
         @. p.T_sfc = orog_adjusted_T_surface
     end
+    # Initialize the surface emissivity so the atmosphere can compute radiation.
+    # Otherwise, it's initialized to 0 which causes NaNs in the radiation calculation.
+    @. p.ϵ_sfc = FT(1)
 
     # Update cos(zenith angle) within land model every hour
     update_dt = dt isa ITime ? ITime(3600) : 3600

experiments/ClimaEarth/components/ocean/climaocean.jl

@@ -252,6 +252,7 @@ Interfacer.get_field(sim::OceananigansSimulation, ::Val{:roughness_buoyancy}) =
 Interfacer.get_field(sim::OceananigansSimulation, ::Val{:roughness_momentum}) =
     Float32(5.8e-5)
 Interfacer.get_field(sim::OceananigansSimulation, ::Val{:beta}) = Float32(1)
+Interfacer.get_field(sim::OceananigansSimulation, ::Val{:emissivity}) = Float32(0.97)
 Interfacer.get_field(sim::OceananigansSimulation, ::Val{:surface_direct_albedo}) =
     Float32(0.06)
 Interfacer.get_field(sim::OceananigansSimulation, ::Val{:surface_diffuse_albedo}) =

experiments/ClimaEarth/components/ocean/oceananigans.jl

@@ -129,6 +129,8 @@ Interfacer.get_field(sim::PrescribedOceanSimulation, ::Val{:surface_direct_albed
     sim.cache.α_direct
 Interfacer.get_field(sim::PrescribedOceanSimulation, ::Val{:surface_diffuse_albedo}) =
     sim.cache.α_diffuse
+Interfacer.get_field(sim::PrescribedOceanSimulation, ::Val{:emissivity}) =
+    eltype(sim.cache.T_sfc)(1)
 
 function Interfacer.update_field!(
     sim::PrescribedOceanSimulation,

experiments/ClimaEarth/components/ocean/prescr_ocean.jl

@@ -121,19 +121,34 @@ Updates the coupler with the surface properties. The `Interfacer.get_field`
 functions for (`:surface_temperature`, `:surface_direct_albedo`,
 `:surface_diffuse_albedo`) need to be specified for each surface model.
 
-Note: The calculation of surface humidity uses atmospheric properties stored in
+Note: The calculation of surface humidity done here uses atmospheric properties stored in
 the coupled fields. For these values to be correct, this function should be called
 after `import_atmos_fields!` in a timestep.
 
+Note 2: Not all surface fields are imported here. Some quantities are retrieved
+from each surface model when surface fluxes are computed, in `compute_surface_fluxes!`.
+
 # Arguments
 - `csf`: [NamedTuple] containing coupler fields.
 - `model_sims`: [NamedTuple] containing `ComponentModelSimulation`s.
 - `thermo_params`: [TD.Parameters.ThermodynamicsParameters] the thermodynamic parameters.
 """
 function import_combined_surface_fields!(csf, model_sims, thermo_params)
-    combine_surfaces!(csf.T_sfc, model_sims, Val(:surface_temperature))
-    combine_surfaces!(csf.surface_direct_albedo, model_sims, Val(:surface_direct_albedo))
-    combine_surfaces!(csf.surface_diffuse_albedo, model_sims, Val(:surface_diffuse_albedo))
+    combine_surfaces!(csf.emissivity, model_sims, Val(:emissivity), csf.temp1)
+    combine_surfaces!(
+        csf.surface_direct_albedo,
+        model_sims,
+        Val(:surface_direct_albedo),
+        csf.temp1,
+    )
+    combine_surfaces!(
+        csf.surface_diffuse_albedo,
+        model_sims,
+        Val(:surface_diffuse_albedo),
+        csf.temp1,
+    )
+    # Temperature requires emissivity, so we provide all the coupler fields
+    combine_surfaces!(csf, model_sims, Val(:surface_temperature))
 
     # q_sfc is computed from the atmosphere state and surface temperature, so it's handled differently
     # This is computed on the exchange grid, so there's no need to remap
@@ -181,7 +196,8 @@ end
 """
     update_sim!(atmos_sim::Interfacer.AtmosModelSimulation, csf)
 
-Updates the surface fields for temperature, roughness length, albedo, and specific humidity.
+Updates the atmosphere's fields for surface direct and diffuse albedos, emissivity, and temperature,
+as well as the turbulent fluxes.
 
 # Arguments
 - `atmos_sim`: [Interfacer.AtmosModelSimulation] containing an atmospheric model simulation object.
@@ -198,6 +214,7 @@ function update_sim!(atmos_sim::Interfacer.AtmosModelSimulation, csf)
         Val(:surface_diffuse_albedo),
         csf.surface_diffuse_albedo,
     )
+    Interfacer.update_field!(atmos_sim, Val(:emissivity), csf.emissivity)
     Interfacer.update_field!(atmos_sim, Val(:surface_temperature), csf)
     Interfacer.update_field!(atmos_sim, Val(:turbulent_fluxes), csf)
     return nothing
@@ -270,36 +287,79 @@ function step_model_sims!(cs::Interfacer.CoupledSimulation)
 end
 
 """
-    combine_surfaces!(combined_field::CC.Fields.Field, sims, field_name::Val)
+    combine_surfaces!(combined_field::CC.Fields.Field, sims, field_name::Val, temp1)
 
 Sums the fields, specified by `field_name`, weighted by the respective area fractions of all
 surface simulations. THe result is saved in `combined_field`.
 
+For surface temperature, upward longwave radiation is computed from the temperatures
+of each surface, weighted by their area fractions, and then the combined temperature
+is computed from the combined upward longwave radiation.
+
 # Arguments
 - `combined_field`: [CC.Fields.Field] output object containing weighted values.
+    Note: For the surface temperature, all coupler fields are passed in a NamedTuple.
 - `sims`: [NamedTuple] containing simulations .
 - `field_name`: [Val] containing the name Symbol of the field t be extracted by the `Interfacer.get_field` functions.
+- `temp1`: [CC.Fields.Field] temporary field for intermediate calculations.
+    Omitted for surface temperature method.
 
 # Example
 - `combine_surfaces!(temp_field, cs.model_sims, Val(:surface_temperature))`
 """
-function combine_surfaces!(combined_field, sims, field_name)
+function combine_surfaces!(combined_field, sims, field_name, temp1)
     boundary_space = axes(combined_field)
     combined_field .= 0
     for sim in sims
         if sim isa Interfacer.SurfaceModelSimulation
-            # Zero out the contribution from this surface if the area fraction is zero
+            # Store the area fraction of this simulation in `temp1`
+            Interfacer.get_field!(temp1, sim, Val(:area_fraction))
+            # Zero out the contribution from this surface if the area fraction is zero.
             # Note that multiplying by `area_fraction` is not sufficient in the case of NaNs
-            area_fraction = Interfacer.get_field(sim, Val(:area_fraction))
             combined_field .+=
-                area_fraction .*
+                temp1 .*
                 ifelse.(
-                    area_fraction .≈ 0,
+                    temp1 .≈ 0,
                     zero(combined_field),
                     Interfacer.get_field(sim, field_name, boundary_space),
                 )
         end
     end
+    return nothing
+end
+function combine_surfaces!(csf, sims, field_name::Val{:surface_temperature})
+    # extract the coupler fields we need to get the surface temperature
+    T_sfc = csf.T_sfc
+    emissivity_sfc = csf.emissivity
+
+    boundary_space = axes(T_sfc)
+    FT = CC.Spaces.undertype(boundary_space)
+
+    T_sfc .= FT(0)
+    for sim in sims
+        if sim isa Interfacer.SurfaceModelSimulation
+            # Store the area fraction and emissivity of this simulation in temp fields
+            Interfacer.get_field!(csf.temp1, sim, Val(:area_fraction))
+            area_fraction = csf.temp1
+            Interfacer.get_field!(csf.temp2, sim, Val(:emissivity))
+            emissivity_sim = csf.temp2
+
+            # Zero out the contribution from this surface if the area fraction is zero.
+            # Note that multiplying by `area_fraction` is not sufficient in the case of NaNs
+            # Compute upward longwave radiation from surface temperature for this simulation
+            T_sfc .+=
+                area_fraction .*
+                ifelse.(
+                    area_fraction .≈ 0,
+                    zero(T_sfc),
+                    emissivity_sim .*
+                    Interfacer.get_field(sim, field_name, boundary_space) .^ FT(4),
+                )
+        end
+    end
+    # Convert the combined upward longwave radiation into a surface temperature
+    @. T_sfc = (T_sfc / emissivity_sfc)^FT(1 / 4)
+    return nothing
 end
 
 """

src/FieldExchanger.jl

@@ -52,11 +52,15 @@ Interfacer.get_field(
     sim::Union{TestSurfaceSimulation1, TestSurfaceSimulation2},
     ::Val{:beta},
 ) = sim.cache_field
+Interfacer.get_field(
+    sim::Union{TestSurfaceSimulation1, TestSurfaceSimulation2},
+    ::Val{:emissivity},
+) = eltype(sim.cache_field)(1)
 
 Interfacer.get_field(sim::TestSurfaceSimulation1, ::Val{:area_fraction}) =
-    sim.cache_field .* 0
+    sim.cache_field .* CC.Spaces.undertype(axes(sim.cache_field))(0.25)
 Interfacer.get_field(sim::TestSurfaceSimulation2, ::Val{:area_fraction}) =
-    sim.cache_field .* CC.Spaces.undertype(axes(sim.cache_field))(0.5)
+    sim.cache_field .* CC.Spaces.undertype(axes(sim.cache_field))(0.75)
 
 Interfacer.step!(::TestSurfaceSimulation1, _) = nothing
 
@@ -145,6 +149,8 @@ Interfacer.get_field(sim::TestSurfaceSimulationLand, ::Val{:surface_diffuse_albe
     sim.cache.albedo_diffuse
 Interfacer.get_field(sim::TestSurfaceSimulationLand, ::Val{:surface_temperature}) =
     sim.cache.surface_temperature
+Interfacer.get_field(sim::TestSurfaceSimulationLand, ::Val{:emissivity}) =
+    eltype(sim.cache.surface_temperature)(1)
 function Interfacer.update_field!(
     sim::TestSurfaceSimulationLand,
     ::Val{:turbulent_energy_flux},
@@ -266,8 +272,9 @@ for FT in (Float32, Float64)
             Interfacer.get_field(sims.c, var_name),
             Interfacer.get_field(sims.d, var_name),
         )
+        temp_field = CC.Fields.zeros(test_space)
 
-        FieldExchanger.combine_surfaces!(combined_field, sims, var_name)
+        FieldExchanger.combine_surfaces!(combined_field, sims, var_name, temp_field)
         @test combined_field == fill(FT(sum(fractions .* fields)), test_space)
     end
 
@@ -322,30 +329,31 @@ for FT in (Float32, Float64)
         coupler_names_additional = [:surface_direct_albedo, :surface_diffuse_albedo]
         coupler_names =
             push!(Interfacer.default_coupler_fields(), coupler_names_additional...)
-
-        # coupler cache setup
-        exchanged_fields = (
-            :surface_temperature,
-            :surface_direct_albedo,
-            :surface_diffuse_albedo,
-            :roughness_momentum,
-            :roughness_buoyancy,
-            :beta,
-        )
+        coupler_fields = Interfacer.init_coupler_fields(FT, coupler_names, boundary_space)
 
         sims = (;
             a = TestSurfaceSimulation1(ones(boundary_space)),
             b = TestSurfaceSimulation2(ones(boundary_space)),
         )
 
-        coupler_fields = Interfacer.init_coupler_fields(FT, coupler_names, boundary_space)
-
         thermo_params = TDP.ThermodynamicsParameters(FT)
         FieldExchanger.import_combined_surface_fields!(coupler_fields, sims, thermo_params)
+
+        # Analytically compute expected values and compare
+        expected_field_temp =
+            (
+                Interfacer.get_field(sims.a, Val(:area_fraction)) .*
+                Interfacer.get_field(sims.a, Val(:emissivity)) .*
+                sims.a.cache_field .^ FT(4) .+
+                Interfacer.get_field(sims.b, Val(:area_fraction)) .*
+                Interfacer.get_field(sims.b, Val(:emissivity)) .*
+                sims.b.cache_field .^ FT(4)
+            ) .^ FT(1 / 4)
+        @test coupler_fields.T_sfc == expected_field_temp
+
         expected_field =
             Interfacer.get_field(sims.a, Val(:area_fraction)) .* sims.a.cache_field .+
             Interfacer.get_field(sims.b, Val(:area_fraction)) .* sims.b.cache_field
-        @test coupler_fields.T_sfc == expected_field
         @test coupler_fields.surface_direct_albedo == expected_field
         @test coupler_fields.surface_diffuse_albedo == expected_field
     end

test/field_exchanger_tests.jl

@@ -77,6 +77,7 @@ Interfacer.get_field(
     sim::TestOcean,
     ::Union{Val{:surface_direct_albedo}, Val{:surface_diffuse_albedo}},
 ) = sim.integrator.p.α
+Interfacer.get_field(sim::TestOcean, ::Val{:emissivity}) = eltype(sim.integrator.T)(1)
 
 FieldExchanger.import_atmos_fields!(csf, sim::TestOcean, atmos_sim) = nothing
 
@@ -96,6 +97,8 @@ Interfacer.get_field(sim::DummySurfaceSimulation3, ::Val{:surface_temperature})
 Interfacer.get_field(sim::DummySurfaceSimulation3, ::Val{:area_fraction}) =
     sim.integrator.p.area_fraction
 Interfacer.get_field(sim::DummySurfaceSimulation3, ::Val{:beta}) = sim.integrator.p.beta
+Interfacer.get_field(sim::DummySurfaceSimulation3, ::Val{:emissivity}) =
+    eltype(sim.integrator.T)(1)
 
 for FT in (Float32, Float64)
     @testset "calculate correct fluxes: dry for FT=$FT" begin