Skip to content

Commit abac61c

Browse files
juliasloan25juliasloan25
committed
filling in clima sea ice [skip ci]
1 parent c068af0 commit abac61c

File tree

3 files changed

+359
-139
lines changed

3 files changed

+359
-139
lines changed

‎experiments/ClimaEarth/components/ocean/clima_seaice.jl‎

Lines changed: 216 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,8 @@ import Thermodynamics as TD
77
import ClimaOcean.EN4: download_dataset
88
using KernelAbstractions: @kernel, @index, @inbounds
99

10+
include("climaocean_helpers.jl")
11+
1012
"""
1113
ClimaSeaIceSimulation{SIM, A, OPROP, REMAP}
1214
@@ -19,8 +21,11 @@ It contains the following objects:
1921
- `area_fraction::A`: A ClimaCore Field representing the surface area fraction of this component model on the exchange grid.
2022
- `melting_speed::MS`: An constant characteristic speed for melting/freezing.
2123
- `remapping::REMAP`: Objects needed to remap from the exchange (spectral) grid to Oceananigans spaces.
24+
TODO add skin_temperature field, for now set skin temperature to top layer sea ice temperature
25+
alt: compute skin temp from bulk temp and ocean surface temp, assuming linear profile in the top layer,
26+
might need to limit min sea ice thickness to avoid crazy gradients
2227
"""
23-
struct ClimaSeaIceSimulation{SIM,A,MS,REMAP} <: Interfacer.SeaIceModelSimulation
28+
struct ClimaSeaIceSimulation{SIM, A, MS, REMAP} <: Interfacer.SeaIceModelSimulation
2429
sea_ice::SIM
2530
area_fraction::A
2631
melting_speed::MS
@@ -30,7 +35,7 @@ end
3035
"""
3136
ClimaSeaIceSimulation()
3237
33-
Creates an OceananigansSimulation object containing a model, an integrator, and
38+
Creates an ClimaSeaIceSimulation object containing a model, an integrator, and
3439
a surface area fraction field.
3540
This type is used to indicate that this simulation is an ocean simulation for
3641
dispatch in coupling.
@@ -46,6 +51,8 @@ function ClimaSeaIceSimulation(area_fraction, ocean; output_dir)
4651
sea_ice = CO.sea_ice_simulation(grid, ocean.ocean; advection)
4752

4853
melting_speed = 1e-4
54+
55+
# Since ocean and sea ice share the same grid, we can also share the remapping objects
4956
remapping = ocean.remapping
5057

5158
# Before version 0.96.22, the NetCDFWriter was broken on GPU
@@ -68,3 +75,210 @@ function ClimaSeaIceSimulation(area_fraction, ocean; output_dir)
6875
sim = ClimaSeaIceSimulation(sea_ice, area_fraction, melting_speed, remapping)
6976
return sim
7077
end
78+
79+
###############################################################################
80+
### Functions required by ClimaCoupler.jl for a SurfaceModelSimulation
81+
###############################################################################
82+
83+
# Timestep the simulation forward to time `t`
84+
Interfacer.step!(sim::ClimaSeaIceSimulation, t) =
85+
OC.time_step!(sim.sea_ice, float(t) - sim.sea_ice.model.clock.time)
86+
87+
88+
Interfacer.get_field(sim::ClimaSeaIceSimulation, ::Val{:area_fraction}) = sim.area_fraction
89+
90+
# TODO: Better values for this
91+
# At the moment, we return always Float32. This is because we always want to run
92+
# Oceananingans with Float64, so we have no way to know the float type here. Sticking with
93+
# Float32 ensures that nothing is accidentally promoted to Float64. We will need to change
94+
# this anyway.
95+
Interfacer.get_field(sim::ClimaSeaIceSimulation, ::Val{:roughness_buoyancy}) =
96+
Float32(5.8e-5)
97+
Interfacer.get_field(sim::ClimaSeaIceSimulation, ::Val{:roughness_momentum}) =
98+
Float32(5.8e-5)
99+
Interfacer.get_field(sim::ClimaSeaIceSimulation, ::Val{:beta}) = Float32(1)
100+
Interfacer.get_field(sim::ClimaSeaIceSimulation, ::Val{:surface_direct_albedo}) =
101+
Float32(0.011)
102+
Interfacer.get_field(sim::ClimaSeaIceSimulation, ::Val{:surface_diffuse_albedo}) =
103+
Float32(0.069)
104+
105+
# TODO approximate surface temp from bulk temp and ocean surface temp, assuming linear profile in the top layer
106+
# TODO how to get ocean surface temp here? Is it in the sea ice or do we need to pass the ocean sim too?
107+
Interfacer.get_field(sim::ClimaSeaIceSimulation, ::Val{:surface_temperature}) =
108+
273.15 + sim.ocean.model.tracers.T
109+
110+
"""
111+
FluxCalculator.update_turbulent_fluxes!(sim::ClimaSeaIceSimulation, fields)
112+
113+
Update the turbulent fluxes in the simulation using the values stored in the coupler fields.
114+
These include latent heat flux, sensible heat flux, momentum fluxes, and moisture flux.
115+
116+
A note on sign conventions:
117+
SurfaceFluxes and ClimaSeaIce both use the convention that a positive flux is an upward flux.
118+
No sign change is needed during the exchange, except for moisture/salinity fluxes:
119+
SurfaceFluxes provides moisture moving from atmosphere to ocean as a negative flux at the surface,
120+
and ClimaSeaIce represents moisture moving from atmosphere to ocean as a positive salinity flux,
121+
so a sign change is needed when we convert from moisture to salinity flux.
122+
"""
123+
function FluxCalculator.update_turbulent_fluxes!(sim::ClimaSeaIceSimulation, fields)
124+
# Only LatitudeLongitudeGrid are supported because otherwise we have to rotate the vectors
125+
126+
(; F_lh, F_sh, F_turb_ρτxz, F_turb_ρτyz, F_turb_moisture) = fields
127+
grid = sim.sea_ice.model.grid
128+
129+
# Remap momentum fluxes onto reduced 2D Center, Center fields using scratch arrays and fields
130+
CC.Remapping.interpolate!(
131+
sim.remapping.scratch_arr1,
132+
sim.remapping.remapper_cc,
133+
F_turb_ρτxz,
134+
)
135+
OC.set!(sim.remapping.scratch_cc1, sim.remapping.scratch_arr1) # zonal momentum flux
136+
CC.Remapping.interpolate!(
137+
sim.remapping.scratch_arr2,
138+
sim.remapping.remapper_cc,
139+
F_turb_ρτyz,
140+
)
141+
OC.set!(sim.remapping.scratch_cc2, sim.remapping.scratch_arr2) # meridional momentum flux
142+
143+
# Rename for clarity; these are now Center, Center Oceananigans fields
144+
F_turb_ρτxz_cc = sim.remapping.scratch_cc1
145+
F_turb_ρτyz_cc = sim.remapping.scratch_cc2
146+
147+
# Set the momentum flux BCs at the correct locations using the remapped scratch fields
148+
oc_flux_u = surface_flux(sim.sea_ice.model.velocities.u)
149+
oc_flux_v = surface_flux(sim.sea_ice.model.velocities.v)
150+
set_from_extrinsic_vectors!(
151+
(; u = oc_flux_u, v = oc_flux_v),
152+
grid,
153+
F_turb_ρτxz_cc,
154+
F_turb_ρτyz_cc,
155+
)
156+
157+
# Remap the latent and sensible heat fluxes using scratch arrays
158+
CC.Remapping.interpolate!(sim.remapping.scratch_arr1, sim.remapping.remapper_cc, F_lh) # latent heat flux
159+
CC.Remapping.interpolate!(sim.remapping.scratch_arr2, sim.remapping.remapper_cc, F_sh) # sensible heat flux
160+
161+
# Rename for clarity; recall F_turb_energy = F_lh + F_sh
162+
remapped_F_lh = sim.remapping.scratch_arr1
163+
remapped_F_sh = sim.remapping.scratch_arr2
164+
165+
166+
# TODO what is sim.sea_ice.model.ice_thermodynamics.top_surface_temperature? where is it set?
167+
# TODO ocean_reference_density -> sea_ice.model.ice_density ?
168+
oc_flux_T = surface_flux(sim.ocean.model.tracers.T)
169+
OC.interior(oc_flux_T, :, :, 1) .=
170+
OC.interior(oc_flux_T, :, :, 1) .+
171+
(remapped_F_lh .+ remapped_F_sh) ./ (ocean_reference_density * ocean_heat_capacity)
172+
173+
# Add the part of the salinity flux that comes from the moisture flux. We also need to
174+
# add the component due to precipitation (that was done with the radiative fluxes)
175+
CC.Remapping.interpolate!(
176+
sim.remapping.scratch_arr1,
177+
sim.remapping.remapper_cc,
178+
F_turb_moisture,
179+
)
180+
moisture_fresh_water_flux = sim.remapping.scratch_arr1 ./ ocean_fresh_water_density # TODO do we need this for sea ice too?
181+
oc_flux_S = surface_flux(sim.sea_ice.model.tracers.S)
182+
surface_salinity = OC.interior(sim.sea_ice.model.tracers.S, :, :, 1)
183+
OC.interior(oc_flux_S, :, :, 1) .=
184+
OC.interior(oc_flux_S, :, :, 1) .- surface_salinity .* moisture_fresh_water_flux
185+
return nothing
186+
end
187+
188+
function Interfacer.update_field!(sim::OceananigansSimulation, ::Val{:area_fraction}, field)
189+
sim.area_fraction .= field
190+
return nothing
191+
end
192+
193+
"""
194+
FieldExchanger.update_sim!(sim::OceananigansSimulation, csf, area_fraction)
195+
196+
Update the ocean simulation with the provided fields, which have been filled in
197+
by the coupler.
198+
199+
Update the portion of the surface_fluxes for T and S that is due to radiation and
200+
precipitation. The rest will be updated in `update_turbulent_fluxes!`.
201+
202+
A note on sign conventions:
203+
ClimaAtmos and Oceananigans both use the convention that a positive flux is an upward flux.
204+
No sign change is needed during the exchange, except for precipitation/salinity fluxes.
205+
ClimaAtmos provides precipitation as a negative flux at the surface, and
206+
Oceananigans represents precipitation as a positive salinity flux,
207+
so a sign change is needed when we convert from precipitation to salinity flux.
208+
"""
209+
function FieldExchanger.update_sim!(sim::OceananigansSimulation, csf, area_fraction)
210+
(; ocean_reference_density, ocean_heat_capacity, ocean_fresh_water_density) =
211+
sim.ocean_properties
212+
213+
# Remap radiative flux onto scratch array; rename for clarity
214+
CC.Remapping.interpolate!(
215+
sim.remapping.scratch_arr1,
216+
sim.remapping.remapper_cc,
217+
csf.F_radiative,
218+
)
219+
remapped_F_radiative = sim.remapping.scratch_arr1
220+
221+
# Update only the part due to radiative fluxes. For the full update, the component due
222+
# to latent and sensible heat is missing and will be updated in update_turbulent_fluxes.
223+
oc_flux_T = surface_flux(sim.ocean.model.tracers.T)
224+
OC.interior(oc_flux_T, :, :, 1) .=
225+
remapped_F_radiative ./ (ocean_reference_density * ocean_heat_capacity)
226+
227+
# Remap precipitation fields onto scratch arrays; rename for clarity
228+
CC.Remapping.interpolate!(
229+
sim.remapping.scratch_arr1,
230+
sim.remapping.remapper_cc,
231+
csf.P_liq,
232+
)
233+
CC.Remapping.interpolate!(
234+
sim.remapping.scratch_arr2,
235+
sim.remapping.remapper_cc,
236+
csf.P_snow,
237+
)
238+
remapped_P_liq = sim.remapping.scratch_arr1
239+
remapped_P_snow = sim.remapping.scratch_arr2
240+
241+
# Virtual salt flux
242+
oc_flux_S = surface_flux(sim.ocean.model.tracers.S)
243+
precipitating_fresh_water_flux =
244+
(remapped_P_liq .+ remapped_P_snow) ./ ocean_fresh_water_density
245+
surface_salinity_flux =
246+
OC.interior(sim.ocean.model.tracers.S, :, :, 1) .* precipitating_fresh_water_flux
247+
OC.interior(oc_flux_S, :, :, 1) .= .-surface_salinity_flux
248+
return nothing
249+
end
250+
251+
# TODO add stub to FluxCalculator
252+
# TODO fill this in using ClimaOcean: https://github.com/CliMA/ClimaOcean.jl/pull/627
253+
# TODO instead of cs we can take in the ocean and sea ice sims, dispatch on them to do nothing for other sea ice models
254+
# TODO docstring
255+
function ocean_seaice_fluxes!(cs)
256+
seaice_sim = cs.models.sea_ice
257+
ocean_sim = cs.models.ocean
258+
259+
melting_speed = seaice_sim.melting_speed
260+
ocean_properties = ocean_sim.ocean_properties
261+
262+
# TODO what should fluxes be here?
263+
OC.compute_sea_ice_ocean_fluxes!(
264+
fluxes,
265+
ocean_sim.ocean,
266+
seaice_sim.sea_ice,
267+
melting_speed,
268+
ocean_properties,
269+
)
270+
return nothing
271+
end
272+
273+
"""
274+
get_model_prog_state(sim::ClimaSeaIceSimulation)
275+
276+
Returns the model state of a simulation as a `ClimaCore.FieldVector`.
277+
It's okay to leave this unimplemented for now, but we won't be able to use the
278+
restart system.
279+
280+
TODO extend this for non-ClimaCore states.
281+
"""
282+
function Checkpointer.get_model_prog_state(sim::ClimaSeaIceSimulation)
283+
@warn "get_model_prog_state not implemented for ClimaSeaIceSimulation"
284+
end

‎experiments/ClimaEarth/components/ocean/climaocean_helpers.jl‎

Lines changed: 128 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,128 @@
1+
2+
# We always want the surface, so we always set zero(pt.lat) for z
3+
"""
4+
to_node(pt::CC.Geometry.LatLongPoint)
5+
6+
Transform `LatLongPoint` into a tuple (long, lat, 0), where the 0 is needed because we only
7+
care about the surface.
8+
"""
9+
@inline to_node(pt::CC.Geometry.LatLongPoint) = pt.long, pt.lat, zero(pt.lat)
10+
# This next one is needed if we have "LevelGrid"
11+
@inline to_node(pt::CC.Geometry.LatLongZPoint) = pt.long, pt.lat, zero(pt.lat)
12+
13+
"""
14+
map_interpolate(points, oc_field::OC.Field)
15+
16+
Interpolate the given 3D field onto the target points.
17+
18+
If the underlying grid does not contain a given point, return 0 instead.
19+
20+
TODO: Use a non-allocating version of this function (simply replace `map` with `map!`)
21+
"""
22+
function map_interpolate(points, oc_field::OC.Field)
23+
loc = map(L -> L(), OC.Fields.location(oc_field))
24+
grid = oc_field.grid
25+
data = oc_field.data
26+
27+
# TODO: There has to be a better way
28+
min_lat, max_lat = extrema(OC.φnodes(grid, OC.Center(), OC.Center(), OC.Center()))
29+
30+
map(points) do pt
31+
FT = eltype(pt)
32+
33+
# The oceananigans grid does not cover the entire globe, so we should not
34+
# interpolate outside of its latitude bounds. Instead we return 0
35+
min_lat < pt.lat < max_lat || return FT(0)
36+
37+
fᵢ = OC.Fields.interpolate(to_node(pt), data, loc, grid)
38+
convert(FT, fᵢ)::FT
39+
end
40+
end
41+
42+
"""
43+
surface_flux(f::OC.AbstractField)
44+
45+
Extract the top boundary conditions for the given field.
46+
"""
47+
function surface_flux(f::OC.AbstractField)
48+
top_bc = f.boundary_conditions.top
49+
if top_bc isa OC.BoundaryCondition{<:OC.BoundaryConditions.Flux}
50+
return top_bc.condition
51+
else
52+
return nothing
53+
end
54+
end
55+
56+
function Interfacer.remap(field::OC.Field, target_space)
57+
return map_interpolate(CC.Fields.coordinate_field(target_space), field)
58+
end
59+
60+
function Interfacer.remap(operation::OC.AbstractOperations.AbstractOperation, target_space)
61+
evaluated_field = OC.Field(operation)
62+
OC.compute!(evaluated_field)
63+
return Interfacer.remap(evaluated_field, target_space)
64+
end
65+
66+
"""
67+
set_from_extrinsic_vectors!(vectors, grid, u_cc, v_cc)
68+
69+
Given the extrinsic vector components `u_cc` and `v_cc` as `Center, Center`
70+
fields, rotate them onto the target grid and remap to `Face, Center` and
71+
`Center, Face` fields, respectively.
72+
"""
73+
function set_from_extrinsic_vectors!(vectors, grid, u_cc, v_cc)
74+
arch = grid.architecture
75+
76+
# Rotate vectors onto the grid
77+
OC.Utils.launch!(arch, grid, :xy, _rotate_velocities!, u_cc, v_cc, grid)
78+
79+
# Fill halo regions with the rotated vectors so we can use them to interpolate
80+
OC.fill_halo_regions!(u_cc)
81+
OC.fill_halo_regions!(v_cc)
82+
83+
# Interpolate the vectors to face/center and center/face respectively
84+
OC.Utils.launch!(
85+
arch,
86+
grid,
87+
:xy,
88+
_interpolate_velocities!,
89+
vectors.u,
90+
vectors.v,
91+
grid,
92+
u_cc,
93+
v_cc,
94+
)
95+
return nothing
96+
end
97+
98+
"""
99+
_rotate_velocities!(u, v, grid)
100+
101+
Rotate the velocities from the extrinsic coordinate system to the intrinsic
102+
coordinate system.
103+
"""
104+
@kernel function _rotate_velocities!(u, v, grid)
105+
# Use `k = 1` to index into the reduced Fields
106+
i, j = @index(Global, NTuple)
107+
# Rotate u, v from extrinsic to intrinsic coordinate system
108+
ur, vr = OC.Operators.intrinsic_vector(i, j, 1, grid, u, v)
109+
@inbounds begin
110+
u[i, j, 1] = ur
111+
v[i, j, 1] = vr
112+
end
113+
end
114+
115+
"""
116+
_interpolate_velocities!(u, v, grid, u_cc, v_cc)
117+
118+
Interpolate the input velocities `u_cc` and `v_cc`, which are Center/Center
119+
Fields to Face/Center and Center/Face coordinates, respectively.
120+
"""
121+
@kernel function _interpolate_velocities!(u, v, grid, u_cc, v_cc)
122+
# Use `k = 1` to index into the reduced Fields
123+
i, j = @index(Global, NTuple)
124+
@inbounds begin
125+
u[i, j, 1] = OC.Operators.ℑxyᶠᶜᵃ(i, j, 1, grid, u_cc)
126+
v[i, j, 1] = OC.Operators.ℑxyᶜᶠᵃ(i, j, 1, grid, v_cc)
127+
end
128+
end

0 commit comments

Comments
 (0)