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