Skip to content

Commit 7e4eb70

Browse files
juliasloan25juliasloan25
committed
use individual model surface temps for surface fluxes
1 parent b5ef5e8 commit 7e4eb70

File tree

13 files changed

+348
-236
lines changed

13 files changed

+348
-236
lines changed

‎NEWS.md‎

Lines changed: 18 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,24 @@ ClimaCoupler.jl Release Notes
66

77
### ClimaCoupler features
88

9+
#### Use individual surface model temperatures to compute turbulent fluxes PR[#1498](https://github.com/CliMA/ClimaCoupler.jl/pull/1498)
10+
We use a partitioned approach to computing turbulent fluxes, so we should be using
11+
the individual surface temperature (and humidity, air density, etc) for each component
12+
model. Previously we were using a combined surface temperature across all surface models.
13+
This change doesn't affect simulations using integer area fractions, but will affect simulations
14+
with fractional area fractions.
15+
16+
This PR also adds ensures that the area fraction of an OceananigansSimulation
17+
never exceeds the latitude limits when using a LatitudeLongitudeGrid. Note that
18+
this simulation requires running with an ice model, which is used to fill in
19+
the polar regions.
20+
21+
#### Combine LW fluxes to get surface temperature for radiation PR[#1492](https://github.com/CliMA/ClimaCoupler.jl/issues/1492)
22+
Previously, we combined the temperatures of each surface models directly using an
23+
area-weighted sum. Now, we instead compute the longwave flux for each component, compute
24+
the area fraction-weighted sum of those, and then convert the total flux back to get
25+
surface temperature. This temperature is then provided to the atmosphere for radiation.
26+
927
#### Remove coarse nightly CMIP PR[#1485](https://github.com/CliMA/ClimaCoupler.jl/pull/1485)
1028
To avoid depending on the main branch of too many packages in the nightly pipeline,
1129
we remove the CMIP nightly run and will only test AMIP nightly.

‎docs/src/interfacer.md‎

Lines changed: 48 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -77,31 +77,47 @@ exchange the extra fields only when necessary. All coupler fields are defined on
7777
boundary space.
7878
- Any additional fields specified here will likely also require an `update_field!`
7979
method defined for this component model, so the coupler can update the component.
80+
They may also require a new method of `import_atmos_fields!` or `combine_surfaces!`
81+
to update the coupler fields from the component model that computes the field.
8082

8183
The default coupler exchange fields are the following, defined in
8284
`default_coupler_fields()` in the Interfacer module:
8385

84-
| Coupler name | Description | Units |
85-
|-------------------|-----------------------------------------------------|------------|
86-
| `z0m_sfc` | momentum roughness length | m |
87-
| `z0b_sfc` | buoyancy roughness length | m |
88-
| `beta` | factor to scale evaporation from the surface | - |
89-
| `emissivity` | surface emissivity | - |
90-
| `T_atmos` | atmosphere temperature at the bottom layer | K |
91-
| `q_atmos` | atmosphere humidity at the bottom layer | kg kg⁻¹ |
92-
| `ρ_atmos` | atmosphere air density at the bottom layer | kg m⁻³ |
93-
| `T_sfc` | surface temperature, averaged across components | K |
94-
| `q_sfc` | surface humidity | kg kg⁻¹ |
95-
| `F_lh` | latent heat flux | W m⁻² |
96-
| `F_sh` | sensible heat flux | W m⁻² |
97-
| `F_turb_moisture` | turbulent moisture flux | kg m⁻² s⁻¹ |
98-
| `F_turb_ρτxz` | turbulent momentum flux in the zonal direction | kg m⁻¹ s⁻² |
99-
| `F_turb_ρτyz` | turbulent momentum flux in the meridional direction | kg m⁻¹ s⁻² |
100-
| `F_radiative` | net radiative flux at the surface | W m⁻² |
101-
| `P_liq` | liquid precipitation | kg m⁻² s⁻¹ |
102-
| `P_snow` | snow precipitation | kg m⁻² s⁻¹ |
103-
| `temp1` | a surface field used for intermediate calculations | - |
104-
| `temp2` | a surface field used for intermediate calculations | - |
86+
| Coupler name | Description | Units |
87+
|-------------------|-------------------------------------------------------------|------------|
88+
| `T_atmos` | atmosphere temperature at the bottom layer | K |
89+
| `q_atmos` | atmosphere humidity at the bottom layer | kg kg⁻¹ |
90+
| `ρ_atmos` | atmosphere air density at the bottom layer | kg m⁻³ |
91+
| `z_int` | height of the first internal atmosphere level (center) | m |
92+
| `z_sfc` | height of the bottom atmosphere layer (face) | m |
93+
| `F_lh` | latent heat flux | W m⁻² |
94+
| `F_sh` | sensible heat flux | W m⁻² |
95+
| `F_turb_moisture` | turbulent moisture flux | kg m⁻² s⁻¹ |
96+
| `F_turb_ρτxz` | turbulent momentum flux in the zonal direction | kg m⁻¹ s⁻² |
97+
| `F_turb_ρτyz` | turbulent momentum flux in the meridional direction | kg m⁻¹ s⁻² |
98+
| `F_radiative` | net radiative flux at the surface | W m⁻² |
99+
| `emissivity` | surface emissivity | - |
100+
| `T_sfc` | surface temperature, averaged across components | K |
101+
| `P_liq` | liquid precipitation | kg m⁻² s⁻¹ |
102+
| `P_snow` | snow precipitation | kg m⁻² s⁻¹ |
103+
| `scalar_temp1` | a surface scalar field used for intermediate calculations | - |
104+
| `scalar_temp2` | a surface scalar field used for intermediate calculations | - |
105+
| `scalar_temp3` | a surface scalar field used for intermediate calculations | - |
106+
| `scalar_temp4` | a surface scalar field used for intermediate calculations | - |
107+
108+
!!! note "What should be stored in the coupler exchange fields?"
109+
In general, the coupler fields should contain exchange fields for fluxes, including
110+
turbulent fluxes, radiative fluxes, and precipitation. They also hold any quantities
111+
that a component model requires from another component. For example, the atmosphere needs
112+
surface temperature and emissivity from the surface models to compute radiation, so the
113+
coupler allocates space to exchange them.
114+
The coupler exchange fields may also hold quantities from components that are used to
115+
compute turbulent fluxes. As a general rule, we tend to store such quantities that come
116+
from the atmosphere model, but access them when needed for surface models. This is because we
117+
compute fluxes indvidually for the interface between each surface and the atmosphere
118+
model. As a result, the atmosphere quantities are used for each of these calculations,
119+
so storing them in the coupler fields allows us to avoid regridding them to the coupler
120+
space multiple times per coupling timestep.
105121

106122
- `update_sim!(::ComponentModelSimulation, csf)`: A
107123
function to update each of the fields of the component model simulation
@@ -186,9 +202,10 @@ for the following properties:
186202
| `LW_d` | downwards longwave flux | W m⁻² |
187203
| `SW_d` | downwards shortwave flux | W m⁻² |
188204

189-
Note that `co2`, `diffuse_fraction`, `LW_d` and
190-
`SW_d` will not be present in a `ClimaAtmosSimulation` if the model is setup with no radiation.
191-
Because of this, a `ClimaAtmosSimulation` must have radiation if running with the full `ClimaLand` model.
205+
!!! note
206+
`co2`, `diffuse_fraction`, `LW_d` and `SW_d` will not be present in a `ClimaAtmosSimulation`
207+
if the model is setup with no radiation. Because of this, a `ClimaAtmosSimulation` must have
208+
radiation enabled if running with the full `ClimaLand` model.
192209

193210
### SurfaceModelSimulation - required functions
194211
Analogously to the `AtmosModelSimulation`, a `SurfaceModelSimulation`
@@ -207,8 +224,9 @@ for the following properties:
207224
| `surface_diffuse albedo` | bulk diffuse surface albedo | |
208225
| `surface_temperature` | surface temperature | K |
209226

210-
Note: `area_fraction` is expected to be defined on the boundary space of the simulation,
211-
while all other fields will likely be on the simulation's own space.
227+
!!! note
228+
`area_fraction` is expected to be defined on the boundary space of the simulation,
229+
while all other fields will likely be on the simulation's own space.
212230

213231
- `update_field!(::SurfaceModelSimulation, ::Val{property}, field)`:
214232
A function to update the value of property in the component model
@@ -230,9 +248,10 @@ properties needed by a component model.
230248
| `turbulent_energy_flux` | aerodynamic turbulent surface fluxes of energy (sensible and latent heat) | W m⁻² |
231249
| `turbulent_moisture_flux` | aerodynamic turbulent surface fluxes of energy (evaporation) | kg m⁻² s⁻¹ |
232250

233-
Note: `update_field!(::SurfaceModelSimulation, ::Val{:area_fraction}, field)` is
234-
not required to be extended for land models, since they're assumed to have a
235-
constant area fraction.
251+
!!! note
252+
`update_field!(::SurfaceModelSimulation, ::Val{:area_fraction}, field)` is
253+
not required to be extended for land models, since they're assumed to have a
254+
constant area fraction.
236255

237256
### SurfaceModelSimulation - optional functions
238257
- `get_field(::SurfaceModelSimulation, ::Val{property})`:

‎experiments/ClimaEarth/components/atmosphere/climaatmos.jl‎

Lines changed: 28 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -350,28 +350,49 @@ function Interfacer.update_field!(
350350
temp_field_surface = sim.integrator.p.scratch.ᶠtemp_field_level
351351
@assert axes(temp_field_surface) == atmos_surface_space
352352

353+
# Compute surface humidity on the coupler space, then remap to atmosphere surface space
354+
Interfacer.get_field!(csf.scalar_temp1, sim, Val(:air_temperature))
355+
T_atmos = csf.scalar_temp1
356+
Interfacer.get_field!(csf.scalar_temp2, sim, Val(:specific_humidity))
357+
q_atmos = csf.scalar_temp2
358+
Interfacer.get_field!(csf.scalar_temp3, sim, Val(:air_density))
359+
ρ_atmos = csf.scalar_temp3
360+
361+
thermo_params = get_thermo_params(sim)
362+
FluxCalculator.compute_surface_humidity!(
363+
csf.scalar_temp4,
364+
T_atmos,
365+
q_atmos,
366+
ρ_atmos,
367+
csf.T_sfc,
368+
thermo_params,
369+
)
370+
371+
# Remap surface temperature and humidity to atmosphere surface space
353372
# NOTE: This is allocating! If we had 2 more scratch fields, we could avoid this
354373
T_sfc_atmos = Interfacer.remap(csf.T_sfc, atmos_surface_space)
355-
q_sfc_atmos = Interfacer.remap(csf.q_sfc, atmos_surface_space)
374+
q_sfc_atmos = Interfacer.remap(csf.scalar_temp4, atmos_surface_space)
375+
356376
# Store `ρ_sfc_atmos` in an atmosphere scratch field on the surface space
357377
temp_field_surface =
358378
FluxCalculator.extrapolate_ρ_to_sfc.(
359-
get_thermo_params(sim),
379+
thermo_params,
360380
sim.integrator.p.precomputed.sfc_conditions.ts,
361381
T_sfc_atmos,
362-
) # ρ_sfc_atmos
382+
)
383+
ρ_sfc_atmos = temp_field_surface
363384

364385
if sim.integrator.p.atmos.moisture_model isa CA.DryModel
365386
sim.integrator.p.precomputed.sfc_conditions.ts .=
366-
TD.PhaseDry_ρT.(get_thermo_params(sim), temp_field_surface, T_sfc_atmos) # temp_field_surface = ρ_sfc_atmos
387+
TD.PhaseDry_ρT.(thermo_params, ρ_sfc_atmos, T_sfc_atmos)
367388
else
368389
sim.integrator.p.precomputed.sfc_conditions.ts .=
369390
TD.PhaseNonEquil_ρTq.(
370-
get_thermo_params(sim),
371-
temp_field_surface,
391+
thermo_params,
392+
ρ_sfc_atmos,
372393
T_sfc_atmos,
373394
TD.PhasePartition.(q_sfc_atmos),
374-
) # temp_field_surface = ρ_sfc_atmos
395+
)
375396
end
376397
end
377398
Interfacer.get_field(sim::ClimaAtmosSimulation, ::Val{:height_int}) =

‎experiments/ClimaEarth/components/land/climaland_integrated.jl‎

Lines changed: 43 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -447,16 +447,21 @@ function FieldExchanger.update_sim!(sim::ClimaLandSimulation, csf, area_fraction
447447
Interfacer.update_field!(sim, Val(:snow_precipitation), csf.P_snow)
448448
end
449449

450+
"""
451+
FieldExchanger.import_atmos_fields!(csf, sim::ClimaLandSimulation, atmos_sim)
452+
453+
Import non-default coupler fields from the atmosphere simulation into the coupler fields.
454+
These include the diffuse fraction of light, shortwave and longwave downwelling radiation,
455+
air pressure, and CO2 concentration.
456+
457+
The default coupler fields will be imported by the default method implemented in
458+
FieldExchanger.jl.
459+
"""
450460
function FieldExchanger.import_atmos_fields!(csf, sim::ClimaLandSimulation, atmos_sim)
451461
Interfacer.get_field!(csf.diffuse_fraction, atmos_sim, Val(:diffuse_fraction))
452462
Interfacer.get_field!(csf.SW_d, atmos_sim, Val(:SW_d))
453463
Interfacer.get_field!(csf.LW_d, atmos_sim, Val(:LW_d))
454464
Interfacer.get_field!(csf.P_atmos, atmos_sim, Val(:air_pressure))
455-
Interfacer.get_field!(csf.T_atmos, atmos_sim, Val(:air_temperature))
456-
Interfacer.get_field!(csf.q_atmos, atmos_sim, Val(:specific_humidity))
457-
Interfacer.get_field!(csf.ρ_atmos, atmos_sim, Val(:air_density))
458-
Interfacer.get_field!(csf.P_liq, atmos_sim, Val(:liquid_precipitation))
459-
Interfacer.get_field!(csf.P_snow, atmos_sim, Val(:snow_precipitation))
460465
# CO2 is a scalar for now so it doesn't need remapping
461466
csf.c_co2 .= Interfacer.get_field(atmos_sim, Val(:co2))
462467
return nothing
@@ -596,75 +601,82 @@ function FluxCalculator.compute_surface_fluxes!(
596601
# Combine turbulent energy fluxes from each component of the land model
597602
# Use temporary variables to avoid allocating
598603
Interfacer.remap!(
599-
csf.temp1,
604+
csf.scalar_temp1,
600605
canopy_dest.lhf .+ soil_dest.lhf .* (1 .- p.snow.snow_cover_fraction) .+
601606
p.snow.snow_cover_fraction .* snow_dest.lhf,
602607
)
603608
Interfacer.remap!(
604-
csf.temp2,
609+
csf.scalar_temp2,
605610
canopy_dest.shf .+ soil_dest.shf .* (1 .- p.snow.snow_cover_fraction) .+
606611
p.snow.snow_cover_fraction .* snow_dest.shf,
607612
)
608613

609614
# Zero out the fluxes where the area fraction is zero
610-
@. csf.temp1 = ifelse(area_fraction == 0, zero(csf.temp1), csf.temp1)
611-
@. csf.temp2 = ifelse(area_fraction == 0, zero(csf.temp2), csf.temp2)
615+
@. csf.scalar_temp1 =
616+
ifelse(area_fraction == 0, zero(csf.scalar_temp1), csf.scalar_temp1)
617+
@. csf.scalar_temp2 =
618+
ifelse(area_fraction == 0, zero(csf.scalar_temp2), csf.scalar_temp2)
612619

613620
# Update the coupler field in-place
614-
@. csf.F_lh += csf.temp1 * area_fraction
615-
@. csf.F_sh += csf.temp2 * area_fraction
621+
@. csf.F_lh += csf.scalar_temp1 * area_fraction
622+
@. csf.F_sh += csf.scalar_temp2 * area_fraction
616623

617624
# Combine turbulent moisture fluxes from each component of the land model
618625
# Note that we multiply by ρ_liq to convert from m s-1 to kg m-2 s-1
619626
ρ_liq = (LP.ρ_cloud_liq(sim.model.soil.parameters.earth_param_set))
620627
Interfacer.remap!(
621-
csf.temp1,
628+
csf.scalar_temp1,
622629
(
623630
canopy_dest.transpiration .+
624631
(soil_dest.vapor_flux_liq .+ soil_dest.vapor_flux_ice) .*
625632
(1 .- p.snow.snow_cover_fraction) .+
626633
p.snow.snow_cover_fraction .* snow_dest.vapor_flux
627634
) .* ρ_liq,
628635
)
629-
@. csf.temp1 = ifelse(area_fraction == 0, zero(csf.temp1), csf.temp1)
630-
@. csf.F_turb_moisture += csf.temp1 * area_fraction
636+
@. csf.scalar_temp1 =
637+
ifelse(area_fraction == 0, zero(csf.scalar_temp1), csf.scalar_temp1)
638+
@. csf.F_turb_moisture += csf.scalar_temp1 * area_fraction
631639

632640
# Combine turbulent momentum fluxes from each component of the land model
633641
# Note that we exclude the canopy component here for now, since we can have nonzero momentum fluxes
634642
# where there is zero LAI. This should be fixed in ClimaLand.
635643
Interfacer.remap!(
636-
csf.temp1,
644+
csf.scalar_temp1,
637645
soil_dest.ρτxz .* (1 .- p.snow.snow_cover_fraction) .+
638646
p.snow.snow_cover_fraction .* snow_dest.ρτxz,
639647
)
640-
@. csf.temp1 = ifelse(area_fraction == 0, zero(csf.temp1), csf.temp1)
641-
@. csf.F_turb_ρτxz += csf.temp1 * area_fraction
648+
@. csf.scalar_temp1 =
649+
ifelse(area_fraction == 0, zero(csf.scalar_temp1), csf.scalar_temp1)
650+
@. csf.F_turb_ρτxz += csf.scalar_temp1 * area_fraction
642651

643652
Interfacer.remap!(
644-
csf.temp1,
653+
csf.scalar_temp1,
645654
soil_dest.ρτyz .* (1 .- p.snow.snow_cover_fraction) .+
646655
p.snow.snow_cover_fraction .* snow_dest.ρτyz,
647656
)
648-
@. csf.temp1 = ifelse(area_fraction == 0, zero(csf.temp1), csf.temp1)
649-
@. csf.F_turb_ρτyz += csf.temp1 * area_fraction
657+
@. csf.scalar_temp1 =
658+
ifelse(area_fraction == 0, zero(csf.scalar_temp1), csf.scalar_temp1)
659+
@. csf.F_turb_ρτyz += csf.scalar_temp1 * area_fraction
650660

651661
# Combine the buoyancy flux from each component of the land model
652662
# Note that we exclude the canopy component here for now, since ClimaLand doesn't
653663
# include its extra resistance term in the buoyancy flux calculation.
654664
Interfacer.remap!(
655-
csf.temp1,
665+
csf.scalar_temp1,
656666
soil_dest.buoy_flux .* (1 .- p.snow.snow_cover_fraction) .+
657667
p.snow.snow_cover_fraction .* snow_dest.buoy_flux,
658668
)
659-
@. csf.temp1 = ifelse(area_fraction == 0, zero(csf.temp1), csf.temp1)
660-
@. csf.buoyancy_flux += csf.temp1 * area_fraction
669+
@. csf.scalar_temp1 =
670+
ifelse(area_fraction == 0, zero(csf.scalar_temp1), csf.scalar_temp1)
671+
@. csf.buoyancy_flux += csf.scalar_temp1 * area_fraction
661672

662673
# Compute ustar from the momentum fluxes and surface air density
663674
# ustar = sqrt(ρτ / ρ)
664-
@. csf.temp1 = sqrt(sqrt(csf.F_turb_ρτxz^2 + csf.F_turb_ρτyz^2) / csf.ρ_atmos)
665-
@. csf.temp1 = ifelse(area_fraction == 0, zero(csf.temp1), csf.temp1)
675+
@. csf.scalar_temp1 = sqrt(sqrt(csf.F_turb_ρτxz^2 + csf.F_turb_ρτyz^2) / csf.ρ_atmos)
676+
@. csf.scalar_temp1 =
677+
ifelse(area_fraction == 0, zero(csf.scalar_temp1), csf.scalar_temp1)
666678
# If ustar is zero, set it to eps to avoid division by zero in the atmosphere
667-
@. csf.ustar += max(csf.temp1 * area_fraction, eps(FT))
679+
@. csf.ustar += max(csf.scalar_temp1 * area_fraction, eps(FT))
668680

669681
# Compute the Monin-Obukhov length from ustar and the buoyancy flux
670682
# L_MO = -u^3 / (k * buoyancy_flux)
@@ -674,11 +686,13 @@ function FluxCalculator.compute_surface_fluxes!(
674686
return abs(v) < eps(FT) ? eps(FT) * sign_of_v : v
675687
end
676688
surface_params = LP.surface_fluxes_parameters(sim.model.soil.parameters.earth_param_set)
677-
@. csf.temp1 =
689+
@. csf.scalar_temp1 =
678690
-csf.ustar^3 / SFP.von_karman_const(surface_params) / non_zero(csf.buoyancy_flux)
679-
@. csf.temp1 = ifelse(area_fraction == 0, zero(csf.temp1), csf.temp1)
691+
@. csf.scalar_temp1 =
692+
ifelse(area_fraction == 0, zero(csf.scalar_temp1), csf.scalar_temp1)
680693
# When L_MO is infinite, avoid multiplication by zero to prevent NaN
681-
@. csf.L_MO += ifelse(isinf(csf.temp1), csf.temp1, csf.temp1 * area_fraction)
694+
@. csf.L_MO +=
695+
ifelse(isinf(csf.scalar_temp1), csf.scalar_temp1, csf.scalar_temp1 * area_fraction)
682696

683697
return nothing
684698
end

0 commit comments

Comments
 (0)