Atmospheric-Dynamics-GUF / PinCFlow.jl / 18600325485

Committed 17 Oct 2025 05:35PM UTC coverage: 64.858% (-1.1%) from 65.961%

Build # 18600325485

Build Type

github

Committed by

Commit Message

Merge 6a432eb7e into 2192c9fbf

Pull Request Pull Request #117: Removed density reconstructions and fluxes in Boussinesq mode and restructured the BC functions.

Run Details

60 of 152 new or added lines in 13 files covered. (39.47%)

41 existing lines in 6 files now uncovered.

4889 of 7538 relevant lines covered (64.86%)

1417695.35 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

63.91

/src/FluxCalculator/compute_fluxes!.jl

1	"""
2	```julia
3	compute_fluxes!(state::State, predictands::Predictands)
4	```
5
6	Compute fluxes by dispatching to specialized methods for each prognostic variable.
7
8	```julia
9	compute_fluxes!(state::State, predictands::Predictands, variable::Rho)
10	```
11
12	Compute the density fluxes in all three directions, by dispatching to a model-specific method.
13
14	```julia
15	compute_fluxes!(
16	state::State,
17	predictands::Predictands,
18	variable::Rho,
19	model::Boussinesq,
20	)
21	```
22
23	Return in Boussinesq mode.
24
25	```julia
26	compute_fluxes!(
27	state::State,
28	predictands::Predictands,
29	variable::Rho,
30	model::Union{PseudoIncompressible, Compressible},
31	)
32	```
33
34	Compute the density fluxes in all three directions.
35
36	The fluxes are given by
37
38	```math
39	\\begin{align*}
40	\\mathcal{F}^{\\rho, \\widehat{x}}_{i + 1 / 2} & = \\frac{\\tau_{\\widehat{x}}}{2} \\left\\{\\left[1 + \\mathrm{sgn} \\left(\\tau_{\\widehat{x}}\\right)\\right] {\\widetilde{\\phi}}^\\mathrm{R} + \\left[1 - \\mathrm{sgn} \\left(\\tau_{\\widehat{x}}\\right)\\right] {\\widetilde{\\phi}}_{i + 1}^\\mathrm{L}\\right\\},\\\\
41	\\mathcal{F}^{\\rho, \\widehat{y}}_{j + 1 / 2} & = \\frac{\\tau_{\\widehat{y}}}{2} \\left\\{\\left[1 + \\mathrm{sgn} \\left(\\tau_{\\widehat{y}}\\right)\\right] {\\widetilde{\\phi}}^\\mathrm{F} + \\left[1 - \\mathrm{sgn} \\left(\\tau_{\\widehat{y}}\\right)\\right] {\\widetilde{\\phi}}_{j + 1}^\\mathrm{B}\\right\\},\\\\
42	\\mathcal{F}^{\\rho, \\widehat{z}}_{k + 1 / 2} & = \\frac{\\tau_{\\widehat{z}}}{2} \\left\\{\\left[1 + \\mathrm{sgn} \\left(\\tau_{\\widehat{z}}\\right)\\right] {\\widetilde{\\phi}}^\\mathrm{U} + \\left[1 - \\mathrm{sgn} \\left(\\tau_{\\widehat{z}}\\right)\\right] {\\widetilde{\\phi}}_{k + 1}^\\mathrm{D}\\right\\},
43	\\end{align*}
44	```
45
46	where
47
48
49	```math
50	\\begin{align*}
51	\\tau_{\\widehat{x}} & = \\left(J P_\\mathrm{old}\\right)_{i + 1 / 2} u_{\\mathrm{old}, i + 1 / 2},\\\\
52	\\tau_{\\widehat{y}} & = \\left(J P_\\mathrm{old}\\right)_{j + 1 / 2} v_{\\mathrm{old}, j + 1 / 2},\\\\
53	\\tau_{\\widehat{z}} & = \\left(J P_\\mathrm{old}\\right)_{k + 1 / 2} \\widehat{w}_{\\mathrm{old}, k + 1 / 2}
54	\\end{align*}
55	```
56
57	are the transporting velocities (weighted by the Jacobian) and ``\\widetilde{\\phi}`` is the reconstruction of ``\\rho / P_\\mathrm{old}``. More specifically, the superscripts ``\\mathrm{R}``, ``\\mathrm{L}``, ``\\mathrm{F}``, ``\\mathrm{B}``, ``\\mathrm{U}`` and ``\\mathrm{D}`` indicate reconstructions at the right, left, forward, backward, upward and downward cell interfaces of the respective grid points, respectively. Quantities with the subscript ``\\mathrm{old}`` are obtained from a previous state, which is partially passed to the method via `predictands`.
58
59
60	```julia
61	compute_fluxes!(state::State, predictands::Predictands, variable::RhoP)
62	```
63
64	Compute the density-fluctuations fluxes in all three directions.
65
66	The computation is analogous to that of the density fluxes.
67
68	```julia
69	compute_fluxes!(
70	state::State,
71	predictands::Predictands,
72	model::Union{Boussinesq, PseudoIncompressible},
73	variable::P,
74	)
75	```
76
77	Return in non-compressible modes.
78
79	```julia
80	compute_fluxes!(
81	state::State,
82	predictands::Predictands,
83	model::Compressible,
84	variable::P,
85	)
86	```
87
88	Compute the mass-weighted potential-temperature fluxes in all three directions.
89
90	The fluxes are given by
91
92	```math
93	\\begin{align*}
94	\\mathcal{F}^{P, \\widehat{x}}_{i + 1 / 2} & = \\left(J P_\\mathrm{old}\\right)_{i + 1 / 2} u_{\\mathrm{old}, i + 1 / 2},\\\\
95	\\mathcal{F}^{P, \\widehat{y}}_{j + 1 / 2} & = \\left(J P_\\mathrm{old}\\right)_{j + 1 / 2} v_{\\mathrm{old}, j + 1 / 2},\\\\
96	\\mathcal{F}^{P, \\widehat{z}}_{k + 1 / 2} & = \\left(J P_\\mathrm{old}\\right)_{k + 1 / 2} \\widehat{w}_{\\mathrm{old}, k + 1 / 2}.
97	\\end{align*}
98	```
99
100	```julia
101	compute_fluxes!(state::State, old_predictands::Predictands, variable::U)
102	```
103
104	Compute the zonal-momentum fluxes in all three directions.
105
106	The fluxes are first set to the advective parts
107
108	```math
109	\\begin{align*}
110	\\mathcal{F}^{\\rho u, \\widehat{x}}_{i + 1} & = \\frac{\\tau_{\\widehat{x}}}{2} \\left\\{\\left[1 + \\mathrm{sgn} \\left(\\tau_{\\widehat{x}}\\right)\\right] {\\widetilde{\\phi}}_{i + 1 / 2}^\\mathrm{R} + \\left[1 - \\mathrm{sgn} \\left(\\tau_{\\widehat{x}}\\right)\\right] {\\widetilde{\\phi}}_{i + 3 / 2}^\\mathrm{L}\\right\\},\\\\
111	\\mathcal{F}^{\\rho u, \\widehat{y}}_{i + 1 / 2, j + 1 / 2} & = \\frac{\\tau_{\\widehat{y}}}{2} \\left\\{\\left[1 + \\mathrm{sgn} \\left(\\tau_{\\widehat{y}}\\right)\\right] {\\widetilde{\\phi}}_{i + 1 / 2}^\\mathrm{F} + \\left[1 - \\mathrm{sgn} \\left(\\tau_{\\widehat{y}}\\right)\\right] {\\widetilde{\\phi}}_{i + 1 / 2, j + 1}^\\mathrm{B}\\right\\},\\\\
112	\\mathcal{F}^{\\rho u, \\widehat{z}}_{i + 1 / 2, k + 1 / 2} & = \\frac{\\tau_{\\widehat{z}}}{2} \\left\\{\\left[1 + \\mathrm{sgn} \\left(\\tau_{\\widehat{z}}\\right)\\right] {\\widetilde{\\phi}}_{i + 1 / 2}^\\mathrm{U} + \\left[1 - \\mathrm{sgn} \\left(\\tau_{\\widehat{z}}\\right)\\right] {\\widetilde{\\phi}}_{i + 1 / 2, k + 1}^\\mathrm{D}\\right\\},
113	\\end{align*}
114	```
115
116	with
117
118	```math
119	\\begin{align*}
120	\\tau_{\\widehat{x}} & = \\left[\\left(J P_\\mathrm{old}\\right)_{i + 1 / 2} u_{\\mathrm{old}, i + 1 / 2}\\right]_{i + 1},\\\\
121	\\tau_{\\widehat{y}} & = \\left[\\left(J P_\\mathrm{old}\\right)_{j + 1 / 2} v_{\\mathrm{old}, j + 1 / 2}\\right]_{i + 1 / 2, j + 1 / 2},\\\\
122	\\tau_{\\widehat{z}} & = \\left[\\left(J P_\\mathrm{old}\\right)_{k + 1 / 2} \\widehat{w}_{\\mathrm{old}, k + 1 / 2}\\right]_{i + 1 / 2, k + 1 / 2}
123	\\end{align*}
124	```
125
126	and ``\\widetilde{\\phi}`` being the reconstruction of ``\\rho_{i + 1 / 2} u_{i + 1 / 2} / P_{\\mathrm{old}, i + 1 / 2}``. If the viscosity is nonzero, the viscous parts (weighted by the Jacobian) are then added, i.e.
127
128	```math
129	\\begin{align*}
130	\\mathcal{F}^{\\rho u, \\widehat{x}}_{i + 1} & \\rightarrow \\mathcal{F}^{\\rho u, \\widehat{x}}_{i + 1} - \\eta_{i + 1} \\left(J \\widehat{\\Pi}^{11}\\right)_{i + 1},\\\\
131	\\mathcal{F}^{\\rho u, \\widehat{y}}_{i + 1 / 2, j + 1 / 2} & \\rightarrow \\mathcal{F}^{\\rho u, \\widehat{y}}_{i + 1 / 2, j + 1 / 2} - \\eta_{i + 1 / 2, j + 1 / 2} \\left(J \\widehat{\\Pi}^{12}\\right)_{i + 1 / 2, j + 1 / 2},\\\\
132	\\mathcal{F}^{\\rho u, \\widehat{z}}_{i + 1 / 2, k + 1 / 2} & \\rightarrow \\mathcal{F}^{\\rho u, \\widehat{z}}_{i + 1 / 2, k + 1 / 2} - \\eta_{i + 1 / 2, k + 1 / 2} \\left(J \\widehat{\\Pi}^{13}\\right)_{i + 1 / 2, k + 1 / 2}.
133	\\end{align*}
134	```
135
136	Finally, if the diffusivity ``\\mu`` is nonzero, the diffusive parts (weighted by the Jacobian) are added, i.e.
137
138	```math
139	\\begin{align*}
140	\\mathcal{F}^{\\rho u, \\widehat{x}}_{i + 1} & \\rightarrow \\mathcal{F}^{\\rho u, \\widehat{x}}_{i + 1} - \\mu_{i + 1} \\left[J \\widehat{\\left(\\boldsymbol{\\nabla} u\\right)}^{\\widehat{x}}\\right]_{i + 1},\\\\
141	\\mathcal{F}^{\\rho u, \\widehat{y}}_{i + 1 / 2, j + 1 / 2} & \\rightarrow \\mathcal{F}^{\\rho u, \\widehat{y}}_{i + 1 / 2, j + 1 / 2} - \\mu_{i + 1 / 2, j + 1 / 2} \\left[J \\widehat{\\left(\\boldsymbol{\\nabla} u\\right)}^{\\widehat{y}}\\right]_{i + 1 / 2, j + 1 / 2},\\\\
142	\\mathcal{F}^{\\rho u, \\widehat{z}}_{i + 1 / 2, k + 1 / 2} & \\rightarrow \\mathcal{F}^{\\rho u, \\widehat{z}}_{i + 1 / 2, k + 1 / 2} - \\mu_{i + 1 / 2, k + 1 / 2} \\left[J \\widehat{\\left(\\boldsymbol{\\nabla} u\\right)}^{\\widehat{z}}\\right]_{i + 1 / 2, k + 1 / 2}.
143	\\end{align*}
144	```
145
146	```julia
147	compute_fluxes!(state::State, old_predictands::Predictands, variable::V)
148	```
149
150	Compute the meridional-momentum fluxes in all three directions.
151
152	The fluxes are first set to the advective parts
153
154	```math
155	\\begin{align*}
156	\\mathcal{F}^{\\rho v, \\widehat{x}}_{i + 1 / 2, j + 1 / 2} & = \\frac{\\tau_{\\widehat{x}}}{2} \\left\\{\\left[1 + \\mathrm{sgn} \\left(\\tau_{\\widehat{x}}\\right)\\right] {\\widetilde{\\phi}}_{j + 1 / 2}^\\mathrm{R} + \\left[1 - \\mathrm{sgn} \\left(\\tau_{\\widehat{x}}\\right)\\right] {\\widetilde{\\phi}}_{i + 1, j + 1 / 2}^\\mathrm{L}\\right\\},\\\\
157	\\mathcal{F}^{\\rho v, \\widehat{y}}_{j + 1} & = \\frac{\\tau_{\\widehat{y}}}{2} \\left\\{\\left[1 + \\mathrm{sgn} \\left(\\tau_{\\widehat{y}}\\right)\\right] {\\widetilde{\\phi}}_{j + 1 / 2}^\\mathrm{F} + \\left[1 - \\mathrm{sgn} \\left(\\tau_{\\widehat{y}}\\right)\\right] {\\widetilde{\\phi}}_{j + 3 / 2}^\\mathrm{B}\\right\\},\\\\
158	\\mathcal{F}^{\\rho v, \\widehat{z}}_{j + 1 / 2, k + 1 / 2} & = \\frac{\\tau_{\\widehat{z}}}{2} \\left\\{\\left[1 + \\mathrm{sgn} \\left(\\tau_{\\widehat{z}}\\right)\\right] {\\widetilde{\\phi}}_{j + 1 / 2}^\\mathrm{U} + \\left[1 - \\mathrm{sgn} \\left(\\tau_{\\widehat{z}}\\right)\\right] {\\widetilde{\\phi}}_{j + 1 / 2, k + 1}^\\mathrm{D}\\right\\},
159	\\end{align*}
160	```
161
162	with
163
164	```math
165	\\begin{align*}
166	\\tau_{\\widehat{x}} & = \\left[\\left(J P_\\mathrm{old}\\right)_{i + 1 / 2} u_{\\mathrm{old}, i + 1 / 2}\\right]_{i + 1 / 2, j + 1 / 2},\\\\
167	\\tau_{\\widehat{y}} & = \\left[\\left(J P_\\mathrm{old}\\right)_{j + 1 / 2} v_{\\mathrm{old}, j + 1 / 2}\\right]_{j + 1},\\\\
168	\\tau_{\\widehat{z}} & = \\left[\\left(J P_\\mathrm{old}\\right)_{k + 1 / 2} \\widehat{w}_{\\mathrm{old}, k + 1 / 2}\\right]_{j + 1 / 2, k + 1 / 2}
169	\\end{align*}
170	```
171
172	and ``\\widetilde{\\phi}`` being the reconstruction of ``\\rho_{j + 1 / 2} v_{j + 1 / 2} / P_{\\mathrm{old}, j + 1 / 2}``. If the viscosity is nonzero, the viscous parts (weighted by the Jacobian) are then added, i.e.
173
174	```math
175	\\begin{align*}
176	\\mathcal{F}^{\\rho v, \\widehat{x}}_{i + 1 / 2, j + 1 / 2} & \\rightarrow \\mathcal{F}^{\\rho v, \\widehat{x}}_{i + 1 / 2, j + 1 / 2} - \\eta_{i + 1 / 2, j + 1 / 2} \\left(J \\widehat{\\Pi}^{12}\\right)_{i + 1 / 2, j + 1 / 2},\\\\
177	\\mathcal{F}^{\\rho v, \\widehat{y}}_{j + 1} & \\rightarrow \\mathcal{F}^{\\rho v, \\widehat{y}}_{j + 1} - \\eta_{j + 1} \\left(J \\widehat{\\Pi}^{22}\\right)_{j + 1},\\\\
178	\\mathcal{F}^{\\rho v, \\widehat{z}}_{j + 1 / 2, k + 1 / 2} & \\rightarrow \\mathcal{F}^{\\rho v, \\widehat{z}}_{j + 1 / 2, k + 1 / 2} - \\eta_{j + 1 / 2, k + 1 / 2} \\left(J \\widehat{\\Pi}^{23}\\right)_{j + 1 / 2, k + 1 / 2}.
179	\\end{align*}
180	```
181
182	Finally, if the diffusivity ``\\mu`` is nonzero, the diffusive parts (weighted by the Jacobian) are added, i.e.
183
184	```math
185	\\begin{align*}
186	\\mathcal{F}^{\\rho v, \\widehat{x}}_{i + 1 / 2, j + 1 / 2} & \\rightarrow \\mathcal{F}^{\\rho v, \\widehat{x}}_{i + 1 / 2, j + 1 / 2} - \\mu_{i + 1 / 2, j + 1 / 2} \\left[J \\widehat{\\left(\\boldsymbol{\\nabla} v\\right)}^{\\widehat{x}}\\right]_{i + 1 / 2, j + 1 / 2},\\\\
187	\\mathcal{F}^{\\rho v, \\widehat{y}}_{j + 1} & \\rightarrow \\mathcal{F}^{\\rho v, \\widehat{y}}_{j + 1} - \\mu_{j + 1} \\left[J \\widehat{\\left(\\boldsymbol{\\nabla} v\\right)}^{\\widehat{y}}\\right]_{j + 1},\\\\
188	\\mathcal{F}^{\\rho v, \\widehat{z}}_{j + 1 / 2, k + 1 / 2} & \\rightarrow \\mathcal{F}^{\\rho v, \\widehat{z}}_{j + 1 / 2, k + 1 / 2} - \\mu_{j + 1 / 2, k + 1 / 2} \\left[J \\widehat{\\left(\\boldsymbol{\\nabla} v\\right)}^{\\widehat{z}}\\right]_{j + 1 / 2, k + 1 / 2}.
189	\\end{align*}
190	```
191
192	```julia
193	compute_fluxes!(state::State, old_predictands::Predictands, variable::W)
194	```
195
196	Compute the vertical-momentum fluxes in all three directions.
197
198	The fluxes are first set to the advective parts
199
200	```math
201	\\begin{align*}
202	\\mathcal{F}^{\\rho w, \\widehat{x}}_{i + 1 / 2, k + 1 / 2} & = \\frac{\\tau_{\\widehat{x}}}{2} \\left\\{\\left[1 + \\mathrm{sgn} \\left(\\tau_{\\widehat{x}}\\right)\\right] {\\widetilde{\\phi}}_{k + 1 / 2}^\\mathrm{R} + \\left[1 - \\mathrm{sgn} \\left(\\tau_{\\widehat{x}}\\right)\\right] {\\widetilde{\\phi}}_{i + 1, k + 1 / 2}^\\mathrm{L}\\right\\},\\\\
203	\\mathcal{F}^{\\rho w, \\widehat{y}}_{j + 1 / 2, k + 1 / 2} & = \\frac{\\tau_{\\widehat{y}}}{2} \\left\\{\\left[1 + \\mathrm{sgn} \\left(\\tau_{\\widehat{y}}\\right)\\right] {\\widetilde{\\phi}}_{k + 1 / 2}^\\mathrm{F} + \\left[1 - \\mathrm{sgn} \\left(\\tau_{\\widehat{y}}\\right)\\right] {\\widetilde{\\phi}}_{j + 1, k + 1 / 2}^\\mathrm{B}\\right\\},\\\\
204	\\mathcal{F}^{\\rho w, \\widehat{z}}_{k + 1} & = \\frac{\\tau_{\\widehat{z}}}{2} \\left\\{\\left[1 + \\mathrm{sgn} \\left(\\tau_{\\widehat{z}}\\right)\\right] {\\widetilde{\\phi}}_{k + 1 / 2}^\\mathrm{U} + \\left[1 - \\mathrm{sgn} \\left(\\tau_{\\widehat{z}}\\right)\\right] {\\widetilde{\\phi}}_{k + 3 / 2}^\\mathrm{D}\\right\\},
205	\\end{align*}
206	```
207
208	with
209
210	```math
211	\\begin{align*}
212	\\tau_{\\widehat{x}} & = \\left[\\left(J P_\\mathrm{old}\\right)_{i + 1 / 2} u_{\\mathrm{old}, i + 1 / 2}\\right]_{i + 1 / 2, k + 1 / 2},\\\\
213	\\tau_{\\widehat{y}} & = \\left[\\left(J P_\\mathrm{old}\\right)_{j + 1 / 2} v_{\\mathrm{old}, j + 1 / 2}\\right]_{j + 1 / 2, k + 1 / 2},\\\\
214	\\tau_{\\widehat{z}} & = \\left[\\left(J P_\\mathrm{old}\\right)_{k + 1 / 2} \\widehat{w}_{\\mathrm{old}, k + 1 / 2}\\right]_{k + 1}
215	\\end{align*}
216	```
217
218	and ``\\widetilde{\\phi}`` being the reconstruction of ``\\rho_{k + 1 / 2} w_{k + 1 / 2} / P_{\\mathrm{old}, k + 1 / 2}``. If the viscosity is nonzero, the viscous parts (weighted by the Jacobian) are then added, i.e.
219
220	```math
221	\\begin{align*}
222	\\mathcal{F}^{\\rho w, \\widehat{x}}_{i + 1 / 2, k + 1 / 2} & \\rightarrow \\mathcal{F}^{\\rho w, \\widehat{x}}_{i + 1 / 2, k + 1 / 2} - \\eta_{i + 1 / 2, k + 1 / 2} \\left(J \\Pi^{13}\\right)_{i + 1 / 2, k + 1 / 2},\\\\
223	\\mathcal{F}^{\\rho w, \\widehat{y}}_{j + 1 / 2, k + 1 / 2} & \\rightarrow \\mathcal{F}^{\\rho w, \\widehat{y}}_{j + 1 / 2, k + 1 / 2} - \\eta_{j + 1 / 2, k + 1 / 2} \\left(J \\Pi^{23}\\right)_{j + 1 / 2, k + 1 / 2},\\\\
224	\\mathcal{F}^{\\rho w, \\widehat{z}}_{k + 1} & \\rightarrow \\mathcal{F}^{\\rho w, \\widehat{z}}_{k + 1} - \\eta_{k + 1} \\left[\\left(J G^{13} \\Pi^{13}\\right)_{k + 1} - \\left(J G^{23} \\Pi^{23}\\right)_{k + 1} - \\Pi^{33}_{k + 1}\\right].
225	\\end{align*}
226	```
227
228	Finally, if the diffusivity ``\\mu`` is nonzero, the diffusive parts (weighted by the Jacobian) are added, i.e.
229
230	```math
231	\\begin{align*}
232	\\mathcal{F}^{\\rho w, \\widehat{x}}_{i + 1 / 2, k + 1 / 2} & \\rightarrow \\mathcal{F}^{\\rho w, \\widehat{x}}_{i + 1 / 2, k + 1 / 2} - \\mu_{i + 1 / 2, k + 1 / 2} \\left[J \\widehat{\\left(\\boldsymbol{\\nabla} w\\right)}^{\\widehat{x}}\\right]_{i + 1 / 2, k + 1 / 2},\\\\
233	\\mathcal{F}^{\\rho w, \\widehat{y}}_{j + 1 / 2, k + 1 / 2} & \\rightarrow \\mathcal{F}^{\\rho w, \\widehat{y}}_{j + 1 / 2, k + 1 / 2} - \\mu_{j + 1 / 2, k + 1 / 2} \\left[J \\widehat{\\left(\\boldsymbol{\\nabla} w\\right)}^{\\widehat{y}}\\right]_{j + 1 / 2, k + 1 / 2},\\\\
234	\\mathcal{F}^{\\rho w, \\widehat{z}}_{k + 1} & \\rightarrow \\mathcal{F}^{\\rho w, \\widehat{z}}_{k + 1} - \\mu_{k + 1} \\left[J \\widehat{\\left(\\boldsymbol{\\nabla} w\\right)}^{\\widehat{z}}\\right]_{k + 1}.
235	\\end{align*}
236	```
237
238	```julia
239	compute_fluxes!(state::State, predictands::Predictands, tracer_setup::NoTracer)
240	```
241
242	Return for configurations without tracer transport.
243
244	```julia
245	compute_fluxes!(state::State, predictands::Predictands, tracer_setup::TracerOn)
246	```
247
248	Compute the tracer fluxes in all three directions.
249
250	The computation is analogous to that of the density fluxes.
251
252	```julia
253	compute_fluxes!(state::State, predictands::Predictands, variable::Theta)
254	```
255
256	Compute the potential temperature fluxes due to heat conduction (weighted by the Jacobian).
257
258	The fluxes are given by
259
260	```math
261	\\begin{align*}
262	\\mathcal{F}^{\\theta, \\widehat{x}}_{i + 1 / 2} & = - \\lambda_{i + 1 / 2} \\left\\{\\frac{J_{i + 1 / 2}}{\\Delta \\widehat{x}} \\left[\\left(\\frac{P}{\\rho}\\right)_{i + 1} - \\frac{P}{\\rho}\\right]\\right.\\\\
263	& \\qquad \\qquad \\qquad + \\left.\\frac{\\left(J G^{13}\\right)_{i + 1 / 2}}{2 \\Delta \\widehat{z}} \\left[\\left(\\frac{P}{\\rho}\\right)_{i + 1 / 2, k + 1} - \\left(\\frac{P}{\\rho}\\right)_{i + 1 / 2, k - 1}\\right]\\right\\},\\\\
264	\\mathcal{F}^{\\theta, \\widehat{y}}_{j + 1 / 2} & = - \\lambda_{j + 1 / 2} \\left\\{\\frac{J_{j + 1 / 2}}{\\Delta \\widehat{y}} \\left[\\left(\\frac{P}{\\rho}\\right)_{j + 1} - \\frac{P}{\\rho}\\right]\\right.\\\\
265	& \\qquad \\qquad \\qquad + \\left.\\frac{\\left(J G^{23}\\right)_{j + 1 / 2}}{2 \\Delta \\widehat{z}} \\left[\\left(\\frac{P}{\\rho}\\right)_{j + 1 / 2, k + 1} - \\left(\\frac{P}{\\rho}\\right)_{j + 1 / 2, k - 1}\\right]\\right\\},\\\\
266	\\mathcal{F}^{\\theta, \\widehat{z}}_{k + 1 / 2} & = - \\lambda_{k + 1 / 2} \\left\\{\\frac{\\left(J G^{13}\\right)_{k + 1 / 2}}{2 \\Delta \\widehat{x}} \\left[\\left(\\frac{P}{\\rho}\\right)_{i + 1, k + 1 / 2} - \\left(\\frac{P}{\\rho}\\right)_{i - 1, k + 1 / 2}\\right]\\right.\\\\
267	& \\qquad \\qquad \\qquad + \\frac{\\left(J G^{23}\\right)_{k + 1 / 2}}{2 \\Delta \\widehat{y}} \\left[\\left(\\frac{P}{\\rho}\\right)_{j + 1, k + 1 / 2} - \\left(\\frac{P}{\\rho}\\right)_{j - 1, k + 1 / 2}\\right]\\\\
268	& \\qquad \\qquad \\qquad + \\left.\\frac{\\left(J G^{33}\\right)_{k + 1 / 2}}{\\Delta \\widehat{z}} \\left[\\left(\\frac{P}{\\rho}\\right)_{k + 1} - \\frac{P}{\\rho}\\right]\\right\\},
269	\\end{align*}
270	```
271
272	where ``\\lambda`` is the thermal conductivity (computed from `state.namelists.atmosphere.thermal_conductivity`).
273
274	# Arguments
275
276	- `state`: Model state.
277
278	- `predictands`/`old_predictands`: The predictands that are used to compute the transporting velocities.
279
280	- `model`: Dynamic equations.
281
282	- `variable`: Flux variable.
283
284	- `tracer_setup`: General tracer-transport configuration.
285
286	# See also
287
288	- [`PinCFlow.FluxCalculator.compute_flux`](@ref)
289
290	- [`PinCFlow.Update.compute_stress_tensor`](@ref)
291
292	- [`PinCFlow.Update.conductive_heating`](@ref)
293	"""
294	function compute_fluxes! end
295
296	function compute_fluxes!(state::State, predictands::Predictands)	1,044✔
297	(; model) = state.namelists.atmosphere	1,044✔
298
299	compute_fluxes!(state, predictands, Rho())	1,044✔
300	compute_fluxes!(state, predictands, RhoP())	1,044✔
301	compute_fluxes!(state, predictands, U())	1,044✔
302	compute_fluxes!(state, predictands, V())	1,044✔
303	compute_fluxes!(state, predictands, W())	1,044✔
304
305	compute_fluxes!(state, predictands, model, P())	1,044✔
306	compute_fluxes!(state, predictands, state.namelists.tracer.tracer_setup)	1,044✔
307	return	1,044✔
308	end
309
310	function compute_fluxes!(state::State, predictands::Predictands, variable::Rho)	1,044✔
311	(; model) = state.namelists.atmosphere	1,044✔
312	compute_fluxes!(state, predictands, variable, model)	1,044✔
313	return	1,044✔
314	end
315
NEW 316	function compute_fluxes!(	×
317	state::State,
318	predictands::Predictands,
319	variable::Rho,
320	model::Boussinesq,
321	)
NEW 322	return	×
323	end
324
325	function compute_fluxes!(	1,044✔
326	state::State,
327	predictands::Predictands,
328	variable::Rho,
329	model::Union{PseudoIncompressible, Compressible},
330	)
331	(; i0, i1, j0, j1, k0, k1) = state.domain	1,044✔
332	(; jac) = state.grid	1,044✔
333	(; pbar, rhobar) = state.atmosphere	1,044✔
334	(; rhotilde) = state.variables.reconstructions	1,044✔
335	(; phirho) = state.variables.fluxes	1,044✔
336
337	(u0, v0, w0) = (predictands.u, predictands.v, predictands.w)	1,044✔
338
339	#-----------------------------------------
340	# Zonal fluxes
341	#-----------------------------------------
342
343	@ivy for k in k0:k1, j in j0:j1, i in (i0 - 1):i1	10,440✔
344	rhobaredger = 0.5 * (rhobar[i, j, k] + rhobar[i + 1, j, k])	708,840✔
345	pedger = 0.5 * (pbar[i, j, k] + pbar[i + 1, j, k])	708,840✔
346	rhor = rhotilde[i + 1, j, k, 1, 1] + rhobaredger / pedger	708,840✔
347	rhol = rhotilde[i, j, k, 1, 2] + rhobaredger / pedger	708,840✔
348
349	pedger =	708,840✔
350	0.5 * (
351	jac[i, j, k] * pbar[i, j, k] +
352	jac[i + 1, j, k] * pbar[i + 1, j, k]
353	)
354	usurf = pedger * u0[i, j, k]	708,840✔
355
356	frho = compute_flux(usurf, rhol, rhor)	708,840✔
357
358	phirho[i, j, k, 1] = frho	708,840✔
359	end	×
360
361	#-----------------------------------------
362	# Meridional fluxes
363	#-----------------------------------------
364
365	@ivy for k in k0:k1, j in (j0 - 1):j1, i in i0:i1	10,440✔
366	rhobaredgef = 0.5 * (rhobar[i, j, k] + rhobar[i, j + 1, k])	748,800✔
367	pedgef = 0.5 * (pbar[i, j, k] + pbar[i, j + 1, k])	748,800✔
368	rhof = rhotilde[i, j + 1, k, 2, 1] + rhobaredgef / pedgef	748,800✔
369	rhob = rhotilde[i, j, k, 2, 2] + rhobaredgef / pedgef	748,800✔
370
371	pedgef =	748,800✔
372	0.5 * (
373	jac[i, j, k] * pbar[i, j, k] +
374	jac[i, j + 1, k] * pbar[i, j + 1, k]
375	)
376	vsurf = pedgef * v0[i, j, k]	748,800✔
377
378	grho = compute_flux(vsurf, rhob, rhof)	1,171,887✔
379
380	phirho[i, j, k, 2] = grho	748,800✔
381	end	×
382
383	#-----------------------------------------
384	# Vertical fluxes
385	#-----------------------------------------
386
387	@ivy for k in (k0 - 1):k1, j in j0:j1, i in i0:i1	11,484✔
388	rhobaredgeu =	708,840✔
389	(
390	jac[i, j, k + 1] * rhobar[i, j, k] +
391	jac[i, j, k] * rhobar[i, j, k + 1]
392	) / (jac[i, j, k] + jac[i, j, k + 1])
393	pedgeu =	708,840✔
394	(
395	jac[i, j, k + 1] * pbar[i, j, k] +
396	jac[i, j, k] * pbar[i, j, k + 1]
397	) / (jac[i, j, k] + jac[i, j, k + 1])
398	rhou = rhotilde[i, j, k + 1, 3, 1] + rhobaredgeu / pedgeu	708,840✔
399	rhod = rhotilde[i, j, k, 3, 2] + rhobaredgeu / pedgeu	708,840✔
400
401	pedgeu =	708,840✔
402	jac[i, j, k] *
403	jac[i, j, k + 1] *
404	(pbar[i, j, k] + pbar[i, j, k + 1]) /
405	(jac[i, j, k] + jac[i, j, k + 1])
406	wsurf = pedgeu * w0[i, j, k]	708,840✔
407
408	hrho = compute_flux(wsurf, rhod, rhou)	1,136,094✔
409
410	phirho[i, j, k, 3] = hrho	708,840✔
411	end	×
412
413	return	1,044✔
414	end
415
416	function compute_fluxes!(state::State, predictands::Predictands, variable::RhoP)	1,044✔
417	(; i0, i1, j0, j1, k0, k1) = state.domain	1,044✔
418	(; jac) = state.grid	1,044✔
419	(; pbar) = state.atmosphere	1,044✔
420	(; rhoptilde) = state.variables.reconstructions	1,044✔
421	(; phirhop) = state.variables.fluxes	1,044✔
422
423	(u0, v0, w0) = (predictands.u, predictands.v, predictands.w)	1,044✔
424
425	#-----------------------------------------
426	# Zonal fluxes
427	#-----------------------------------------
428
429	@ivy for k in k0:k1, j in j0:j1, i in (i0 - 1):i1	10,440✔
430	rhor = rhoptilde[i + 1, j, k, 1, 1]	708,840✔
431	rhol = rhoptilde[i, j, k, 1, 2]	708,840✔
432
433	pedger =	708,840✔
434	0.5 * (
435	jac[i, j, k] * pbar[i, j, k] +
436	jac[i + 1, j, k] * pbar[i + 1, j, k]
437	)
438	usurf = pedger * u0[i, j, k]	708,840✔
439
440	frhop = compute_flux(usurf, rhol, rhor)	708,840✔
441
442	phirhop[i, j, k, 1] = frhop	708,840✔
443	end	×
444
445	#-----------------------------------------
446	# Meridional fluxes
447	#-----------------------------------------
448
449	@ivy for k in k0:k1, j in (j0 - 1):j1, i in i0:i1	10,440✔
450	rhof = rhoptilde[i, j + 1, k, 2, 1]	748,800✔
451	rhob = rhoptilde[i, j, k, 2, 2]	748,800✔
452
453	pedgef =	748,800✔
454	0.5 * (
455	jac[i, j, k] * pbar[i, j, k] +
456	jac[i, j + 1, k] * pbar[i, j + 1, k]
457	)
458	vsurf = pedgef * v0[i, j, k]	748,800✔
459
460	grhop = compute_flux(vsurf, rhob, rhof)	1,171,887✔
461
462	phirhop[i, j, k, 2] = grhop	748,800✔
463	end	×
464
465	#-----------------------------------------
466	# Vertical fluxes
467	#-----------------------------------------
468
469	@ivy for k in (k0 - 1):k1, j in j0:j1, i in i0:i1	11,484✔
470	rhou = rhoptilde[i, j, k + 1, 3, 1]	708,840✔
471	rhod = rhoptilde[i, j, k, 3, 2]	708,840✔
472
473	pedgeu =	708,840✔
474	jac[i, j, k] *
475	jac[i, j, k + 1] *
476	(pbar[i, j, k] + pbar[i, j, k + 1]) /
477	(jac[i, j, k] + jac[i, j, k + 1])
478	wsurf = pedgeu * w0[i, j, k]	708,840✔
479
480	hrhop = compute_flux(wsurf, rhod, rhou)	1,136,094✔
481
482	phirhop[i, j, k, 3] = hrhop	708,840✔
483	end	×
484
485	return	1,044✔
486	end
487
488	function compute_fluxes!(	1,044✔
489	state::State,
490	predictands::Predictands,
491	model::Union{Boussinesq, PseudoIncompressible},
492	variable::P,
493	)
494	return	1,044✔
495	end
496
497	function compute_fluxes!(	×
498	state::State,
499	predictands::Predictands,
500	model::Compressible,
501	variable::P,
502	)
503	(; i0, i1, j0, j1, k0, k1) = state.domain	×
504	(; jac) = state.grid	×
505	(; pbar) = state.atmosphere	×
506	(; phip) = state.variables.fluxes	×
507
508	(u0, v0, w0) = (predictands.u, predictands.v, predictands.w)	×
509
510	#-----------------------------------------
511	# Zonal fluxes
512	#-----------------------------------------
513
514	@ivy for k in k0:k1, j in j0:j1, i in (i0 - 1):i1	×
515	phip[i, j, k, 1] =	×
516	0.5 *
517	(
518	jac[i, j, k] * pbar[i, j, k] +
519	jac[i + 1, j, k] * pbar[i + 1, j, k]
520	) *
521	u0[i, j, k]
522	end	×
523
524	#-----------------------------------------
525	# Meridional fluxes
526	#-----------------------------------------
527
528	@ivy for k in k0:k1, j in (j0 - 1):j1, i in i0:i1	×
529	phip[i, j, k, 2] =	×
530	0.5 *
531	(
532	jac[i, j, k] * pbar[i, j, k] +
533	jac[i, j + 1, k] * pbar[i, j + 1, k]
534	) *
535	v0[i, j, k]
536	end	×
537
538	#-----------------------------------------
539	# Vertical fluxes
540	#-----------------------------------------
541
542	@ivy for k in (k0 - 1):k1, j in j0:j1, i in i0:i1	×
543	phip[i, j, k, 3] =	×
544	jac[i, j, k] *
545	jac[i, j, k + 1] *
546	(pbar[i, j, k] + pbar[i, j, k + 1]) /
547	(jac[i, j, k] + jac[i, j, k + 1]) * w0[i, j, k]
548	end	×
549
550	return	×
551	end
552
553	function compute_fluxes!(	1,044✔
554	state::State,
555	old_predictands::Predictands,
556	variable::U,
557	)
558	(; grid) = state	1,044✔
559	(; re, uref, lref) = state.constants	1,044✔
560	(; zz_size, nzz, ko, i0, i1, j0, j1, k0, k1) = state.domain	1,044✔
561	(; jac, met) = grid	1,044✔
562	(; pbar, rhobar) = state.atmosphere	1,044✔
563	(; utilde) = state.variables.reconstructions	1,044✔
564	(; phiu) = state.variables.fluxes	1,044✔
565	(; kinematic_diffusivity) = state.namelists.atmosphere	1,044✔
566
567	(u0, v0, w0) = (old_predictands.u, old_predictands.v, old_predictands.w)	1,044✔
568
569	kmin = k0	1,044✔
570	kmax = ko + nzz == zz_size ? k1 : k1 + 1	1,044✔
571
572	#-----------------------------------------
573	# Zonal fluxes
574	#-----------------------------------------
575
576	@ivy for k in kmin:kmax, j in j0:j1, i in (i0 - 2):i1	10,440✔
577	ur = utilde[i + 1, j, k, 1, 1]	773,280✔
578	ul = utilde[i, j, k, 1, 2]	773,280✔
579
580	pedger =	773,280✔
581	0.5 * (
582	jac[i, j, k] * pbar[i, j, k] +
583	jac[i + 1, j, k] * pbar[i + 1, j, k]
584	)
585	predger =	773,280✔
586	0.5 * (
587	jac[i + 1, j, k] * pbar[i + 1, j, k] +
588	jac[i + 2, j, k] * pbar[i + 2, j, k]
589	)
590	usurf = 0.5 * (pedger * u0[i, j, k] + predger * u0[i + 1, j, k])	773,280✔
591
592	frhou = compute_flux(usurf, ul, ur)	773,280✔
593
594	phiu[i, j, k, 1] = frhou	773,280✔
595	end	×
596
597	#-----------------------------------------
598	# Meridional fluxes
599	#-----------------------------------------
600
601	@ivy for k in kmin:kmax, j in (j0 - 1):j1, i in (i0 - 1):i1	10,440✔
602	uf = utilde[i, j + 1, k, 2, 1]	823,680✔
603	ub = utilde[i, j, k, 2, 2]	823,680✔
604
605	pedgef =	823,680✔
606	0.5 * (
607	jac[i, j, k] * pbar[i, j, k] +
608	jac[i, j + 1, k] * pbar[i, j + 1, k]
609	)
610	predgef =	823,680✔
611	0.5 * (
612	jac[i + 1, j, k] * pbar[i + 1, j, k] +
613	jac[i + 1, j + 1, k] * pbar[i + 1, j + 1, k]
614	)
615	vsurf = 0.5 * (pedgef * v0[i, j, k] + predgef * v0[i + 1, j, k])	823,680✔
616
617	grhou = compute_flux(vsurf, ub, uf)	1,289,853✔
618
619	phiu[i, j, k, 2] = grhou	823,680✔
620	end	×
621
622	#-----------------------------------------
623	# Vertical fluxes
624	#-----------------------------------------
625
626	@ivy for k in (kmin - 1):kmax, j in j0:j1, i in (i0 - 1):i1	11,484✔
627	uu = utilde[i, j, k + 1, 3, 1]	779,724✔
628	ud = utilde[i, j, k, 3, 2]	779,724✔
629
630	pedgeu =	779,724✔
631	jac[i, j, k] *
632	jac[i, j, k + 1] *
633	(pbar[i, j, k] + pbar[i, j, k + 1]) /
634	(jac[i, j, k] + jac[i, j, k + 1])
635	predgeu =	779,724✔
636	jac[i + 1, j, k] *
637	jac[i + 1, j, k + 1] *
638	(pbar[i + 1, j, k] + pbar[i + 1, j, k + 1]) /
639	(jac[i + 1, j, k] + jac[i + 1, j, k + 1])
640	wsurf = 0.5 * (pedgeu * w0[i, j, k] + predgeu * w0[i + 1, j, k])	779,724✔
641
642	hrhou = compute_flux(wsurf, ud, uu)	1,251,156✔
643
644	phiu[i, j, k, 3] = hrhou	779,724✔
645	end	×
646
647	#-------------------------------------------------------------------
648	# Viscous fluxes
649	#-------------------------------------------------------------------
650
651	if 1 / re <= eps() && kinematic_diffusivity == 0.0	1,044✔
652	return	×
653	end
654
655	#-----------------------------------------
656	# Zonal fluxes
657	#-----------------------------------------
658
659	@ivy for k in kmin:kmax, j in j0:j1, i in (i0 - 2):i1	10,440✔
660	coef_v = 1 / re * rhobar[i + 1, j, k0]	773,280✔
661
662	frhou_visc =	773,280✔
663	coef_v *
664	jac[i + 1, j, k] *
665	compute_stress_tensor(i + 1, j, k, 1, 1, state)
666
667	phiu[i, j, k, 1] -= frhou_visc	773,280✔
668	end	×
669
670	#-----------------------------------------
671	# Meridional fluxes
672	#-----------------------------------------
673
674	@ivy for k in kmin:kmax, j in (j0 - 1):j1, i in (i0 - 1):i1	10,440✔
675	coef_v =	823,680✔
676	1 / re *
677	0.25 *
678	(
679	rhobar[i, j, k0] +
680	rhobar[i + 1, j, k0] +
681	rhobar[i, j + 1, k0] +
682	rhobar[i + 1, j + 1, k0]
683	)
684
685	grhou_visc =	823,680✔
686	coef_v *
687	0.25 *
688	(
689	jac[i, j, k] * compute_stress_tensor(i, j, k, 1, 2, state) +
690	jac[i + 1, j, k] *
691	compute_stress_tensor(i + 1, j, k, 1, 2, state) +
692	jac[i, j + 1, k] *
693	compute_stress_tensor(i, j + 1, k, 1, 2, state) +
694	jac[i + 1, j + 1, k] *
695	compute_stress_tensor(i + 1, j + 1, k, 1, 2, state)
696	)
697
698	phiu[i, j, k, 2] -= grhou_visc	823,680✔
699	end	×
700
701	#-----------------------------------------
702	# Vertical fluxes
703	#-----------------------------------------
704
705	@ivy for k in (kmin - 1):kmax, j in j0:j1, i in (i0 - 1):i1	11,484✔
706	coef_v = 1 / re * 0.5 * (rhobar[i, j, k0] + rhobar[i + 1, j, k0])	779,724✔
707
708	stresstens13 =	779,724✔
709	met[i, j, k, 1, 3] * compute_stress_tensor(i, j, k, 1, 1, state) +
710	met[i, j, k, 2, 3] * compute_stress_tensor(i, j, k, 1, 2, state) +
711	compute_stress_tensor(i, j, k, 1, 3, state) / jac[i, j, k]
712	stresstens13r =	779,724✔
713	met[i + 1, j, k, 1, 3] *
714	compute_stress_tensor(i + 1, j, k, 1, 1, state) +
715	met[i + 1, j, k, 2, 3] *
716	compute_stress_tensor(i + 1, j, k, 1, 2, state) +
717	compute_stress_tensor(i + 1, j, k, 1, 3, state) / jac[i + 1, j, k]
718	stresstens13u =	779,724✔
719	met[i, j, k + 1, 1, 3] *
720	compute_stress_tensor(i, j, k + 1, 1, 1, state) +
721	met[i, j, k + 1, 2, 3] *
722	compute_stress_tensor(i, j, k + 1, 1, 2, state) +
723	compute_stress_tensor(i, j, k + 1, 1, 3, state) / jac[i, j, k + 1]
724	stresstens13ru =	779,724✔
725	met[i + 1, j, k + 1, 1, 3] *
726	compute_stress_tensor(i + 1, j, k + 1, 1, 1, state) +
727	met[i + 1, j, k + 1, 2, 3] *
728	compute_stress_tensor(i + 1, j, k + 1, 1, 2, state) +
729	compute_stress_tensor(i + 1, j, k + 1, 1, 3, state) /
730	jac[i + 1, j, k + 1]
731	hrhou_visc =	779,724✔
732	coef_v *
733	0.5 *
734	(
735	jac[i, j, k] *
736	jac[i, j, k + 1] *
737	(stresstens13 + stresstens13u) /
738	(jac[i, j, k] + jac[i, j, k + 1]) +
739	jac[i + 1, j, k] *
740	jac[i + 1, j, k + 1] *
741	(stresstens13r + stresstens13ru) /
742	(jac[i + 1, j, k] + jac[i + 1, j, k + 1])
743	)
744
745	phiu[i, j, k, 3] -= hrhou_visc	779,724✔
746	end	×
747
748	#-------------------------------------------------------------------
749	# Diffusion fluxes
750	#-------------------------------------------------------------------
751
752	if kinematic_diffusivity == 0.0	1,044✔
753	return	1,044✔
754	end
755
756	mu_mom_diff = kinematic_diffusivity / uref / lref	×
757
758	#-----------------------------------------
759	# Zonal fluxes
760	#-----------------------------------------
761
762	@ivy for k in kmin:kmax, j in j0:j1, i in (i0 - 2):i1	×
763	coef_d = mu_mom_diff * rhobar[i + 1, j, k0]	×
764
765	frhou_diff =	×
766	coef_d *
767	jac[i + 1, j, k] *
768	compute_momentum_diffusion_terms(state, i + 1, j, k, U(), X())
769
770	phiu[i, j, k, 1] -= frhou_diff	×
771	end	×
772
773	#-----------------------------------------
774	# Meridional fluxes
775	#-----------------------------------------
776
777	@ivy for k in kmin:kmax, j in (j0 - 1):j1, i in (i0 - 1):i1	×
778	coef_d =	×
779	mu_mom_diff *
780	0.25 *
781	(
782	rhobar[i, j, k0] +
783	rhobar[i + 1, j, k0] +
784	rhobar[i, j + 1, k0] +
785	rhobar[i + 1, j + 1, k0]
786	)
787
788	grhou_diff =	×
789	coef_d *
790	0.25 *
791	(
792	jac[i, j, k] *
793	compute_momentum_diffusion_terms(state, i, j, k, U(), Y()) +
794	jac[i + 1, j, k] *
795	compute_momentum_diffusion_terms(state, i + 1, j, k, U(), Y()) +
796	jac[i, j + 1, k] *
797	compute_momentum_diffusion_terms(state, i, j + 1, k, U(), Y()) +
798	jac[i + 1, j + 1, k] * compute_momentum_diffusion_terms(
799	state,
800	i + 1,
801	j + 1,
802	k,
803	U(),
804	Y(),
805	)
806	)
807
808	phiu[i, j, k, 2] -= grhou_diff	×
809	end	×
810
811	#-----------------------------------------
812	# Vertical fluxes
813	#-----------------------------------------
814
815	@ivy for k in (kmin - 1):kmax, j in j0:j1, i in (i0 - 1):i1	×
816	coef_dr = mu_mom_diff * rhobar[i + 1, j, k0]	×
817
818	coef_dl = mu_mom_diff * rhobar[i, j, k0]	×
819
820	coef_d = 0.5 * (coef_dr + coef_dl)	×
821
822	mom_diff =	×
823	jac[i, j, k] *
824	compute_momentum_diffusion_terms(state, i, j, k, U(), Z())
825
826	mom_diff_r =	×
827	jac[i + 1, j, k] *
828	compute_momentum_diffusion_terms(state, i + 1, j, k, U(), Z())
829
830	mom_diff_u =	×
831	jac[i, j, k + 1] *
832	compute_momentum_diffusion_terms(state, i, j, k + 1, U(), Z())
833
834	mom_diff_ru =	×
835	jac[i + 1, j, k + 1] *
836	compute_momentum_diffusion_terms(state, i + 1, j, k + 1, U(), Z())
837
838	hrhou_diff =	×
839	coef_d *
840	0.5 *
841	(
842	jac[i, j, k] * jac[i, j, k + 1] * (mom_diff + mom_diff_u) /
843	(jac[i, j, k] + jac[i, j, k + 1]) +
844	jac[i + 1, j, k] *
845	jac[i + 1, j, k + 1] *
846	(mom_diff_r + mom_diff_ru) /
847	(jac[i + 1, j, k] + jac[i + 1, j, k + 1])
848	)
849
850	phiu[i, j, k, 3] -= hrhou_diff	×
851	end	×
852
853	return	×
854	end
855
856	function compute_fluxes!(	1,044✔
857	state::State,
858	old_predictands::Predictands,
859	variable::V,
860	)
861	(; grid) = state	1,044✔
862	(; re, uref, lref) = state.constants	1,044✔
863	(; zz_size, nzz, ko, i0, i1, j0, j1, k0, k1) = state.domain	1,044✔
864	(; jac, met) = grid	1,044✔
865	(; pbar, rhobar) = state.atmosphere	1,044✔
866	(; vtilde) = state.variables.reconstructions	1,044✔
867	(; phiv) = state.variables.fluxes	1,044✔
868	(; kinematic_diffusivity) = state.namelists.atmosphere	1,044✔
869
870	(u0, v0, w0) = (old_predictands.u, old_predictands.v, old_predictands.w)	1,044✔
871
872	kmin = k0	1,044✔
873	kmax = ko + nzz == zz_size ? k1 : k1 + 1	1,044✔
874
875	#-----------------------------------------
876	# Zonal fluxes
877	#-----------------------------------------
878
879	@ivy for k in kmin:kmax, j in (j0 - 1):j1, i in (i0 - 1):i1	10,440✔
880	vr = vtilde[i + 1, j, k, 1, 1]	823,680✔
881	vl = vtilde[i, j, k, 1, 2]	823,680✔
882
883	pedger =	823,680✔
884	0.5 * (
885	jac[i, j, k] * pbar[i, j, k] +
886	jac[i + 1, j, k] * pbar[i + 1, j, k]
887	)
888	pfedger =	823,680✔
889	0.5 * (
890	jac[i, j + 1, k] * pbar[i, j + 1, k] +
891	jac[i + 1, j + 1, k] * pbar[i + 1, j + 1, k]
892	)
893	usurf = 0.5 * (pedger * u0[i, j, k] + pfedger * u0[i, j + 1, k])	823,680✔
894
895	frhov = compute_flux(usurf, vl, vr)	823,680✔
896
897	phiv[i, j, k, 1] = frhov	823,680✔
898	end	×
899
900	#-----------------------------------------
901	# Meridional fluxes
902	#-----------------------------------------
903
904	@ivy for k in kmin:kmax, j in (j0 - 2):j1, i in i0:i1	10,440✔
905	vf = vtilde[i, j + 1, k, 2, 1]	853,200✔
906	vb = vtilde[i, j, k, 2, 2]	853,200✔
907
908	pedgef =	853,200✔
909	0.5 * (
910	jac[i, j, k] * pbar[i, j, k] +
911	jac[i, j + 1, k] * pbar[i, j + 1, k]
912	)
913	pfedgef =	853,200✔
914	0.5 * (
915	jac[i, j + 1, k] * pbar[i, j + 1, k] +
916	jac[i, j + 2, k] * pbar[i, j + 2, k]
917	)
918	vsurf = 0.5 * (pedgef * v0[i, j, k] + pfedgef * v0[i, j + 1, k])	853,200✔
919
920	grhov = compute_flux(vsurf, vb, vf)	1,346,400✔
921
922	phiv[i, j, k, 2] = grhov	853,200✔
923	end	×
924
925	#-----------------------------------------
926	# Vertical fluxes
927	#-----------------------------------------
928
929	@ivy for k in (kmin - 1):kmax, j in (j0 - 1):j1, i in i0:i1	11,484✔
930	vu = vtilde[i, j, k + 1, 3, 1]	823,680✔
931	vd = vtilde[i, j, k, 3, 2]	823,680✔
932
933	pedgeu =	823,680✔
934	jac[i, j, k] *
935	jac[i, j, k + 1] *
936	(pbar[i, j, k] + pbar[i, j, k + 1]) /
937	(jac[i, j, k] + jac[i, j, k + 1])
938	pfedgeu =	823,680✔
939	jac[i, j + 1, k] *
940	jac[i, j + 1, k + 1] *
941	(pbar[i, j + 1, k] + pbar[i, j + 1, k + 1]) /
942	(jac[i, j + 1, k] + jac[i, j + 1, k + 1])
943	wsurf = 0.5 * (pedgeu * w0[i, j, k] + pfedgeu * w0[i, j + 1, k])	823,680✔
944
945	hrhov = compute_flux(wsurf, vd, vu)	1,318,446✔
946
947	phiv[i, j, k, 3] = hrhov	823,680✔
948	end	×
949
950	#-------------------------------------------------------------------
951	# Viscous fluxes
952	#-------------------------------------------------------------------
953
954	if 1 / re <= eps() && kinematic_diffusivity == 0.0	1,044✔
955	return	×
956	end
957
958	#-----------------------------------------
959	# Zonal fluxes
960	#-----------------------------------------
961
962	@ivy for k in kmin:kmax, j in (j0 - 1):j1, i in (i0 - 1):i1	10,440✔
963	coef_v =	823,680✔
964	1 / re *
965	0.25 *
966	(
967	rhobar[i, j, k0] +
968	rhobar[i + 1, j, k0] +
969	rhobar[i, j + 1, k0] +
970	rhobar[i + 1, j + 1, k0]
971	)
972
973	frhov_visc =	823,680✔
974	coef_v *
975	0.25 *
976	(
977	jac[i, j, k] * compute_stress_tensor(i, j, k, 2, 1, state) +
978	jac[i + 1, j, k] *
979	compute_stress_tensor(i + 1, j, k, 2, 1, state) +
980	jac[i, j + 1, k] *
981	compute_stress_tensor(i, j + 1, k, 2, 1, state) +
982	jac[i + 1, j + 1, k] *
983	compute_stress_tensor(i + 1, j + 1, k, 2, 1, state)
984	)
985
986	phiv[i, j, k, 1] -= frhov_visc	823,680✔
987	end	×
988
989	#-----------------------------------------
990	# Meridional fluxes
991	#-----------------------------------------
992
993	@ivy for k in kmin:kmax, j in (j0 - 2):j1, i in i0:i1	10,440✔
994	coef_v = 1 / re * rhobar[i, j + 1, k0]	853,200✔
995
996	grhov_visc =	853,200✔
997	coef_v *
998	jac[i, j + 1, k] *
999	compute_stress_tensor(i, j + 1, k, 2, 2, state)
1000
1001	phiv[i, j, k, 2] -= grhov_visc	853,200✔
1002	end	×
1003
1004	#-----------------------------------------
1005	# Vertical fluxes
1006	#-----------------------------------------
1007
1008	@ivy for k in (kmin - 1):kmax, j in (j0 - 1):j1, i in i0:i1	11,484✔
1009	coef_v = 1 / re * 0.5 * (rhobar[i, j, k0] + rhobar[i, j + 1, k0])	823,680✔
1010
1011	stresstens23 =	823,680✔
1012	met[i, j, k, 1, 3] * compute_stress_tensor(i, j, k, 2, 1, state) +
1013	met[i, j, k, 2, 3] * compute_stress_tensor(i, j, k, 2, 2, state) +
1014	compute_stress_tensor(i, j, k, 2, 3, state) / jac[i, j, k]
1015	stresstens23f =	823,680✔
1016	met[i, j + 1, k, 1, 3] *
1017	compute_stress_tensor(i, j + 1, k, 2, 1, state) +
1018	met[i, j + 1, k, 2, 3] *
1019	compute_stress_tensor(i, j + 1, k, 2, 2, state) +
1020	compute_stress_tensor(i, j + 1, k, 2, 3, state) / jac[i, j + 1, k]
1021	stresstens23u =	823,680✔
1022	met[i, j, k + 1, 1, 3] *
1023	compute_stress_tensor(i, j, k + 1, 2, 1, state) +
1024	met[i, j, k + 1, 2, 3] *
1025	compute_stress_tensor(i, j, k + 1, 2, 2, state) +
1026	compute_stress_tensor(i, j, k + 1, 2, 3, state) / jac[i, j, k + 1]
1027	stresstens23fu =	823,680✔
1028	met[i, j + 1, k + 1, 1, 3] *
1029	compute_stress_tensor(i, j + 1, k + 1, 2, 1, state) +
1030	met[i, j + 1, k + 1, 2, 3] *
1031	compute_stress_tensor(i, j + 1, k + 1, 2, 2, state) +
1032	compute_stress_tensor(i, j + 1, k + 1, 2, 3, state) /
1033	jac[i, j + 1, k + 1]
1034	hrhov_visc =	823,680✔
1035	coef_v *
1036	0.5 *
1037	(
1038	jac[i, j, k] *
1039	jac[i, j, k + 1] *
1040	(stresstens23 + stresstens23u) /
1041	(jac[i, j, k] + jac[i, j, k + 1]) +
1042	jac[i, j + 1, k] *
1043	jac[i, j + 1, k + 1] *
1044	(stresstens23f + stresstens23fu) /
1045	(jac[i, j + 1, k] + jac[i, j + 1, k + 1])
1046	)
1047
1048	phiv[i, j, k, 3] -= hrhov_visc	823,680✔
1049	end	×
1050
1051	#-------------------------------------------------------------------
1052	# Diffusion fluxes
1053	#-------------------------------------------------------------------
1054
1055	if kinematic_diffusivity == 0.0	1,044✔
1056	return	1,044✔
1057	end
1058
1059	mu_mom_diff = kinematic_diffusivity / uref / lref	×
1060
1061	#-----------------------------------------
1062	# Zonal fluxes
1063	#-----------------------------------------
1064
1065	@ivy for k in kmin:kmax, j in (j0 - 1):j1, i in (i0 - 1):i1	×
1066	coef_d =	×
1067	mu_mom_diff *
1068	0.25 *
1069	(
1070	rhobar[i, j, k0] +
1071	rhobar[i + 1, j, k0] +
1072	rhobar[i, j + 1, k0] +
1073	rhobar[i + 1, j + 1, k0]
1074	)
1075
1076	frhov_diff =	×
1077	coef_d *
1078	0.25 *
1079	(
1080	jac[i, j, k] *
1081	compute_momentum_diffusion_terms(state, i, j, k, V(), X()) +
1082	jac[i + 1, j, k] *
1083	compute_momentum_diffusion_terms(state, i + 1, j, k, V(), X()) +
1084	jac[i, j + 1, k] *
1085	compute_momentum_diffusion_terms(state, i, j + 1, k, V(), X()) +
1086	jac[i + 1, j + 1, k] * compute_momentum_diffusion_terms(
1087	state,
1088	i + 1,
1089	j + 1,
1090	k,
1091	V(),
1092	X(),
1093	)
1094	)
1095
1096	phiv[i, j, k, 1] -= frhov_diff	×
1097	end	×
1098
1099	#-----------------------------------------
1100	# Meridional fluxes
1101	#-----------------------------------------
1102
1103	@ivy for k in kmin:kmax, j in (j0 - 2):j1, i in i0:i1	×
1104	coef_d = mu_mom_diff * rhobar[i, j + 1, k0]	×
1105
1106	grhov_diff =	×
1107	coef_d *
1108	jac[i, j + 1, k] *
1109	compute_momentum_diffusion_terms(state, i, j + 1, k, V(), Y())
1110
1111	phiv[i, j, k, 2] -= grhov_diff	×
1112	end	×
1113
1114	#-----------------------------------------
1115	# Vertical fluxes
1116	#-----------------------------------------
1117
1118	@ivy for k in (kmin - 1):kmax, j in (j0 - 1):j1, i in i0:i1	×
1119	coef_dr = mu_mom_diff * rhobar[i, j + 1, k0]	×
1120
1121	coef_dl = mu_mom_diff * rhobar[i, j, k0]	×
1122
1123	coef_d = 0.5 * (coef_dr + coef_dl)	×
1124
1125	u_diff =	×
1126	jac[i, j, k] *
1127	compute_momentum_diffusion_terms(state, i, j, k, V(), Z())
1128
1129	u_diff_f =	×
1130	jac[i, j + 1, k] *
1131	compute_momentum_diffusion_terms(state, i, j + 1, k, V(), Z())
1132
1133	u_diff_u =	×
1134	jac[i, j, k + 1] *
1135	compute_momentum_diffusion_terms(state, i, j, k + 1, V(), Z())
1136
1137	u_diff_fu =	×
1138	jac[i, j + 1, k + 1] *
1139	compute_momentum_diffusion_terms(state, i, j + 1, k + 1, V(), Z())
1140
1141	hrhov_diff =	×
1142	coef_d *
1143	0.5 *
1144	(
1145	jac[i, j, k] * jac[i, j, k + 1] * (u_diff + u_diff_u) /
1146	(jac[i, j, k] + jac[i, j, k + 1]) +
1147	jac[i, j + 1, k] *
1148	jac[i, j + 1, k + 1] *
1149	(u_diff_f + u_diff_fu) /
1150	(jac[i, j + 1, k] + jac[i, j + 1, k + 1])
1151	)
1152
1153	phiv[i, j, k, 3] -= hrhov_diff	×
1154	end	×
1155
1156	return	×
1157	end
1158
1159	function compute_fluxes!(	1,044✔
1160	state::State,
1161	old_predictands::Predictands,
1162	variable::W,
1163	)
1164	(; grid) = state	1,044✔
1165	(; re, uref, lref) = state.constants	1,044✔
1166	(; i0, i1, j0, j1, k0, k1) = state.domain	1,044✔
1167	(; jac, met) = grid	1,044✔
1168	(; pbar, rhobar) = state.atmosphere	1,044✔
1169	(; wtilde) = state.variables.reconstructions	1,044✔
1170	(; phiw) = state.variables.fluxes	1,044✔
1171	(; kinematic_diffusivity) = state.namelists.atmosphere	1,044✔
1172
1173	(u0, v0, w0) = (old_predictands.u, old_predictands.v, old_predictands.w)	1,044✔
1174
1175	#-----------------------------------------
1176	# Zonal fluxes
1177	#-----------------------------------------
1178
1179	@ivy for k in (k0 - 1):k1, j in j0:j1, i in (i0 - 1):i1	11,484✔
1180	wr = wtilde[i + 1, j, k, 1, 1]	779,724✔
1181	wl = wtilde[i, j, k, 1, 2]	779,724✔
1182
1183	pedger =	779,724✔
1184	0.5 * (
1185	jac[i, j, k] * pbar[i, j, k] +
1186	jac[i + 1, j, k] * pbar[i + 1, j, k]
1187	)
1188	puedger =	779,724✔
1189	0.5 * (
1190	jac[i, j, k + 1] * pbar[i, j, k + 1] +
1191	jac[i + 1, j, k + 1] * pbar[i + 1, j, k + 1]
1192	)
1193	usurf =	779,724✔
1194	(
1195	(jac[i, j, k + 1] + jac[i + 1, j, k + 1]) *
1196	pedger *
1197	u0[i, j, k] +
1198	(jac[i, j, k] + jac[i + 1, j, k]) * puedger * u0[i, j, k + 1]
1199	) / (
1200	jac[i, j, k] +
1201	jac[i + 1, j, k] +
1202	jac[i, j, k + 1] +
1203	jac[i + 1, j, k + 1]
1204	)
1205
1206	frhow = compute_flux(usurf, wl, wr)	779,724✔
1207
1208	phiw[i, j, k, 1] = frhow	779,724✔
1209	end	×
1210
1211	#-----------------------------------------
1212	# Meridional fluxes
1213	#-----------------------------------------
1214
1215	@ivy for k in (k0 - 1):k1, j in (j0 - 1):j1, i in i0:i1	11,484✔
1216	wf = wtilde[i, j + 1, k, 2, 1]	823,680✔
1217	wb = wtilde[i, j, k, 2, 2]	823,680✔
1218
1219	pedgef =	823,680✔
1220	0.5 * (
1221	jac[i, j, k] * pbar[i, j, k] +
1222	jac[i, j + 1, k] * pbar[i, j + 1, k]
1223	)
1224	puedgef =	823,680✔
1225	0.5 * (
1226	jac[i, j, k + 1] * pbar[i, j, k + 1] +
1227	jac[i, j + 1, k + 1] * pbar[i, j + 1, k + 1]
1228	)
1229	vsurf =	823,680✔
1230	(
1231	(jac[i, j, k + 1] + jac[i, j + 1, k + 1]) *
1232	pedgef *
1233	v0[i, j, k] +
1234	(jac[i, j, k] + jac[i, j + 1, k]) * puedgef * v0[i, j, k + 1]
1235	) / (
1236	jac[i, j, k] +
1237	jac[i, j + 1, k] +
1238	jac[i, j, k + 1] +
1239	jac[i, j + 1, k + 1]
1240	)
1241
1242	grhow = compute_flux(vsurf, wb, wf)	1,289,529✔
1243
1244	phiw[i, j, k, 2] = grhow	823,680✔
1245	end	×
1246
1247	#-----------------------------------------
1248	# Vertical fluxes
1249	#-----------------------------------------
1250
1251	@ivy for k in (k0 - 2):k1, j in j0:j1, i in i0:i1	12,528✔
1252	wu = wtilde[i, j, k + 1, 3, 1]	773,280✔
1253	wd = wtilde[i, j, k, 3, 2]	773,280✔
1254
1255	pedgeu =	773,280✔
1256	jac[i, j, k] *
1257	jac[i, j, k + 1] *
1258	(pbar[i, j, k] + pbar[i, j, k + 1]) /
1259	(jac[i, j, k] + jac[i, j, k + 1])
1260	puedgeu =	773,280✔
1261	jac[i, j, k + 1] *
1262	jac[i, j, k + 2] *
1263	(pbar[i, j, k + 1] + pbar[i, j, k + 2]) /
1264	(jac[i, j, k + 1] + jac[i, j, k + 2])
1265	wsurf = 0.5 * (pedgeu * w0[i, j, k] + puedgeu * w0[i, j, k + 1])	773,280✔
1266
1267	hrhow = compute_flux(wsurf, wd, wu)	1,171,536✔
1268
1269	phiw[i, j, k, 3] = hrhow	773,280✔
1270	end	×
1271
1272	#-------------------------------------------------------------------
1273	# Viscous fluxes
1274	#-------------------------------------------------------------------
1275
1276	if 1 / re <= eps() && kinematic_diffusivity == 0.0	1,044✔
1277	return	×
1278	end
1279
1280	#-----------------------------------------
1281	# Zonal fluxes
1282	#-----------------------------------------
1283
1284	@ivy for k in (k0 - 1):k1, j in j0:j1, i in (i0 - 1):i1	11,484✔
1285	coef_v = 1 / re * 0.5 * (rhobar[i, j, k0] + rhobar[i + 1, j, k0])	779,724✔
1286
1287	frhow_visc =	779,724✔
1288	coef_v *
1289	0.5 *
1290	(
1291	jac[i, j, k] *
1292	jac[i, j, k + 1] *
1293	(
1294	compute_stress_tensor(i, j, k, 3, 1, state) +
1295	compute_stress_tensor(i, j, k + 1, 3, 1, state)
1296	) / (jac[i, j, k] + jac[i, j, k + 1]) +
1297	jac[i + 1, j, k] *
1298	jac[i + 1, j, k + 1] *
1299	(
1300	compute_stress_tensor(i + 1, j, k, 3, 1, state) +
1301	compute_stress_tensor(i + 1, j, k + 1, 3, 1, state)
1302	) / (jac[i + 1, j, k] + jac[i + 1, j, k + 1])
1303	)
1304
1305	phiw[i, j, k, 1] -= frhow_visc	779,724✔
1306	end	×
1307
1308	#-----------------------------------------
1309	# Meridional fluxes
1310	#-----------------------------------------
1311
1312	@ivy for k in (k0 - 1):k1, j in (j0 - 1):j1, i in i0:i1	11,484✔
1313	coef_v = 1 / re * 0.5 * (rhobar[i, j, k0] + rhobar[i, j + 1, k0])	823,680✔
1314
1315	grhow_visc =	823,680✔
1316	coef_v *
1317	0.5 *
1318	(
1319	jac[i, j, k] *
1320	jac[i, j, k + 1] *
1321	(
1322	compute_stress_tensor(i, j, k, 3, 1, state) +
1323	compute_stress_tensor(i, j, k + 1, 3, 1, state)
1324	) / (jac[i, j, k] + jac[i, j, k + 1]) +
1325	jac[i, j + 1, k] *
1326	jac[i, j + 1, k + 1] *
1327	(
1328	compute_stress_tensor(i, j + 1, k, 3, 1, state) +
1329	compute_stress_tensor(i, j + 1, k + 1, 3, 1, state)
1330	) / (jac[i, j + 1, k] + jac[i, j + 1, k + 1])
1331	)
1332
1333	phiw[i, j, k, 2] -= grhow_visc	823,680✔
1334	end	×
1335
1336	#-----------------------------------------
1337	# Vertical fluxes
1338	#-----------------------------------------
1339
1340	@ivy for k in (k0 - 2):k1, j in j0:j1, i in i0:i1	12,528✔
1341	coef_v = 1 / re * rhobar[i, j, k0]	773,280✔
1342
1343	hrhow_visc =	773,280✔
1344	coef_v * (
1345	jac[i, j, k + 1] *
1346	met[i, j, k + 1, 1, 3] *
1347	compute_stress_tensor(i, j, k + 1, 3, 1, state) +
1348	jac[i, j, k + 1] *
1349	met[i, j, k + 1, 2, 3] *
1350	compute_stress_tensor(i, j, k + 1, 3, 2, state) +
1351	compute_stress_tensor(i, j, k + 1, 3, 3, state)
1352	)
1353
1354	phiw[i, j, k, 3] -= hrhow_visc	773,280✔
1355	end	×
1356
1357	#-------------------------------------------------------------------
1358	# Diffusion fluxes
1359	#-------------------------------------------------------------------
1360
1361	if kinematic_diffusivity == 0.0	1,044✔
1362	return	1,044✔
1363	end
1364
1365	mu_mom_diff = kinematic_diffusivity / uref / lref	×
1366
1367	#-----------------------------------------
1368	# Zonal fluxes
1369	#-----------------------------------------
1370
1371	@ivy for k in (k0 - 1):k1, j in j0:j1, i in (i0 - 1):i1	×
1372	coef_dr = mu_mom_diff * rhobar[i + 1, j, k0]	×
1373
1374	coef_dl = mu_mom_diff * rhobar[i, j, k0]	×
1375
1376	coef_d = 0.5 * (coef_dr + coef_dl)	×
1377
1378	w_diff =	×
1379	jac[i, j, k] *
1380	compute_momentum_diffusion_terms(state, i, j, k, W(), X())
1381
1382	w_diff_r =	×
1383	jac[i + 1, j, k] *
1384	compute_momentum_diffusion_terms(state, i + 1, j, k, W(), X())
1385
1386	w_diff_u =	×
1387	jac[i, j, k + 1] *
1388	compute_momentum_diffusion_terms(state, i, j, k + 1, W(), X())
1389
1390	w_diff_ru =	×
1391	jac[i + 1, j, k + 1] *
1392	compute_momentum_diffusion_terms(state, i + 1, j, k + 1, W(), X())
1393
1394	frhow_diff =	×
1395	coef_d *
1396	0.5 *
1397	(
1398	jac[i, j, k] * jac[i, j, k + 1] * (w_diff + w_diff_u) /
1399	(jac[i, j, k] + jac[i, j, k + 1]) +
1400	jac[i + 1, j, k] *
1401	jac[i + 1, j, k + 1] *
1402	(w_diff_r + w_diff_ru) /
1403	(jac[i + 1, j, k] + jac[i + 1, j, k + 1])
1404	)
1405
1406	phiw[i, j, k, 1] -= frhow_diff	×
1407	end	×
1408
1409	#-----------------------------------------
1410	# Meridional fluxes
1411	#-----------------------------------------
1412
1413	@ivy for k in (k0 - 1):k1, j in (j0 - 1):j1, i in i0:i1	×
1414	coef_dr = mu_mom_diff * rhobar[i, j + 1, k0]	×
1415
1416	coef_dl = mu_mom_diff * rhobar[i, j, k0]	×
1417
1418	coef_d = 0.5 * (coef_dr + coef_dl)	×
1419
1420	w_diff =	×
1421	jac[i, j, k] *
1422	compute_momentum_diffusion_terms(state, i, j, k, W(), Y())
1423
1424	w_diff_f =	×
1425	jac[i, j + 1, k] *
1426	compute_momentum_diffusion_terms(state, i, j + 1, k, W(), Y())
1427
1428	w_diff_u =	×
1429	jac[i, j, k + 1] *
1430	compute_momentum_diffusion_terms(state, i, j, k + 1, W(), Y())
1431
1432	w_diff_fu =	×
1433	jac[i, j + 1, k + 1] *
1434	compute_momentum_diffusion_terms(state, i, j + 1, k + 1, W(), Y())
1435
1436	grhow_diff =	×
1437	coef_d *
1438	0.5 *
1439	(
1440	jac[i, j, k] * jac[i, j, k + 1] * (w_diff + w_diff_u) /
1441	(jac[i, j, k] + jac[i, j, k + 1]) +
1442	jac[i, j + 1, k] *
1443	jac[i, j + 1, k + 1] *
1444	(w_diff_f + w_diff_fu) /
1445	(jac[i, j + 1, k] + jac[i, j + 1, k + 1])
1446	)
1447
1448	phiw[i, j, k, 2] -= grhow_diff	×
1449	end	×
1450
1451	#-----------------------------------------
1452	# Vertical fluxes
1453	#-----------------------------------------
1454
1455	@ivy for k in (k0 - 2):k1, j in j0:j1, i in i0:i1	×
1456	coef_d = mu_mom_diff * rhobar[i, j, k0]	×
1457
1458	hrhow_visc =	×
1459	coef_d *
1460	jac[i, j, k + 1] *
1461	compute_momentum_diffusion_terms(state, i, j, k + 1, W(), Z())
1462
1463	phiw[i, j, k, 3] -= hrhow_visc	×
1464	end	×
1465
1466	return	×
1467	end
1468
1469	function compute_fluxes!(	1,044✔
1470	state::State,
1471	predictands::Predictands,
1472	tracer_setup::NoTracer,
1473	)
1474	return	1,044✔
1475	end
1476
1477	function compute_fluxes!(	×
1478	state::State,
1479	predictands::Predictands,
1480	tracer_setup::TracerOn,
1481	)
1482	(; i0, i1, j0, j1, k0, k1) = state.domain	×
1483	(; jac) = state.grid	×
1484	(; pbar) = state.atmosphere	×
1485	(; tracerreconstructions, tracerfluxes) = state.tracer	×
1486
1487	(u0, v0, w0) = (predictands.u, predictands.v, predictands.w)	×
1488
1489	@ivy for field in 1:fieldcount(TracerPredictands)	×
1490	chir = getfield(tracerreconstructions, field)[2:end, :, :, 1, 1]	×
1491	chil = getfield(tracerreconstructions, field)[:, :, :, 1, 2]	×
1492	fchi = getfield(tracerfluxes, field)[:, :, :, 1]	×
1493	for k in k0:k1, j in j0:j1, i in (i0 - 1):i1	×
1494	pedger =	×
1495	0.5 * (
1496	jac[i, j, k] * pbar[i, j, k] +
1497	jac[i + 1, j, k] * pbar[i + 1, j, k]
1498	)
1499	usurf = pedger * u0[i, j, k]	×
1500
1501	fchi[i, j, k] = compute_flux(usurf, chil[i, j, k], chir[i, j, k])	×
1502	end	×
1503
1504	chif = getfield(tracerreconstructions, field)[:, 2:end, :, 2, 1]	×
1505	chib = getfield(tracerreconstructions, field)[:, :, :, 2, 2]	×
1506	gchi = getfield(tracerfluxes, field)[:, :, :, 2]	×
1507	for k in k0:k1, j in (j0 - 1):j1, i in i0:i1	×
1508	pedgef =	×
1509	0.5 * (
1510	jac[i, j, k] * pbar[i, j, k] +
1511	jac[i, j + 1, k] * pbar[i, j + 1, k]
1512	)
1513	vsurf = pedgef * v0[i, j, k]	×
1514
1515	gchi[i, j, k] = compute_flux(vsurf, chib[i, j, k], chif[i, j, k])	×
1516	end	×
1517
1518	chiu = getfield(tracerreconstructions, field)[:, :, 2:end, 3, 1]	×
1519	chid = getfield(tracerreconstructions, field)[:, :, :, 3, 2]	×
1520	hchi = getfield(tracerfluxes, field)[:, :, :, 3]	×
1521	for k in (k0 - 1):k1, j in j0:j1, i in i0:i1	×
1522	pedgeu =	×
1523	jac[i, j, k] *
1524	jac[i, j, k + 1] *
1525	(pbar[i, j, k] + pbar[i, j, k + 1]) /
1526	(jac[i, j, k] + jac[i, j, k + 1])
1527	wsurf = pedgeu * w0[i, j, k]	×
1528
1529	hchi[i, j, k] = compute_flux(wsurf, chid[i, j, k], chiu[i, j, k])	×
1530	end	×
1531	end	×
1532
1533	return	×
1534	end
1535
1536	function compute_fluxes!(	174✔
1537	state::State,
1538	predictands::Predictands,
1539	variable::Theta,
1540	)
1541	(; i0, i1, j0, j1, k0, k1) = state.domain	174✔
1542	(; jac, dx, dy, dz, met) = state.grid	174✔
1543	(; pbar, rhobar) = state.atmosphere	174✔
1544	(; phitheta) = state.variables.fluxes	174✔
1545	(; thermal_conductivity) = state.namelists.atmosphere	174✔
1546	(; uref, lref) = state.constants	174✔
1547	(; rho) = predictands	174✔
1548
1549	if thermal_conductivity == 0.0	174✔
1550	return	×
1551	end
1552
1553	mu_conduct = thermal_conductivity / uref / lref	174✔
1554
1555	#-----------------------------------------
1556	# Zonal fluxes
1557	#-----------------------------------------
1558
1559	@ivy for k in k0:k1, j in j0:j1, i in (i0 - 1):i1	1,740✔
1560	coef_t =	118,140✔
1561	mu_conduct *
1562	0.5 *
1563	(
1564	rhobar[i, j, k0] / rhobar[i, j, k] +
1565	rhobar[i + 1, j, k0] / rhobar[i + 1, j, k]
1566	)
1567
1568	thetal = pbar[i, j, k] / (rho[i, j, k] + rhobar[i, j, k])	118,140✔
1569	thetar = pbar[i + 1, j, k] / (rho[i + 1, j, k] + rhobar[i + 1, j, k])	118,140✔
1570
1571	thetad =	118,140✔
1572	0.5 * (
1573	pbar[i, j, k - 1] / (rho[i, j, k - 1] + rhobar[i, j, k - 1]) +
1574	pbar[i + 1, j, k - 1] /
1575	(rho[i + 1, j, k - 1] + rhobar[i + 1, j, k - 1])
1576	)
1577	thetau =	118,140✔
1578	0.5 * (
1579	pbar[i, j, k + 1] / (rho[i, j, k + 1] + rhobar[i, j, k + 1]) +
1580	pbar[i + 1, j, k + 1] /
1581	(rho[i + 1, j, k + 1] + rhobar[i + 1, j, k + 1])
1582	)
1583
1584	dtht_dxi =	118,140✔
1585	0.5 * (jac[i, j, k] + jac[i + 1, j, k]) * (thetar - thetal) / dx +
1586	0.5 *
1587	(
1588	jac[i, j, k] * met[i, j, k, 1, 3] +
1589	jac[i + 1, j, k] * met[i + 1, j, k, 1, 3]
1590	) *
1591	(thetau - thetad) / (2.0 * dz)
1592
1593	phitheta[i, j, k, 1] = -coef_t * dtht_dxi	118,140✔
1594	end	×
1595
1596	#-----------------------------------------
1597	# Meridional fluxes
1598	#-----------------------------------------
1599
1600	@ivy for k in k0:k1, j in (j0 - 1):j1, i in i0:i1	1,740✔
1601	coef_t =	124,800✔
1602	mu_conduct *
1603	0.5 *
1604	(
1605	rhobar[i, j, k0] / rhobar[i, j, k] +
1606	rhobar[i, j + 1, k0] / rhobar[i, j + 1, k]
1607	)
1608
1609	thetab = pbar[i, j, k] / (rho[i, j, k] + rhobar[i, j, k])	124,800✔
1610	thetaf = pbar[i, j + 1, k] / (rho[i, j + 1, k] + rhobar[i, j + 1, k])	124,800✔
1611
1612	thetad =	124,800✔
1613	0.5 * (
1614	pbar[i, j, k - 1] / (rho[i, j, k - 1] + rhobar[i, j, k - 1]) +
1615	pbar[i, j + 1, k - 1] /
1616	(rho[i, j + 1, k - 1] + rhobar[i, j + 1, k - 1])
1617	)
1618	thetau =	124,800✔
1619	0.5 * (
1620	pbar[i, j, k + 1] / (rho[i, j, k + 1] + rhobar[i, j, k + 1]) +
1621	pbar[i, j + 1, k + 1] /
1622	(rho[i, j + 1, k + 1] + rhobar[i, j + 1, k + 1])
1623	)
1624
1625	dtht_dyi =	124,800✔
1626	0.5 * (jac[i, j, k] + jac[i, j + 1, k]) * (thetaf - thetab) / dy +
1627	0.5 *
1628	(
1629	jac[i, j, k] * met[i, j, k, 2, 3] +
1630	jac[i, j + 1, k] * met[i, j + 1, k, 2, 3]
1631	) *
1632	(thetau - thetad) / (2.0 * dz)
1633
1634	phitheta[i, j, k, 2] = -coef_t * dtht_dyi	124,800✔
1635	end	×
1636
1637	#-----------------------------------------
1638	# Vertical fluxes
1639	#-----------------------------------------
1640
1641	@ivy for k in (k0 - 1):k1, j in j0:j1, i in i0:i1	1,914✔
1642	coef_t =	118,140✔
1643	mu_conduct * (
1644	jac[i, j, k + 1] * rhobar[i, j, 1] / rhobar[i, j, k] +
1645	jac[i, j, k] * rhobar[i, j, 1] / rhobar[i, j, k + 1]
1646	) / (jac[i, j, k + 1] + jac[i, j, k])
1647
1648	thetal =	118,140✔
1649	(
1650	jac[i - 1, j, k + 1] * pbar[i - 1, j, k] /
1651	(rho[i - 1, j, k] + rhobar[i - 1, j, k]) +
1652	jac[i - 1, j, k] * pbar[i - 1, j, k + 1] /
1653	(rho[i - 1, j, k + 1] + rhobar[i - 1, j, k + 1])
1654	) / (jac[i - 1, j, k + 1] + jac[i - 1, j, k])
1655
1656	thetar =	118,140✔
1657	(
1658	jac[i + 1, j, k + 1] * pbar[i + 1, j, k] /
1659	(rho[i + 1, j, k] + rhobar[i + 1, j, k]) +
1660	jac[i + 1, j, k] * pbar[i + 1, j, k + 1] /
1661	(rho[i + 1, j, k + 1] + rhobar[i + 1, j, k + 1])
1662	) / (jac[i + 1, j, k + 1] + jac[i + 1, j, k])
1663
1664	thetab =	118,140✔
1665	(
1666	jac[i, j - 1, k + 1] * pbar[i, j - 1, k] /
1667	(rho[i, j - 1, k] + rhobar[i, j - 1, k]) +
1668	jac[i, j - 1, k] * pbar[i, j - 1, k + 1] /
1669	(rho[i, j - 1, k + 1] + rhobar[i, j - 1, k + 1])
1670	) / (jac[i, j - 1, k + 1] + jac[i, j - 1, k])
1671
1672	thetaf =	118,140✔
1673	(
1674	jac[i, j + 1, k + 1] * pbar[i, j + 1, k] /
1675	(rho[i, j + 1, k] + rhobar[i, j + 1, k]) +
1676	jac[i, j + 1, k] * pbar[i, j + 1, k + 1] /
1677	(rho[i, j + 1, k + 1] + rhobar[i, j + 1, k + 1])
1678	) / (jac[i, j + 1, k + 1] + jac[i, j + 1, k])
1679
1680	thetad = pbar[i, j, k] / (rho[i, j, k] + rhobar[i, j, k])	118,140✔
1681	thetau = pbar[i, j, k + 1] / (rho[i, j, k + 1] + rhobar[i, j, k + 1])	118,140✔
1682
1683	dtht_dzi =	118,140✔
1684	jac[i, j, k] *
1685	jac[i, j, k + 1] *
1686	(
1687	(met[i, j, k, 1, 3] + met[i, j, k + 1, 1, 3]) *
1688	(thetar - thetal) / (2.0 * dx) +
1689	(met[i, j, k, 2, 3] + met[i, j, k + 1, 2, 3]) *
1690	(thetaf - thetab) / (2.0 * dy) +
1691	(met[i, j, k, 3, 3] + met[i, j, k + 1, 3, 3]) *
1692	(thetau - thetad) / dz
1693	) / (jac[i, j, k] + jac[i, j, k + 1])
1694
1695	phitheta[i, j, k, 3] = -coef_t * dtht_dzi	118,140✔
1696	end	×
1697
1698	return	174✔
1699	end