REMORA
Regional Modeling of Oceans Refined Adaptively
Loading...
Searching...
No Matches
REMORA_FennelBiology.cpp
Go to the documentation of this file.
1#include <REMORA_FennelBiology.H>
2#include <REMORA_DataStruct.H>
3#include <REMORA_Constants.H>
4#include <REMORA_DateClock.H>
5#include <AMReX_FArrayBox.H>
6#include <AMReX_MFIter.H>
7#include <AMReX_ParmParse.H>
8#include <AMReX_Utility.H>
9#include <algorithm>
10#include <cmath>
11
12using namespace amrex;
13
14namespace {
15
16#ifdef REMORA_USE_BIOLOGY_DIAG
17/**
18 * Emit one Path B diagnostic record.
19 *
20 * Format is contractual and must stay byte-identical to fennel_tag_line in
21 * Source/Biology/Fortran/REMORA_fennel_roms.F, whose Fortran edit descriptor
22 * is ('FENNEL-FORT ',a10,' iter=',i3,' i=',i5,' j=',i5,' k=',i4,' ',a4,' ',
23 * ES26.17E2). Both paths print k in REMORA's 0-based convention and iter
24 * 1-based, so the two logs diff without any index map.
25 *
26 * Compiled only under REMORA_USE_BIOLOGY_DIAG (GNUmake USE_BIOLOGY_DIAG=TRUE,
27 * CMake REMORA_ENABLE_BIOLOGY_DIAG=ON). This is called from inside the device
28 * lambda, so in a default build the printf and its device-side runtime are
29 * absent rather than merely unreached.
30 */
31AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
32void remora_fennel_tag (const char* tag, int iter, int i, int j, int k,
33 const char* var, Real val) noexcept
34{
35 printf("FENNEL-CPP %-10s iter=%3d i=%5d j=%5d k=%4d %-4s %26.17E\n",
36 tag, iter, i, j, k, var, double(val));
37}
38#endif
39
40AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
41Real fennel_pco2_water (Real T, Real S, Real TIC, Real TAlk) noexcept
42{
43 constexpr Real zero = Real(0.0);
44 constexpr int IbrackMax = 30;
45
46 // Determine coefficients for surface carbon chemistry.
47 const Real Tk = T + Real(273.15);
48 const Real centiTk = Real(0.01) * Tk;
49 const Real invTk = Real(1.0) / Tk;
50 const Real logTk = std::log(Tk);
51 const Real sqrtS = std::sqrt(S);
52 const Real SO4 = Real(19.924) * S / (Real(1000.0) - Real(1.005) * S);
53 const Real sqrtSO4 = std::sqrt(SO4);
54 const Real scl = S / Real(1.80655);
55
56 const Real alk = TAlk * Real(0.000001);
57 const Real dic = TIC * Real(0.000001);
58 const Real phos = zero;
59 const Real sili = zero;
60
61 // Correction term for non-ideality, ff=k0*(1-pH2O). Equation 13 with
62 // table 6 values from Weiss and Price (1980, Mar. Chem., 8, 347-359).
63 const Real ff = std::exp(-Real(162.8301) +
64 Real(218.2968) / centiTk +
65 std::log(centiTk) * Real(90.9241) -
66 centiTk * centiTk * Real(1.47696) +
67 S * (Real(0.025695) -
68 centiTk * (Real(0.025225) -
69 centiTk * Real(0.0049867))));
70
71 // Compute first (K1) and second (K2) dissociation constant of carboinic
72 // acid:
73 // K1 = [H][HCO3]/[H2CO3]
74 // K2 = [H][CO3]/[HCO3]
75 // From Millero (1995; page 664) using Mehrbach et al. (1973) data on
76 // seawater scale.
77 const Real K1 = std::pow(Real(10.0), Real(62.008) -
78 invTk * Real(3670.7) -
79 logTk * Real(9.7944) +
80 S * (Real(0.0118) - S * Real(0.000116)));
81 const Real K2 = std::pow(Real(10.0), -Real(4.777) -
82 invTk * Real(1394.7) +
83 S * (Real(0.0184) - S * Real(0.000118)));
84
85 // Compute dissociation constant of boric acid, Kb=[H][BO2]/[HBO2].
86 // From Millero (1995; page 669) using data from Dickson (1990).
87 const Real Kb = std::exp(-invTk * (Real(8966.90) +
88 sqrtS * (Real(2890.53) +
89 sqrtS * (Real(77.942) -
90 sqrtS * (Real(1.728) -
91 sqrtS * Real(0.0996))))) -
92 logTk * (Real(24.4344) +
93 sqrtS * (Real(25.085) + sqrtS * Real(0.2474))) +
94 Tk * (sqrtS * Real(0.053105)) +
95 Real(148.0248) +
96 sqrtS * (Real(137.1942) + sqrtS * Real(1.62142)));
97
98 // Compute first (K1p), second (K2p), and third (K3p) dissociation
99 // constant of phosphoric acid:
100 // K1p = [H][H2PO4]/[H3PO4]
101 // K2p = [H][HPO4]/[H2PO4]
102 // K3p = [H][PO4]/[HPO4]
103 // From DOE (1994) equations 7.2.20, 7.2.23, and 7.2.26, respectively.
104 // With footnote using data from Millero (1974).
105 const Real K1p = std::exp(Real(115.525) -
106 invTk * Real(4576.752) -
107 logTk * Real(18.453) +
108 sqrtS * (Real(0.69171) - invTk * Real(106.736)) -
109 S * (Real(0.01844) + invTk * Real(0.65643)));
110 const Real K2p = std::exp(Real(172.0883) -
111 invTk * Real(8814.715) -
112 logTk * Real(27.927) +
113 sqrtS * (Real(1.3566) - invTk * Real(160.340)) -
114 S * (Real(0.05778) - invTk * Real(0.37335)));
115 const Real K3p = std::exp(-Real(18.141) -
116 invTk * Real(3070.75) +
117 sqrtS * (Real(2.81197) + invTk * Real(17.27039)) -
118 S * (Real(0.09984) + invTk * Real(44.99486)));
119
120 // Compute dissociation constant of silica, Ksi=[H][SiO(OH)3]/[Si(OH)4].
121 // From Millero (1995; page 671) using data from Yao and Millero (1995).
122 const Real Ksi = std::exp(Real(117.385) -
123 invTk * Real(8904.2) -
124 logTk * Real(19.334) +
125 sqrtSO4 * (Real(3.5913) - invTk * Real(458.79)) -
126 SO4 * (Real(1.5998) - invTk * Real(188.74) -
127 SO4 * (Real(0.07871) - invTk * Real(12.1652))) +
128 std::log(Real(1.0) - Real(0.001005) * S));
129
130 // Compute ion product of whater, Kw = [H][OH].
131 // From Millero (1995; page 670) using composite data.
132 const Real Kw = std::exp(Real(148.9652) -
133 invTk * Real(13847.26) -
134 logTk * Real(23.6521) -
135 sqrtS * (Real(5.977) - invTk * Real(118.67) -
136 logTk * Real(1.0495)) -
137 S * Real(0.01615));
138
139 // Compute salinity constant of hydrogen sulfate, Ks = [H][SO4]/[HSO4].
140 // From Dickson (1990, J. chem. Thermodynamics 22, 113).
141 const Real Ks = std::exp(Real(141.328) -
142 invTk * Real(4276.1) -
143 logTk * Real(23.093) +
144 sqrtSO4 * (Real(324.57) - invTk * Real(13856.0) -
145 logTk * Real(47.986) - SO4 * invTk * Real(2698.0)) -
146 SO4 * (Real(771.54) - invTk * Real(35474.0) -
147 logTk * Real(114.723) - SO4 * invTk * Real(1776.0)) +
148 std::log(Real(1.0) - Real(0.001005) * S));
149
150 // Compute stability constant of hydrogen fluorid, Kf = [H][F]/[HF].
151 // From Dickson and Riley (1979) -- change pH scale to total.
152 const Real Kf = std::exp(-Real(12.641) +
153 invTk * Real(1590.2) +
154 sqrtSO4 * Real(1.525) +
155 std::log(Real(1.0) - Real(0.001005) * S) +
156 std::log(Real(1.0) + Real(0.1400) * scl / (Real(96.062) * Ks)));
157
158 // Calculate concentrations for borate (Uppstrom, 1974), sulfate (Morris
159 // and Riley, 1966), and fluoride (Riley, 1965).
160 const Real borate = Real(0.000232) * scl / Real(10.811);
161 const Real sulfate = Real(0.14) * scl / Real(96.062);
162 const Real fluoride = Real(0.000067) * scl / Real(18.9984);
163
164 // Bracket and bisection method.
165 Real X_lo = std::pow(Real(10.0), -Real(10.0));
166 Real X_hi = std::pow(Real(10.0), -Real(5.0));
167 Real X_mid = Real(0.5) * (X_lo + X_hi);
168 Real X = X_mid;
169 const Real K12 = K1 * K2;
170 const Real K12p = K1p * K2p;
171 const Real K123p = K12p * K3p;
172 const Real invKb = Real(1.0) / Kb;
173 const Real invKs = Real(1.0) / Ks;
174 const Real invKsi = Real(1.0) / Ksi;
175 Real fni1 = zero;
176 Real fni3 = zero;
177
178 for (int Ibrack = 0; Ibrack < IbrackMax; ++Ibrack) {
179 for (int Hstep = 1; Hstep <= 3; ++Hstep) {
180 if (Hstep == 1) { X = X_hi; }
181 if (Hstep == 2) { X = X_lo; }
182 if (Hstep == 3) { X = X_mid; }
183
184 // Set some common combinations of parameters used in the iterative
185 // [H+] solver.
186 const Real X2 = X * X;
187 const Real X3 = X2 * X;
188 const Real invX = Real(1.0) / X;
189 const Real A = X * (K12p + X * (K1p + X)) + K123p;
190 const Real B = X * (K1 + X) + K12;
191 const Real C = Real(1.0) / (Real(1.0) + sulfate * invKs);
192
193 // Evaluate f([H+]) for bracketing and mid-value cases.
194 const Real fni = dic * (K1 * X + Real(2.0) * K12) / B +
195 borate / (Real(1.0) + X * invKb) +
196 Kw * invX +
197 phos * (K12p * X + Real(2.0) * K123p - X3) / A +
198 sili / (Real(1.0) + X * invKsi) -
199 X * C -
200 sulfate / (Real(1.0) + Ks * invX * C) -
201 fluoride / (Real(1.0) + Kf * invX) - alk;
202 if (Hstep == 1) { fni1 = fni; }
203 if (Hstep == 3) { fni3 = fni; }
204 }
205
206 // Now, bracket solution within two of three.
207 if (fni3 == zero) {
208 break;
209 }
210 const Real ftest = fni1 / fni3;
211 if (ftest > zero) {
212 X_hi = X_mid;
213 } else {
214 X_lo = X_mid;
215 }
216 X_mid = Real(0.5) * (X_lo + X_hi);
217 }
218
219 // Last iteration gives value.
220 X = X_mid;
221
222 // Determine pCO2. Total Hydrogen ion concentration, Htotal = [H+].
223 const Real Htotal = X;
224 const Real Htotal2 = Htotal * Htotal;
225
226 // Calculate [CO2*] (mole/m3) as defined in DOE Methods Handbook 1994
227 // Version 2, ORNL/CDIAC-74, Dickson and Goyet, Eds. (Chapter 2,
228 // page 10, Eq A.49).
229 const Real CO2star = dic * Htotal2 / (Htotal2 + K1 * Htotal + K1 * K2);
230
231 // Add output argument for storing pCO2surf.
232 return CO2star * Real(1000000.0) / ff;
233}
234
235/**
236 * Atmospheric pCO2 (ppmv) for the surface CO2 gas exchange.
237 *
238 * \p time_seconds is the current time on the model clock and \p time_ref the
239 * reference date and calendar it is measured against, both handed to
240 * remora_caldate. Both are double, as the calendar is. The two
241 * time-dependent forms are ROMS's PCO2AIR_DATA and PCO2AIR_SECULAR; neither
242 * varies in space, so this is evaluated once per call on the host and the
243 * result handed to the kernel.
244 */
245Real
246fennel_pco2_air (REMORABiology::PCO2AirType type, double time_ref, double time_seconds,
247 Real pco2air_constant) noexcept
248{
249 if (type == REMORABiology::PCO2AirType::constant) {
250 return pco2air_constant;
251 }
252
253 constexpr Real pi2 = Real(6.2831853071796);
254
255 // The calendar works in double: a modern day number is too large for a
256 // float to resolve within a day.
257 int year = 0;
258 double yday_d = 0.0;
259 remora_caldate(time_ref, time_seconds / 86400.0, year, yday_d);
260 const Real yday = static_cast<Real>(yday_d);
261
262 if (type == REMORABiology::PCO2AirType::data) {
263 // Annual climatology of Laurent et al. (2017).
264 return Real(380.464) + Real(9.321) *
265 std::sin(pi2 * yday / Real(365.25) + Real(1.068));
266 }
267
268 // Secular trend. ROMS names this pmonth but computes years since 1951 and
269 // multiplies by 12 in the linear term; keep both, so the coefficients are
270 // the published ones.
271 constexpr Real D0 = Real(282.6);
272 constexpr Real D1 = Real(0.125);
273 constexpr Real D2 = Real(-7.18);
274 constexpr Real D3 = Real(0.86);
275 constexpr Real D4 = Real(-0.99);
276 constexpr Real D5 = Real(0.28);
277 constexpr Real D6 = Real(-0.80);
278 constexpr Real D7 = Real(0.06);
279
280 const Real pmonth = Real(year) - Real(1951.0) + yday / Real(365.0);
281 return D0 + D1 * pmonth * Real(12.0) +
282 D2 * std::sin(pi2 * pmonth + D3) +
283 D4 * std::sin(pi2 * pmonth + D5) +
284 D6 * std::sin(pi2 * pmonth + D7);
285}
286
287}
288
289namespace REMORABiology {
290
291void
292FennelParameters::init_params (const std::string& remora_prefix)
293{
294 ParmParse pp(remora_prefix + ".fennel");
295
296 pp.queryAdd("BioIter", BioIter);
297 pp.queryAdd("AttSW", AttSW);
298 pp.queryAdd("AttChl", AttChl);
299 pp.queryAdd("PARfrac", PARfrac);
300 pp.queryAdd("Vp0", Vp0);
301 pp.queryAdd("I_thNH4", I_thNH4);
302 pp.queryAdd("D_p5NH4", D_p5NH4);
303 pp.queryAdd("NitriR", NitriR);
304 pp.queryAdd("K_NO3", K_NO3);
305 pp.queryAdd("K_NH4", K_NH4);
306 pp.queryAdd("K_PO4", K_PO4);
307 pp.queryAdd("K_Phy", K_Phy);
308 pp.queryAdd("Chl2C_m", Chl2C_m);
309 pp.queryAdd("ChlMin", ChlMin);
310 pp.queryAdd("PhyCN", PhyCN);
311 pp.queryAdd("R_P2N", R_P2N);
312 pp.queryAdd("PhyIP", PhyIP);
313 pp.queryAdd("PhyIS", PhyIS);
314 pp.queryAdd("PhyMin", PhyMin);
315 pp.queryAdd("PhyMR", PhyMR);
316 pp.queryAdd("ZooAE_N", ZooAE_N);
317 pp.queryAdd("ZooCN", ZooCN);
318 pp.queryAdd("ZooBM", ZooBM);
319 pp.queryAdd("ZooER", ZooER);
320 pp.queryAdd("ZooGR", ZooGR);
321 pp.queryAdd("ZooMin", ZooMin);
322 pp.queryAdd("ZooMR", ZooMR);
323 pp.queryAdd("LDeRRN", LDeRRN);
324 pp.queryAdd("LDeRRC", LDeRRC);
325 pp.queryAdd("CoagR", CoagR);
326 pp.queryAdd("SDeRRN", SDeRRN);
327 pp.queryAdd("SDeRRC", SDeRRC);
328 pp.queryAdd("RDeRRN", RDeRRN);
329 pp.queryAdd("RDeRRC", RDeRRC);
330 pp.queryAdd("wPhy", wPhy);
331 pp.queryAdd("wLDet", wLDet);
332 pp.queryAdd("wSDet", wSDet);
333 pp.queryAdd("pCO2air", pCO2air);
334 pp.queryAdd("po4", po4);
335 pp.queryAdd("carbon", carbon);
336 pp.queryAdd("oxygen", oxygen);
337 pp.queryAdd("odu", odu);
338 pp.queryAdd("denitrification", denitrification);
339 pp.queryAdd("bio_sediment", bio_sediment);
340 pp.queryAdd("river_don", river_don);
341 pp.queryAdd("talk_nonconserv", talk_nonconserv);
342
343 // Alkalinity only exists as a tracer under carbon, so every term this
344 // option adds would have nowhere to go. ROMS would not compile in this
345 // combination; say so rather than silently ignoring the request.
346 if (talk_nonconserv && !carbon) {
347 amrex::Abort("remora.fennel.talk_nonconserv requires remora.fennel.carbon: "
348 "alkalinity is a carbon-block tracer");
349 }
350
351 static std::string pco2air_type_string = "constant";
352 pp.queryAdd("pco2air_type", pco2air_type_string);
353 const std::string pco2air_type_ci = amrex::toLower(pco2air_type_string);
354 if (pco2air_type_ci == "constant") {
355 pco2air_type = PCO2AirType::constant;
356 } else if (pco2air_type_ci == "data") {
357 pco2air_type = PCO2AirType::data;
358 } else if (pco2air_type_ci == "secular") {
359 pco2air_type = PCO2AirType::secular;
360 } else {
361 amrex::Abort("Unknown remora.fennel.pco2air_type: " + pco2air_type_string +
362 ". Expected constant, data, or secular.");
363 }
364
365 static std::string co2_schmidt_string = "wanninkhof1992";
366 pp.queryAdd("co2_schmidt", co2_schmidt_string);
367 const std::string co2_schmidt_ci = amrex::toLower(co2_schmidt_string);
368 if (co2_schmidt_ci == "wanninkhof1992" || co2_schmidt_ci == "w92") {
369 co2_schmidt = CO2SchmidtType::wanninkhof1992;
370 } else if (co2_schmidt_ci == "wanninkhof2014" || co2_schmidt_ci == "rw14") {
371 co2_schmidt = CO2SchmidtType::wanninkhof2014;
372 } else {
373 amrex::Abort("Unknown remora.fennel.co2_schmidt: " + co2_schmidt_string +
374 ". Expected wanninkhof1992 or wanninkhof2014.");
375 }
376
377 static std::string o2_schmidt_string = "wanninkhof1992";
378 pp.queryAdd("oxygen_schmidt", o2_schmidt_string);
379 const std::string o2_schmidt_ci = amrex::toLower(o2_schmidt_string);
380 if (o2_schmidt_ci == "wanninkhof1992" || o2_schmidt_ci == "w92") {
381 o2_schmidt = O2SchmidtType::wanninkhof1992;
382 } else if (o2_schmidt_ci == "wanninkhof2014" || o2_schmidt_ci == "rw14") {
383 o2_schmidt = O2SchmidtType::wanninkhof2014;
384 } else if (o2_schmidt_ci == "ocmip") {
385 o2_schmidt = O2SchmidtType::ocmip;
386 } else {
387 amrex::Abort("Unknown remora.fennel.oxygen_schmidt: " + o2_schmidt_string +
388 ". Expected wanninkhof1992, wanninkhof2014, or ocmip.");
389 }
390}
391
392void
393FennelBiology::configure (const std::string& remora_prefix)
394{
395 // Keep the name remora_prefix: Exec/Generic/build_schema_remora.py
396 // resolves it to "remora" when generating inputs_generic, and any other
397 // name is misattributed.
398 m_parameters.init_params(remora_prefix);
399 ParmParse pp(remora_prefix);
400 pp.queryAdd("use_biology_cpp_answer", m_use_cpp);
401 pp.queryAdd("biology_debug", m_debug);
402 pp.queryAdd("biology_debug_i", m_debug_i);
403 pp.queryAdd("biology_debug_j", m_debug_j);
404#ifndef REMORA_USE_FENNEL_FORT
405 if (m_use_cpp == 0) {
406 amrex::Abort("remora.use_biology_cpp_answer = 0 selects the ROMS "
407 "Fennel Fortran bridge, which is not compiled in. "
408 "Rebuild with USE_FENNEL_FORT=TRUE (GNUmake) or "
409 "-DREMORA_ENABLE_FENNEL_FORT=ON (CMake).");
410 }
411#endif
412#ifndef REMORA_USE_BIOLOGY_DIAG
413 if (m_debug > 0) {
414 amrex::Abort("remora.biology_debug > 0 requests the Fennel parity "
415 "diagnostics, which are not compiled in. Rebuild with "
416 "USE_BIOLOGY_DIAG=TRUE (GNUmake) or "
417 "-DREMORA_ENABLE_BIOLOGY_DIAG=ON (CMake).");
418 }
419#endif
420 Vector<std::string> names = {"NO3", "NH4", "chlorophyll", "phytoplankton",
421 "zooplankton", "LdetritusN", "SdetritusN"};
422 if (m_parameters.river_don) {
423 names.emplace_back("RdetritusN");
424 }
425 if (m_parameters.po4) {
426 names.emplace_back("PO4");
427 }
428 if (m_parameters.carbon) {
429 names.emplace_back("LdetritusC");
430 names.emplace_back("SdetritusC");
431 names.emplace_back("TIC");
432 names.emplace_back("alkalinity");
433 if (m_parameters.river_don) {
434 names.emplace_back("RdetritusC");
435 }
436 }
437 if (m_parameters.oxygen) {
438 names.emplace_back("oxygen");
439 }
440 if (m_parameters.odu) {
441 names.emplace_back("ODU");
442 }
443 m_declaration.tracer_names = std::move(names);
444 m_declaration.plot_aliases = {"fennel"};
445}
446
447}
448
449void
450REMORABiology::FennelBiology::advance (const BiologyAdvanceContext& context) const
451{
452 if (m_parameters.BioIter <= 0) {
453 return;
454 }
455#ifdef REMORA_USE_FENNEL_FORT
456 if (m_use_cpp == 0) {
457 advance_fortran(context);
458 return;
459 }
460#endif
461 auto const& mf_cons_old = context.old_state;
462 auto& mf_cons_new = context.new_state;
463 auto const& solverChoice = context.solver_choice;
464 const int N = context.top_layer;
465 const int Bio_comp = context.component_start;
466 const Real dt_lev = context.dt;
467 const auto parms = m_parameters;
468
469 const bool need_wind = (parms.oxygen || parms.carbon) && solverChoice.bulk_fluxes;
470 const bool need_stress = (parms.oxygen || parms.carbon) && !solverChoice.bulk_fluxes;
471
472 if (context.hz == nullptr || context.z_w == nullptr || context.mask == nullptr ||
473 context.shortwave == nullptr ||
474 (need_stress && (context.stress_u == nullptr || context.stress_v == nullptr)) ||
475 (need_wind && (context.wind_u == nullptr || context.wind_v == nullptr))) {
476 amrex::Abort("Fennel biology requires Hz, z_w, rmask, srflx, and, when carbon or "
477 "oxygen is active, either uwind/vwind (bulk_fluxes) or sustr/svstr");
478 }
479
480 // Set time-stepping according to the number of iterations.
481 const Real dtdays = dt_lev / Real(86400.0) / static_cast<Real>(parms.BioIter);
482 const Real rho0 = solverChoice.rho0;
483 const bool use_po4 = parms.po4;
484 const bool use_carbon = parms.carbon;
485 const bool use_oxygen = parms.oxygen;
486 const bool use_odu = parms.odu;
487 const bool use_denitrification = parms.denitrification;
488 const bool use_river_don = parms.river_don;
489 const bool use_river_don_c = parms.river_don && parms.carbon;
490 const bool use_talk_nonconserv = parms.talk_nonconserv;
491 const bool use_salt = use_oxygen || use_carbon;
492 const bool do_bulk_flux = solverChoice.bulk_fluxes;
493
494 // Set biological tracer component identifiers. The biology block starts after the
495 // passive scalars, so this is Tracer_comp only when the run carries no dye.
496 const auto bio_comp = REMORABiology::Fennel::components(parms, Bio_comp);
497
498 // Schmidt-number coefficients and the leading rate coefficient for the gas
499 // transfer velocity. Each pair belongs to one published relation, so they
500 // are selected together rather than independently; ROMS makes the same
501 // pairings with RW14_OXYGEN_SC, OCMIP_OXYGEN_SC and RW14_CO2_SC.
502 Real A_O2 = Real(1953.4); // Schmidt number
503 Real B_O2 = Real(128.0); // coefficients from
504 Real C_O2 = Real(3.9918); // Wanninkhof (1992)
505 Real D_O2 = Real(0.050091);
506 Real E_O2 = Real(0.0);
507 Real o2_rate = Real(0.31);
508 if (parms.o2_schmidt == REMORABiology::O2SchmidtType::wanninkhof2014) {
509 A_O2 = Real(1920.4);
510 B_O2 = Real(135.6);
511 C_O2 = Real(5.2122);
512 D_O2 = Real(0.10939);
513 E_O2 = Real(0.00093777);
514 o2_rate = Real(0.251);
515 } else if (parms.o2_schmidt == REMORABiology::O2SchmidtType::ocmip) {
516 // Keeling et al. (1998); Sc is slightly smaller up to about 35C. ROMS
517 // pairs this set with the 1992 rate coefficient, not the 2014 one.
518 A_O2 = Real(1638.0);
519 B_O2 = Real(81.83);
520 C_O2 = Real(1.483);
521 D_O2 = Real(0.008004);
522 E_O2 = Real(0.0);
523 }
524
525 Real A_CO2 = Real(2073.1); // Schmidt number
526 Real B_CO2 = Real(125.62); // coefficients from
527 Real C_CO2 = Real(3.6276); // Wanninkhof (1992)
528 Real D_CO2 = Real(0.043219);
529 Real E_CO2 = Real(0.0);
530 Real co2_rate = Real(0.31);
531 if (parms.co2_schmidt == REMORABiology::CO2SchmidtType::wanninkhof2014) {
532 A_CO2 = Real(2116.8);
533 B_CO2 = Real(136.25);
534 C_CO2 = Real(4.7353);
535 D_CO2 = Real(0.092307);
536 E_CO2 = Real(0.0007555);
537 co2_rate = Real(0.251);
538 }
539
540 // Atmospheric pCO2 varies in time but not in space, so evaluate it once
541 // here rather than once per column. t_old is the time at the start of the
542 // step, which is the state the kernel reads.
543 const Real pco2air = fennel_pco2_air(parms.pco2air_type, solverChoice.time_ref,
544 context.model_time, parms.pCO2air);
545
546 // A non-positive atmospheric pCO2 reverses the sign of the air-sea flux for the
547 // whole run. Only the surface CO2 exchange reads pco2air, so a run without carbon
548 // is unaffected by its value and must not be stopped on account of it.
549 if (use_carbon && pco2air <= zero) {
550 if (parms.pco2air_type == REMORABiology::PCO2AirType::constant) {
551 amrex::Abort("remora.fennel.pCO2air must be a positive partial pressure; got "
552 + std::to_string(pco2air) + " ppmv.");
553 } else {
554 // The secular fit is anchored to 1951 and extrapolates to a negative partial
555 // pressure well before it, so a run whose clock sits near the origin of its
556 // calendar silently draws CO2 out of the ocean for its whole length. That is
557 // the default with remora.time_ref = 0, whose epoch is 0001-01-01.
558 amrex::Abort("remora.fennel.pco2air_type gives a non-positive atmospheric pCO2 ("
559 + std::to_string(pco2air) + " ppmv) at this model time. The secular"
560 " trend is fitted around 1951, so put the run in a real year: set"
561 " remora.time_ref to the reference date and remora.start_time to the"
562 " offset from it.");
563 }
564 }
565
566#ifdef REMORA_USE_BIOLOGY_DIAG
567 // Path B half of the frozen diagnostic contract. Tag names, field order
568 // and numeric format must stay byte-identical to the Fortran emitters in
569 // Source/Biology/Fortran/REMORA_fennel_roms.F; see tag_map.md.
570 const int dbg_level = m_debug;
571 const int dbg_i = m_debug_i;
572 const int dbg_j = m_debug_j;
573#endif
574
575 // Scratch arrays mirror the ROMS Bio, qc, bR, bL, WR, WL, and
576 // inverse-thickness work arrays, but use per-tile FArrayBox storage.
577 // Optional tracer scratch is only allocated for active runtime options.
578 int ncell = 0;
579 const int sc_no3 = ncell++;
580 const int sc_nh4 = ncell++;
581 const int sc_chlo = ncell++;
582 const int sc_phyt = ncell++;
583 const int sc_zoop = ncell++;
584 const int sc_lden = ncell++;
585 const int sc_sden = ncell++;
586 const int sc_rden = use_river_don ? ncell++ : -1;
587 const int sc_po4 = use_po4 ? ncell++ : -1;
588 const int sc_ldec = use_carbon ? ncell++ : -1;
589 const int sc_sdec = use_carbon ? ncell++ : -1;
590 const int sc_tic = use_carbon ? ncell++ : -1;
591 const int sc_talk = use_carbon ? ncell++ : -1;
592 const int sc_rdec = use_river_don_c ? ncell++ : -1;
593 const int sc_oxyg = use_oxygen ? ncell++ : -1;
594 const int sc_odu = use_odu ? ncell++ : -1;
595 const int sc_temp = ncell++;
596 const int sc_salt = use_salt ? ncell++ : -1;
597 const int sc_inv_hz = ncell++;
598 const int sc_inv_hz2 = ncell++;
599 const int sc_inv_hz3 = ncell++;
600 const int sc_qc = ncell++;
601 const int sc_bR = ncell++;
602 const int sc_bL = ncell++;
603 const int sc_WR = ncell++;
604 const int sc_WL = ncell++;
605
606 int nw = 0;
607 const int sw_FC = nw++;
608
609 for (MFIter mfi(mf_cons_new, TilingIfNotGPU()); mfi.isValid(); ++mfi) {
610 Box bx = mfi.tilebox();
611 Box bx2d = bx;
612 bx2d.makeSlab(2, 0);
613 Box wbx = surroundingNodes(bx, 2);
614
615 FArrayBox fab_cell(bx, ncell, amrex::The_Async_Arena());
616 fab_cell.template setVal<RunOn::Device>(zero);
617 FArrayBox fab_w(wbx, nw, amrex::The_Async_Arena());
618 fab_w.template setVal<RunOn::Device>(zero);
619
620 auto no3 = fab_cell.array(sc_no3);
621 auto nh4 = fab_cell.array(sc_nh4);
622 auto chlo = fab_cell.array(sc_chlo);
623 auto phyt = fab_cell.array(sc_phyt);
624 auto zoop = fab_cell.array(sc_zoop);
625 auto lden = fab_cell.array(sc_lden);
626 auto sden = fab_cell.array(sc_sden);
627 Array4<Real> rden;
628 if (use_river_don) {
629 rden = fab_cell.array(sc_rden);
630 }
631 Array4<Real> po4;
632 if (use_po4) {
633 po4 = fab_cell.array(sc_po4);
634 }
635 Array4<Real> ldec;
636 Array4<Real> sdec;
637 Array4<Real> tic;
638 Array4<Real> talk;
639 if (use_carbon) {
640 ldec = fab_cell.array(sc_ldec);
641 sdec = fab_cell.array(sc_sdec);
642 tic = fab_cell.array(sc_tic);
643 talk = fab_cell.array(sc_talk);
644 }
645 Array4<Real> rdec;
646 if (use_river_don_c) {
647 rdec = fab_cell.array(sc_rdec);
648 }
649 Array4<Real> oxyg;
650 if (use_oxygen) {
651 oxyg = fab_cell.array(sc_oxyg);
652 }
653 Array4<Real> odu;
654 if (use_odu) {
655 odu = fab_cell.array(sc_odu);
656 }
657 auto temp = fab_cell.array(sc_temp);
658 Array4<Real> salt;
659 if (use_salt) {
660 salt = fab_cell.array(sc_salt);
661 }
662 auto inv_hz = fab_cell.array(sc_inv_hz);
663 auto inv_hz2 = fab_cell.array(sc_inv_hz2);
664 auto inv_hz3 = fab_cell.array(sc_inv_hz3);
665 auto qc = fab_cell.array(sc_qc);
666 auto bR = fab_cell.array(sc_bR);
667 auto bL = fab_cell.array(sc_bL);
668 auto WR = fab_cell.array(sc_WR);
669 auto WL = fab_cell.array(sc_WL);
670 auto FC = fab_w.array(sw_FC);
671
672 Array4<Real const> const& state_old = mf_cons_old.const_array(mfi);
673 Array4<Real> const& state_new = mf_cons_new.array(mfi);
674 Array4<Real const> const& Hz = context.hz->const_array(mfi);
675 Array4<Real const> const& z_w = context.z_w->const_array(mfi);
676 Array4<Real const> const& srflx = context.shortwave->const_array(mfi);
677 Array4<Real const> const& mskr = context.mask->const_array(mfi);
678 // ROMS takes the gas-exchange transfer velocity from Uwind/Vwind
679 // under BULK_FLUXES and from the surface stress otherwise; only the
680 // selected pair is read.
681 Array4<Real const> sustr;
682 Array4<Real const> svstr;
683 Array4<Real const> uwind;
684 Array4<Real const> vwind;
685 if (use_oxygen || use_carbon) {
686 if (do_bulk_flux) {
687 uwind = context.wind_u->const_array(mfi);
688 vwind = context.wind_v->const_array(mfi);
689 } else {
690 sustr = context.stress_u->const_array(mfi);
691 svstr = context.stress_v->const_array(mfi);
692 }
693 }
694
695 ParallelFor(bx2d, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept
696 {
697 constexpr Real zero = Real(0.0);
698 // Land columns contribute nothing: the increment below is scaled by rmask.
699 // Skip them outright rather than evaluating the chemistry and multiplying the
700 // result away -- a masked column may hold a NetCDF fill value, and the logs
701 // and square roots in the pCO2 and oxygen-saturation blocks would turn that
702 // into a NaN that survives the rmask factor (NaN * 0 is NaN) and lands in
703 // cons_new. ROMS guards pCO2_water with #ifdef MASKING for the same reason.
704 if (mskr(i,j,0) == zero) return;
705
706 constexpr Real one = Real(1.0);
707 constexpr Real two = Real(2.0);
708 constexpr Real eps = Real(1.0e-20);
709 constexpr Real minval = Real(1.0e-6);
710 constexpr Real cff_weno = Real(1.0e-14);
711
712 constexpr Real OA0 = Real(2.00907); // Oxygen
713 constexpr Real OA1 = Real(3.22014); // saturation
714 constexpr Real OA2 = Real(4.05010); // coefficients
715 constexpr Real OA3 = Real(4.94457);
716 constexpr Real OA4 = Real(-0.256847);
717 constexpr Real OA5 = Real(3.88767);
718 constexpr Real OB0 = Real(-0.00624523);
719 constexpr Real OB1 = Real(-0.00737614);
720 constexpr Real OB2 = Real(-0.0103410);
721 constexpr Real OB3 = Real(-0.00817083);
722 constexpr Real OC0 = Real(-0.000000488682);
723 constexpr Real rOxNO3 = Real(8.625); // 138/16
724 constexpr Real rOxNH4 = Real(6.625); // 106/16
725 constexpr Real rOxNH4Denit = Real(115.0) / Real(16.0);
726 constexpr Real denitrification_NH4_fraction = Real(4.0) / Real(16.0);
727 constexpr Real l2mol = Real(1000.0) / Real(22.3916); // liter to mol
728
729 constexpr Real A1 = Real(-60.2409); // surface
730 constexpr Real A2 = Real(93.4517); // CO2
731 constexpr Real A3 = Real(23.3585); // solubility
732 constexpr Real B1 = Real(0.023517); // coefficients
733 constexpr Real B2 = Real(-0.023656);
734 constexpr Real B3 = Real(0.0047036);
735
736 // Compute inverse thickness to avoid repeated divisions.
737 //
738 // Extract biological variables from tracer arrays, place them
739 // into scratch arrays, and restrict their values to be positive
740 // definite. At input, tracers are read from the nstp state
741 // (cons_old); the final increment is applied to nnew
742 // (cons_new) below.
743 for (int k = 0; k <= N; ++k) {
744 inv_hz(i,j,k) = one / Hz(i,j,k);
745 no3(i,j,k) = amrex::max(zero, state_old(i,j,k,bio_comp.no3));
746 nh4(i,j,k) = amrex::max(zero, state_old(i,j,k,bio_comp.nh4));
747 chlo(i,j,k) = amrex::max(zero, state_old(i,j,k,bio_comp.chlo));
748 phyt(i,j,k) = amrex::max(zero, state_old(i,j,k,bio_comp.phyt));
749 zoop(i,j,k) = amrex::max(zero, state_old(i,j,k,bio_comp.zoop));
750 lden(i,j,k) = amrex::max(zero, state_old(i,j,k,bio_comp.lden));
751 sden(i,j,k) = amrex::max(zero, state_old(i,j,k,bio_comp.sden));
752 if (use_river_don) {
753 rden(i,j,k) = amrex::max(zero, state_old(i,j,k,bio_comp.rden));
754 }
755 if (use_po4) {
756 po4(i,j,k) = amrex::max(zero, state_old(i,j,k,bio_comp.po4));
757 }
758 if (use_carbon) {
759 ldec(i,j,k) = amrex::max(zero, state_old(i,j,k,bio_comp.ldec));
760 sdec(i,j,k) = amrex::max(zero, state_old(i,j,k,bio_comp.sdec));
761 tic(i,j,k) = amrex::max(zero, state_old(i,j,k,bio_comp.tic));
762 tic(i,j,k) = amrex::min(tic(i,j,k), Real(3000.0));
763 tic(i,j,k) = amrex::max(tic(i,j,k), Real(400.0));
764 talk(i,j,k) = amrex::max(zero, state_old(i,j,k,bio_comp.talk));
765 if (use_river_don) {
766 rdec(i,j,k) = amrex::max(zero, state_old(i,j,k,bio_comp.rdec));
767 }
768 }
769 if (use_oxygen) {
770 oxyg(i,j,k) = amrex::max(zero, state_old(i,j,k,bio_comp.oxyg));
771 }
772 if (use_odu) {
773 odu(i,j,k) = amrex::max(zero, state_old(i,j,k,bio_comp.odu));
774 }
775 // Extract potential temperature and salinity.
776 temp(i,j,k) = amrex::min(state_old(i,j,k,Temp_comp), Real(35.0));
777 if (use_salt) {
778 salt(i,j,k) = amrex::max(state_old(i,j,k,Salt_comp), zero);
779 }
780 }
781 for (int k = 0; k < N; ++k) {
782 inv_hz2(i,j,k) = one / (Hz(i,j,k) + Hz(i,j,k+1));
783 }
784 for (int k = 1; k < N; ++k) {
785 inv_hz3(i,j,k) = one / (Hz(i,j,k-1) + Hz(i,j,k) + Hz(i,j,k+1));
786 }
787
788 // Calculate surface Photosynthetically Available Radiation
789 // (PAR). REMORA stores srflx in Watts/m2, so only PARfrac is
790 // applied here.
791 const Real PARsur = parms.PARfrac * srflx(i,j,0);
792
793 // Emit one Tier-1 tag for this column: every level, every active
794 // tracer, in the frozen field order. Inactive tracers are skipped
795 // on both paths so the two logs stay line-aligned.
796 //
797 // Two definitions rather than an #ifdef'd body: the disabled form
798 // captures nothing, so the closure and its ~14 Array4 captures
799 // disappear from the device lambda entirely instead of surviving
800 // as dead weight. Call sites are identical either way.
801#ifndef REMORA_USE_BIOLOGY_DIAG
802 auto tag_state = [] (const char*, int) noexcept {};
803#else
804 auto tag_state = [=] (const char* tag, int it) noexcept
805 {
806 if (dbg_level <= 0) return;
807 if (dbg_level == 1 && (i != dbg_i || j != dbg_j)) return;
808 for (int k = 0; k <= N; ++k) {
809 remora_fennel_tag(tag, it, i, j, k, "NO3 ", no3(i,j,k));
810 remora_fennel_tag(tag, it, i, j, k, "NH4 ", nh4(i,j,k));
811 remora_fennel_tag(tag, it, i, j, k, "CHLO", chlo(i,j,k));
812 remora_fennel_tag(tag, it, i, j, k, "PHYT", phyt(i,j,k));
813 remora_fennel_tag(tag, it, i, j, k, "ZOOP", zoop(i,j,k));
814 remora_fennel_tag(tag, it, i, j, k, "LDEN", lden(i,j,k));
815 remora_fennel_tag(tag, it, i, j, k, "SDEN", sden(i,j,k));
816 if (use_river_don) {
817 remora_fennel_tag(tag, it, i, j, k, "RDEN", rden(i,j,k));
818 }
819 if (use_po4) {
820 remora_fennel_tag(tag, it, i, j, k, "PO4 ", po4(i,j,k));
821 }
822 if (use_carbon) {
823 remora_fennel_tag(tag, it, i, j, k, "LDEC", ldec(i,j,k));
824 remora_fennel_tag(tag, it, i, j, k, "SDEC", sdec(i,j,k));
825 remora_fennel_tag(tag, it, i, j, k, "TIC ", tic(i,j,k));
826 remora_fennel_tag(tag, it, i, j, k, "TALK", talk(i,j,k));
827 if (use_river_don) {
828 remora_fennel_tag(tag, it, i, j, k, "RDEC", rdec(i,j,k));
829 }
830 }
831 if (use_oxygen) {
832 remora_fennel_tag(tag, it, i, j, k, "OXYG", oxyg(i,j,k));
833 }
834 if (use_odu) {
835 remora_fennel_tag(tag, it, i, j, k, "ODU ", odu(i,j,k));
836 }
837 }
838 };
839#endif
840
841 tag_state("G00_PRE ", 0);
842
843 // Start internal iterations to achieve convergence of the
844 // nonlinear backward-implicit solution.
845 //
846 // During the iterative procedure a series of fractional time
847 // steps are performed in a chained mode, splitting different
848 // biological conversion processes in sequence of the food chain.
849 // In all stages the concentration of the component being consumed
850 // is treated implicitly, so the algorithm guarantees
851 // non-negative values. The overall algorithm, as well as any
852 // stage of it, is formulated in conservative form except for
853 // explicit sinking.
854 //
855 // The iterative loop is to iterate toward a backward-Euler
856 // treatment of all terms; it damps numerical oscillations but
857 // does not improve formal accuracy.
858 for (int iter = 0; iter < parms.BioIter; ++iter) {
859
860 // Light-limited computations.
861 //
862 // Compute attenuation coefficient based on chlorophyll-a in
863 // each grid box. Then attenuate surface PAR down into the
864 // water column. Thus, PAR at a depth depends on the whole
865 // distribution of chlorophyll-a above.
866 if (PARsur > zero) {
867 Real PAR = PARsur;
868 for (int k = N; k >= 0; --k) {
869 // Compute average light attenuation for each grid
870 // cell. Other attenuation contributions like
871 // suspended sediment or CDOM can be added here.
872 const Real thickness = z_w(i,j,k+1) - z_w(i,j,k);
873 const Real Att = (parms.AttSW + parms.AttChl * chlo(i,j,k)) * thickness;
874 const Real ExpAtt = std::exp(-Att);
875 const Real Itop = PAR;
876 const Real PARavg = Itop * (one - ExpAtt) / Att;
877
878 // Compute Chlorophyll-a phytoplankton ratio,
879 // [mg Chla / mg C].
880 const Real cff = parms.PhyCN * Real(12.0);
881 const Real Chl2C = amrex::min(chlo(i,j,k) / (phyt(i,j,k) * cff + eps), parms.Chl2C_m);
882
883 // Temperature-limited and light-limited growth rate
884 // (Eppley, R.W., 1972, Fishery Bulletin,
885 // 70: 1063-1085; here 0.59=ln(2)*0.851).
886 const Real Vp = parms.Vp0 * Real(0.59) * std::pow(Real(1.066), temp(i,j,k));
887 const Real fac1_light = PARavg * parms.PhyIS;
888 const Real Epp = Vp / std::sqrt(Vp * Vp + fac1_light * fac1_light);
889 const Real t_PPmax = Epp * fac1_light;
890
891 // Nutrient-limitation terms (Parker 1993 Ecol
892 // Mod., 66, 113-120).
893 const Real cff1 = nh4(i,j,k) * parms.K_NH4;
894 const Real cff2 = no3(i,j,k) * parms.K_NO3;
895 const Real inhNH4 = one / (one + cff1);
896 const Real L_NH4 = cff1 / (one + cff1);
897 const Real L_NO3 = cff2 * inhNH4 / (one + cff2);
898 const Real LTOT = L_NO3 + L_NH4;
899 Real LMIN = LTOT;
900 if (use_po4) {
901 const Real cff3 = po4(i,j,k) * parms.K_PO4;
902 const Real L_PO4 = cff3 / (one + cff3);
903 LMIN = amrex::min(LTOT, L_PO4);
904 }
905
906 // Nitrate and ammonium uptake by Phytoplankton.
907 Real cff4 = zero;
908 Real cff5 = zero;
909 Real cff6 = zero;
910 if (use_po4) {
911 const Real mu = dtdays * t_PPmax * LMIN;
912 cff4 = mu * phyt(i,j,k) * L_NO3 /
913 amrex::max(minval, LTOT) / amrex::max(minval, no3(i,j,k));
914 cff5 = mu * phyt(i,j,k) * L_NH4 /
915 amrex::max(minval, LTOT) / amrex::max(minval, nh4(i,j,k));
916 cff6 = parms.R_P2N * mu * phyt(i,j,k) /
917 amrex::max(minval, po4(i,j,k));
918 } else {
919 const Real fac1 = dtdays * t_PPmax;
920 cff4 = fac1 * parms.K_NO3 * inhNH4 / (one + cff2) * phyt(i,j,k);
921 cff5 = fac1 * parms.K_NH4 / (one + cff1) * phyt(i,j,k);
922 }
923 no3(i,j,k) = no3(i,j,k) / (one + cff4);
924 nh4(i,j,k) = nh4(i,j,k) / (one + cff5);
925 if (use_po4) {
926 po4(i,j,k) = po4(i,j,k) / (one + cff6);
927 }
928 const Real N_Flux_NewProd = no3(i,j,k) * cff4;
929 const Real N_Flux_RegProd = nh4(i,j,k) * cff5;
930 phyt(i,j,k) += N_Flux_NewProd + N_Flux_RegProd;
931 chlo(i,j,k) += (dtdays * t_PPmax * t_PPmax * LMIN * LMIN *
932 parms.Chl2C_m * chlo(i,j,k)) /
933 (parms.PhyIS * amrex::max(Chl2C, eps) * PARavg + eps);
934 if (use_oxygen) {
935 oxyg(i,j,k) += N_Flux_NewProd * rOxNO3 + N_Flux_RegProd * rOxNH4;
936 }
937 if (use_carbon) {
938 // Total inorganic carbon (CO2) uptake during
939 // phytoplankton growth.
940 const Real cff1_carbon = parms.PhyCN *
941 (N_Flux_NewProd + N_Flux_RegProd);
942 tic(i,j,k) -= cff1_carbon;
943 if (use_talk_nonconserv) {
944 // Account for the uptake of NO3 on total alkalinity.
945 talk(i,j,k) += N_Flux_NewProd - N_Flux_RegProd;
946 }
947 }
948
949 // The Nitrification of NH4 ==> NO3 is thought to
950 // occur only in dark and only in aerobic water
951 // (Olson, R. J., 1981, JMR: 39, 227-238).
952 //
953 // NH4+ + 3/2 O2 ==> NO2- + H2O, via Nitrosomonas
954 // bacteria; NO2- + 1/2 O2 ==> NO3-, via Nitrobacter
955 // bacteria.
956 //
957 // Note that the entire process has a total loss of
958 // two moles of O2 per mole of NH4. If we were to
959 // resolve NO2 profiles, this is where we would change
960 // the code to split out the differential effects of
961 // the two different bacteria types. If OXYGEN is
962 // defined, nitrification is inhibited at low oxygen
963 // concentrations using a Michaelis-Menten term.
964 Real nitri_fac = dtdays * parms.NitriR;
965 if (use_oxygen) {
966 const Real fac2 = amrex::max(oxyg(i,j,k), zero); // O2 max
967 const Real fac3 = amrex::max(fac2 / (Real(3.0) + fac2), zero); // MM for O2 dependence
968 nitri_fac *= fac3;
969 }
970 const Real light_ratio = (PARavg - parms.I_thNH4) /
971 (parms.D_p5NH4 + PARavg - two * parms.I_thNH4);
972 const Real light_inhibit = one - amrex::max(zero, light_ratio);
973 const Real nitri = nitri_fac * light_inhibit;
974 nh4(i,j,k) = nh4(i,j,k) / (one + nitri);
975 const Real N_Flux_Nitrifi = nh4(i,j,k) * nitri;
976 no3(i,j,k) += N_Flux_Nitrifi;
977 if (use_oxygen) {
978 oxyg(i,j,k) -= two * N_Flux_Nitrifi;
979 }
980 if (use_talk_nonconserv) {
981 talk(i,j,k) -= two * N_Flux_Nitrifi;
982 }
983
984 // Light attenuation at the bottom of the grid cell.
985 // It is the starting PAR value for the next deeper
986 // vertical grid cell.
987 PAR = Itop * ExpAtt;
988 }
989 } else {
990 // If PARsur=0, nitrification occurs at the maximum rate
991 // (NitriR).
992 const Real nitri = dtdays * parms.NitriR;
993 for (int k = N; k >= 0; --k) {
994 nh4(i,j,k) = nh4(i,j,k) / (one + nitri);
995 const Real N_Flux_Nitrifi = nh4(i,j,k) * nitri;
996 no3(i,j,k) += N_Flux_Nitrifi;
997 if (use_oxygen) {
998 oxyg(i,j,k) -= two * N_Flux_Nitrifi;
999 }
1000 if (use_talk_nonconserv) {
1001 talk(i,j,k) -= two * N_Flux_Nitrifi;
1002 }
1003 }
1004 }
1005
1006 tag_state("G01_LIGHT ", iter + 1);
1007
1008 // Phytoplankton grazing by zooplankton (rate: ZooGR),
1009 // phytoplankton assimilated to zooplankton (fraction:
1010 // ZooAE_N) and egested to small detritus, and
1011 // phytoplankton mortality (rate: PhyMR) to small detritus
1012 // [Landry 1993 L and O 38:468-472].
1013 const Real grazing_fac = dtdays * parms.ZooGR;
1014 const Real phyt_mort_fac = dtdays * parms.PhyMR;
1015 for (int k = 0; k <= N; ++k) {
1016 // Phytoplankton grazing by zooplankton.
1017 const Real phy2 = phyt(i,j,k) * phyt(i,j,k);
1018 const Real graze = grazing_fac * zoop(i,j,k) * phyt(i,j,k) /
1019 (parms.K_Phy + phy2);
1020 const Real consume = one / (one + graze);
1021 phyt(i,j,k) *= consume;
1022 chlo(i,j,k) *= consume;
1023 // Phytoplankton assimilated to zooplankton and
1024 // egested to small detritus.
1025 const Real N_Flux_Assim = graze * phyt(i,j,k) * parms.ZooAE_N;
1026 const Real N_Flux_Egest = phyt(i,j,k) * graze * (one - parms.ZooAE_N);
1027 zoop(i,j,k) += N_Flux_Assim;
1028 sden(i,j,k) += N_Flux_Egest;
1029
1030 // Phytoplankton mortality, limited by a
1031 // phytoplankton minimum.
1032 const Real N_Flux_Pmortal = phyt_mort_fac * amrex::max(phyt(i,j,k) - parms.PhyMin, zero);
1033 phyt(i,j,k) -= N_Flux_Pmortal;
1034 chlo(i,j,k) -= phyt_mort_fac * amrex::max(chlo(i,j,k) - parms.ChlMin, zero);
1035 sden(i,j,k) += N_Flux_Pmortal;
1036 if (use_carbon) {
1037 sdec(i,j,k) += parms.PhyCN * (N_Flux_Egest + N_Flux_Pmortal) +
1038 (parms.PhyCN - parms.ZooCN) * N_Flux_Assim;
1039 }
1040 }
1041
1042 tag_state("G02_GRAZE ", iter + 1);
1043
1044 // Zooplankton basal metabolism to NH4 (rate: ZooBM),
1045 // zooplankton mortality to small detritus (rate: ZooMR),
1046 // and zooplankton ingestion-related excretion (rate: ZooER).
1047 const Real zoo_metab_fac = dtdays * parms.ZooBM;
1048 const Real zoo_mort_fac = dtdays * parms.ZooMR;
1049 const Real zoo_excrete_fac = dtdays * parms.ZooER;
1050 for (int k = 0; k <= N; ++k) {
1051 const Real phy2 = phyt(i,j,k) * phyt(i,j,k);
1052 const Real ingestion_excretion = zoo_excrete_fac * phy2 /
1053 (parms.K_Phy + phy2);
1054 const Real cff2 = zoo_mort_fac * zoop(i,j,k);
1055 const Real cff3 = ingestion_excretion * parms.ZooAE_N;
1056 zoop(i,j,k) = zoop(i,j,k) / (one + cff2 + cff3);
1057 // Zooplankton mortality and excretion.
1058 const Real N_Flux_Zmortal = cff2 * zoop(i,j,k);
1059 const Real N_Flux_Zexcret = cff3 * zoop(i,j,k);
1060 nh4(i,j,k) += N_Flux_Zexcret;
1061 if (use_po4) {
1062 po4(i,j,k) += parms.R_P2N * N_Flux_Zexcret;
1063 }
1064 sden(i,j,k) += N_Flux_Zmortal;
1065
1066 // Zooplankton basal metabolism, limited by a
1067 // zooplankton minimum.
1068 const Real N_Flux_Zmetabo = zoo_metab_fac * amrex::max(zoop(i,j,k) - parms.ZooMin, zero);
1069 zoop(i,j,k) -= N_Flux_Zmetabo;
1070 nh4(i,j,k) += N_Flux_Zmetabo;
1071 if (use_po4) {
1072 po4(i,j,k) += parms.R_P2N * N_Flux_Zmetabo;
1073 }
1074 if (use_oxygen) {
1075 oxyg(i,j,k) -= rOxNH4 * (N_Flux_Zmetabo + N_Flux_Zexcret);
1076 }
1077 if (use_carbon) {
1078 sdec(i,j,k) += parms.ZooCN * N_Flux_Zmortal;
1079 tic(i,j,k) += parms.ZooCN * (N_Flux_Zmetabo + N_Flux_Zexcret);
1080 if (use_talk_nonconserv) {
1081 talk(i,j,k) += N_Flux_Zmetabo + N_Flux_Zexcret;
1082 }
1083 }
1084 }
1085
1086 tag_state("G03_ZMETAB", iter + 1);
1087
1088 // Coagulation of phytoplankton and small detritus to
1089 // large detritus.
1090 const Real coag_fac = dtdays * parms.CoagR;
1091 for (int k = 0; k <= N; ++k) {
1092 const Real coag = coag_fac * (sden(i,j,k) + phyt(i,j,k));
1093 const Real consume = one / (one + coag);
1094 phyt(i,j,k) *= consume;
1095 chlo(i,j,k) *= consume;
1096 sden(i,j,k) *= consume;
1097 const Real N_Flux_CoagP = phyt(i,j,k) * coag;
1098 const Real N_Flux_CoagD = sden(i,j,k) * coag;
1099 lden(i,j,k) += N_Flux_CoagP + N_Flux_CoagD;
1100 if (use_carbon) {
1101 sdec(i,j,k) -= parms.PhyCN * N_Flux_CoagD;
1102 ldec(i,j,k) += parms.PhyCN * (N_Flux_CoagP + N_Flux_CoagD);
1103 }
1104 }
1105
1106 tag_state("G04_COAG ", iter + 1);
1107
1108 // Detritus recycling to NH4, remineralization.
1109 if (use_oxygen) {
1110 for (int k = 0; k <= N; ++k) {
1111 const Real fac1 = amrex::max(oxyg(i,j,k) - Real(6.0), zero); // O2 off max
1112 const Real fac2 = amrex::max(fac1 / (Real(3.0) + fac1), zero); // MM for O2 dependence
1113 const Real remin_s = dtdays * parms.SDeRRN * fac2;
1114 const Real remin_s_consume = one / (one + remin_s);
1115 const Real remin_l = dtdays * parms.LDeRRN * fac2;
1116 const Real remin_l_consume = one / (one + remin_l);
1117 sden(i,j,k) *= remin_s_consume;
1118 lden(i,j,k) *= remin_l_consume;
1119 const Real N_Flux_RemineS = sden(i,j,k) * remin_s;
1120 const Real N_Flux_RemineL = lden(i,j,k) * remin_l;
1121 nh4(i,j,k) += N_Flux_RemineS + N_Flux_RemineL;
1122 if (use_po4) {
1123 po4(i,j,k) += parms.R_P2N * (N_Flux_RemineS + N_Flux_RemineL);
1124 }
1125 oxyg(i,j,k) -= (N_Flux_RemineS + N_Flux_RemineL) * rOxNH4;
1126 if (use_talk_nonconserv) {
1127 talk(i,j,k) += N_Flux_RemineS + N_Flux_RemineL;
1128 }
1129 if (use_river_don) {
1130 const Real remin_r = dtdays * parms.RDeRRN * fac2;
1131 const Real remin_r_consume = one / (one + remin_r);
1132 rden(i,j,k) *= remin_r_consume;
1133 const Real N_Flux_RemineR = rden(i,j,k) * remin_r;
1134 nh4(i,j,k) += N_Flux_RemineR;
1135 if (use_po4) {
1136 po4(i,j,k) += parms.R_P2N * N_Flux_RemineR;
1137 }
1138 oxyg(i,j,k) -= N_Flux_RemineR * rOxNH4;
1139 if (use_talk_nonconserv) {
1140 talk(i,j,k) += N_Flux_RemineR;
1141 }
1142 }
1143 }
1144 } else {
1145 const Real remin_s = dtdays * parms.SDeRRN;
1146 const Real remin_l = dtdays * parms.LDeRRN;
1147 const Real remin_s_consume = one / (one + remin_s);
1148 const Real remin_l_consume = one / (one + remin_l);
1149 const Real remin_r = dtdays * parms.RDeRRN;
1150 const Real remin_r_consume = one / (one + remin_r);
1151 for (int k = 0; k <= N; ++k) {
1152 sden(i,j,k) *= remin_s_consume;
1153 lden(i,j,k) *= remin_l_consume;
1154 const Real N_Flux_RemineS = sden(i,j,k) * remin_s;
1155 const Real N_Flux_RemineL = lden(i,j,k) * remin_l;
1156 nh4(i,j,k) += N_Flux_RemineS + N_Flux_RemineL;
1157 if (use_po4) {
1158 po4(i,j,k) += parms.R_P2N * (N_Flux_RemineS + N_Flux_RemineL);
1159 }
1160 if (use_talk_nonconserv) {
1161 talk(i,j,k) += N_Flux_RemineS + N_Flux_RemineL;
1162 }
1163 if (use_river_don) {
1164 rden(i,j,k) *= remin_r_consume;
1165 const Real N_Flux_RemineR = rden(i,j,k) * remin_r;
1166 nh4(i,j,k) += N_Flux_RemineR;
1167 if (use_po4) {
1168 po4(i,j,k) += parms.R_P2N * N_Flux_RemineR;
1169 }
1170 if (use_talk_nonconserv) {
1171 talk(i,j,k) += N_Flux_RemineR;
1172 }
1173 }
1174 }
1175 }
1176
1177 tag_state("G05_REMINN", iter + 1);
1178
1179 if (use_oxygen) {
1180 // Surface O2 gas exchange.
1181 //
1182 // Compute surface O2 gas exchange.
1183 const Real cff1 = rho0 * Real(550.0);
1184 const Real cff2 = dtdays * o2_rate * Real(24.0) / Real(100.0);
1185 const int k = N;
1186
1187 Real u10squ;
1188 // Compute O2 transfer velocity: u10squared (u10 in m/s).
1189 if (do_bulk_flux) {
1190 u10squ = uwind(i,j,0) * uwind(i,j,0) + vwind(i,j,0) * vwind(i,j,0);
1191 } else {
1192 u10squ = cff1 *
1193 std::sqrt((Real(0.5) * (sustr(i,j,0) + sustr(i+1,j,0))) *
1194 (Real(0.5) * (sustr(i,j,0) + sustr(i+1,j,0))) +
1195 (Real(0.5) * (svstr(i,j,0) + svstr(i,j+1,0))) *
1196 (Real(0.5) * (svstr(i,j,0) + svstr(i,j+1,0))));
1197 }
1198 const Real SchmidtN_Ox = A_O2 - temp(i,j,k) *
1199 (B_O2 - temp(i,j,k) *
1200 (C_O2 - temp(i,j,k) *
1201 (D_O2 - temp(i,j,k) * E_O2)));
1202 const Real cff3 = cff2 * u10squ * std::sqrt(Real(660.0) / SchmidtN_Ox);
1203
1204 // Calculate O2 saturation concentration using Garcia and
1205 // Gordon L and O (1992) formula, (EXP(AA) is in ml/l).
1206 const Real TS = std::log((Real(298.15) - temp(i,j,k)) /
1207 (Real(273.15) + temp(i,j,k)));
1208 const Real AA = OA0 + TS * (OA1 + TS * (OA2 + TS * (OA3 + TS * (OA4 + TS * OA5)))) +
1209 salt(i,j,k) * (OB0 + TS * (OB1 + TS * (OB2 + TS * OB3))) +
1210 OC0 * salt(i,j,k) * salt(i,j,k);
1211
1212 // Convert from ml/l to mmol/m3.
1213 const Real O2satu = l2mol * std::exp(AA);
1214
1215 // Add in O2 gas exchange.
1216 const Real O2_Flux = cff3 * (O2satu - oxyg(i,j,k));
1217 oxyg(i,j,k) += O2_Flux * inv_hz(i,j,k);
1218
1219 tag_state("G06_O2FLX ", iter + 1);
1220 }
1221
1222 if (use_carbon) {
1223 // Allow different remineralization rates for detrital C
1224 // and detrital N.
1225 const Real remin_sc = dtdays * parms.SDeRRC;
1226 const Real remin_sc_consume = one / (one + remin_sc);
1227 const Real remin_lc = dtdays * parms.LDeRRC;
1228 const Real remin_lc_consume = one / (one + remin_lc);
1229 const Real remin_rc = dtdays * parms.RDeRRC;
1230 const Real remin_rc_consume = one / (one + remin_rc);
1231 for (int k = 0; k <= N; ++k) {
1232 sdec(i,j,k) *= remin_sc_consume;
1233 ldec(i,j,k) *= remin_lc_consume;
1234 const Real C_Flux_RemineS = sdec(i,j,k) * remin_sc;
1235 const Real C_Flux_RemineL = ldec(i,j,k) * remin_lc;
1236 tic(i,j,k) += C_Flux_RemineS + C_Flux_RemineL;
1237 if (use_river_don_c) {
1238 rdec(i,j,k) *= remin_rc_consume;
1239 tic(i,j,k) += rdec(i,j,k) * remin_rc;
1240 }
1241 }
1242
1243 if (!use_talk_nonconserv) {
1244 // Alkalinity is treated as a diagnostic variable. TAlk =
1245 // f(S[PSU]) following Brewer et al. (1986). Under
1246 // talk_nonconserv it is prognostic instead, and the
1247 // increments accumulated above are what carry it, so
1248 // overwriting here would discard every one of them.
1249 for (int k = 0; k <= N; ++k) {
1250 talk(i,j,k) = Real(587.05) + Real(50.56) * salt(i,j,k);
1251 }
1252 }
1253
1254 tag_state("G07_REMINC", iter + 1);
1255
1256 // Surface CO2 gas exchange.
1257 //
1258 // Compute equilibrium partial pressure inorganic carbon
1259 // (ppmv) at the surface.
1260 const int k = N;
1261 const Real pCO2 = fennel_pco2_water(temp(i,j,k), salt(i,j,k),
1262 tic(i,j,k), talk(i,j,k));
1263
1264 // Compute surface CO2 gas exchange.
1265 const Real cff1 = rho0 * Real(550.0);
1266 const Real cff2 = dtdays * co2_rate * Real(24.0) / Real(100.0);
1267
1268 // Compute CO2 transfer velocity: u10squared (u10 in m/s).
1269 Real u10squ;
1270 if (do_bulk_flux) {
1271 u10squ = uwind(i,j,0) * uwind(i,j,0) + vwind(i,j,0) * vwind(i,j,0);
1272 } else {
1273 u10squ = cff1 *
1274 std::sqrt((Real(0.5) * (sustr(i,j,0) + sustr(i+1,j,0))) *
1275 (Real(0.5) * (sustr(i,j,0) + sustr(i+1,j,0))) +
1276 (Real(0.5) * (svstr(i,j,0) + svstr(i,j+1,0))) *
1277 (Real(0.5) * (svstr(i,j,0) + svstr(i,j+1,0))));
1278 }
1279 const Real SchmidtN = A_CO2 - temp(i,j,k) *
1280 (B_CO2 - temp(i,j,k) *
1281 (C_CO2 - temp(i,j,k) *
1282 (D_CO2 - temp(i,j,k) * E_CO2)));
1283 const Real cff3 = cff2 * u10squ * std::sqrt(Real(660.0) / SchmidtN);
1284
1285 // Calculate CO2 solubility [mol/(kg.atm)] using Weiss
1286 // (1974) formula.
1287 const Real TempK = Real(0.01) * (temp(i,j,k) + Real(273.15));
1288 const Real CO2_sol = std::exp(A1 +
1289 A2 / TempK +
1290 A3 * std::log(TempK) +
1291 salt(i,j,k) * (B1 + TempK * (B2 + B3 * TempK)));
1292
1293 // Add in CO2 gas exchange.
1294 const Real CO2_Flux = cff3 * CO2_sol * (pco2air - pCO2);
1295 tic(i,j,k) += CO2_Flux * inv_hz(i,j,k);
1296
1297 tag_state("G08_CO2FLX", iter + 1);
1298 }
1299
1300 // Vertical sinking terms.
1301 //
1302 // Set vertical sinking identification and sinking velocity
1303 // vectors in the same order: phytoplankton, chlorophyll,
1304 // small nitrogen-detritus, large nitrogen-detritus, and,
1305 // when CARBON is active, small and large carbon-detritus.
1306 const int nsink = use_carbon ? 6 : 4;
1307 for (int isink = 0; isink < nsink; ++isink) {
1308 Array4<Real> bio = phyt;
1309 Real wbio = parms.wPhy;
1310 int sediment_type = 1;
1311 const bool phyt_sink = (isink == 0);
1312 if (isink == 1) {
1313 bio = chlo;
1314 wbio = parms.wPhy;
1315 sediment_type = 0;
1316 } else if (isink == 2) {
1317 bio = sden;
1318 wbio = parms.wSDet;
1319 } else if (isink == 3) {
1320 bio = lden;
1321 wbio = parms.wLDet;
1322 } else if (isink == 4) {
1323 bio = sdec;
1324 wbio = parms.wSDet;
1325 sediment_type = 2;
1326 } else if (isink == 5) {
1327 bio = ldec;
1328 wbio = parms.wLDet;
1329 sediment_type = 2;
1330 }
1331
1332 // Copy concentration of biological particulates into
1333 // scratch array qc, restricted to positive values.
1334 for (int k = 0; k <= N; ++k) {
1335 qc(i,j,k) = bio(i,j,k);
1336 bR(i,j,k) = qc(i,j,k);
1337 bL(i,j,k) = qc(i,j,k);
1338 WR(i,j,k) = zero;
1339 WL(i,j,k) = zero;
1340 }
1341 for (int k = 0; k <= N + 1; ++k) {
1342 FC(i,j,k) = zero;
1343 }
1344
1345 // Reconstruct vertical profiles as parabolic segments
1346 // within each grid box.
1347 for (int iface = 1; iface <= N; ++iface) {
1348 FC(i,j,iface) = (qc(i,j,iface) - qc(i,j,iface-1)) * inv_hz2(i,j,iface-1);
1349 }
1350 for (int k = 1; k < N; ++k) {
1351 Real dltR = Hz(i,j,k) * FC(i,j,k+1);
1352 Real dltL = Hz(i,j,k) * FC(i,j,k);
1353 const Real cff = Hz(i,j,k-1) + two * Hz(i,j,k) + Hz(i,j,k+1);
1354 const Real cffR = cff * FC(i,j,k+1);
1355 const Real cffL = cff * FC(i,j,k);
1356
1357 // Apply PPM monotonicity constraint to prevent
1358 // oscillations within the grid box.
1359 if ((dltR * dltL) <= zero) {
1360 dltR = zero;
1361 dltL = zero;
1362 } else if (std::abs(dltR) > std::abs(cffL)) {
1363 dltR = cffL;
1364 } else if (std::abs(dltL) > std::abs(cffR)) {
1365 dltL = cffR;
1366 }
1367
1368 // Compute right and left side values of parabolic
1369 // segments; WR and WL are measures of quadratic
1370 // variations.
1371 const Real curv = (dltR - dltL) * inv_hz3(i,j,k);
1372 dltR -= curv * Hz(i,j,k+1);
1373 dltL += curv * Hz(i,j,k-1);
1374 bR(i,j,k) = qc(i,j,k) + dltR;
1375 bL(i,j,k) = qc(i,j,k) - dltL;
1376 WR(i,j,k) = (two * dltR - dltL) * (two * dltR - dltL);
1377 WL(i,j,k) = (dltR - two * dltL) * (dltR - two * dltL);
1378 }
1379 // Reconcile interfacial values using a WENO procedure
1380 // so the whole profile remains monotonic.
1381 for (int k = 1; k < N - 1; ++k) {
1382 const Real dltL = amrex::max(cff_weno, WL(i,j,k));
1383 const Real dltR = amrex::max(cff_weno, WR(i,j,k+1));
1384 bR(i,j,k) = (dltR * bR(i,j,k) + dltL * bL(i,j,k+1)) / (dltR + dltL);
1385 bL(i,j,k+1) = bR(i,j,k);
1386 }
1387 FC(i,j,N + 1) = zero;
1388 bR(i,j,N) = qc(i,j,N);
1389 bL(i,j,N) = qc(i,j,N);
1390 bR(i,j,N - 1) = qc(i,j,N);
1391 bL(i,j,1) = qc(i,j,0);
1392 bR(i,j,0) = qc(i,j,0);
1393 bL(i,j,0) = qc(i,j,0);
1394
1395 // Apply monotonicity constraint again, since
1396 // reconciled interfacial values may cause non-monotonic
1397 // behavior inside the grid box.
1398 for (int k = 0; k <= N; ++k) {
1399 Real dltR = bR(i,j,k) - qc(i,j,k);
1400 Real dltL = qc(i,j,k) - bL(i,j,k);
1401 const Real cffR = two * dltR;
1402 const Real cffL = two * dltL;
1403 if ((dltR * dltL) < zero) {
1404 dltR = zero;
1405 dltL = zero;
1406 } else if (std::abs(dltR) > std::abs(cffL)) {
1407 dltR = cffL;
1408 } else if (std::abs(dltL) > std::abs(cffR)) {
1409 dltL = cffR;
1410 }
1411 bR(i,j,k) = qc(i,j,k) + dltR;
1412 bL(i,j,k) = qc(i,j,k) - dltL;
1413 }
1414
1415 // After reconstruction, compute vertical advective
1416 // fluxes. The algorithm is free of a CFL criterion by
1417 // allowing semi-Lagrangian integration bounds to use as
1418 // many upstream grid boxes as necessary.
1419 //
1420 // WL is the z-coordinate of the departure point for grid
1421 // box interface z_w. FC is the finite-volume flux.
1422 const Real sink_distance = dtdays * std::abs(wbio);
1423 for (int k = 0; k <= N; ++k) {
1424 FC(i,j,k) = zero;
1425 WL(i,j,k) = z_w(i,j,k) + sink_distance;
1426 WR(i,j,k) = Hz(i,j,k) * qc(i,j,k);
1427 }
1428 FC(i,j,N + 1) = zero;
1429 for (int k = 0; k <= N; ++k) {
1430 int ksource = k;
1431 Real flux = zero;
1432 for (int ks = k; ks < N; ++ks) {
1433 if (WL(i,j,k) > z_w(i,j,ks+1)) {
1434 ksource = ks + 1;
1435 flux += WR(i,j,ks);
1436 }
1437 }
1438 // Finalize flux computation by adding the
1439 // fractional part from the source grid box.
1440 const Real cu = amrex::min(one, (WL(i,j,k) - z_w(i,j,ksource)) *
1441 inv_hz(i,j,ksource));
1442 flux += Hz(i,j,ksource) * cu *
1443 (bL(i,j,ksource) + cu * (Real(0.5) * (bR(i,j,ksource) - bL(i,j,ksource)) -
1444 (Real(1.5) - cu) * (bR(i,j,ksource) + bL(i,j,ksource) -
1445 two * qc(i,j,ksource))));
1446 FC(i,j,k) = flux;
1447 }
1448 for (int k = 0; k <= N; ++k) {
1449 bio(i,j,k) = qc(i,j,k) +
1450 (FC(i,j,k+1) - FC(i,j,k)) * inv_hz(i,j,k);
1451 }
1452 // Particulate flux reaching the seafloor is
1453 // remineralized and returned to the dissolved nitrate
1454 // pool. Without this conversion, particulate material
1455 // falls out of the system. This is a temporary fix to
1456 // restore total nitrogen conservation. It will be replaced
1457 // later by a parameterization that includes the time delay
1458 // of remineralization and dissolved oxygen.
1459 if (parms.bio_sediment && sediment_type != 0) {
1460 const Real cff1 = FC(i,j,0) * inv_hz(i,j,0);
1461 if (sediment_type == 1) {
1462 if (use_denitrification) {
1463 nh4(i,j,0) += cff1 * denitrification_NH4_fraction;
1464 if (use_oxygen) {
1465 oxyg(i,j,0) -= cff1 * rOxNH4Denit;
1466 }
1467 } else {
1468 nh4(i,j,0) += cff1;
1469 if (use_oxygen) {
1470 oxyg(i,j,0) -= cff1 * rOxNH4;
1471 }
1472 if (use_talk_nonconserv) {
1473 talk(i,j,0) += cff1;
1474 }
1475 }
1476 if (use_po4) {
1477 po4(i,j,0) += cff1 * parms.R_P2N;
1478 }
1479 if (use_carbon && phyt_sink) {
1480 tic(i,j,0) += cff1 * parms.PhyCN;
1481 }
1482 } else if (use_carbon) {
1483 tic(i,j,0) += cff1;
1484 }
1485 }
1486 }
1487
1488 tag_state("G09_SINK ", iter + 1);
1489 }
1490
1491 // Emitted at the same semantic point as the Fortran path: right
1492 // after ITER_LOOP and before anything else touches the scratch
1493 // state, so it brackets the post-loop TIC clamp below. ROMS applies
1494 // that same clamp in the same place (fennel.h, just before "Update
1495 // global tracer variables"), so a divergence that appears between
1496 // G10 and G11 is in the clamp or the increment, not in whether the
1497 // clamp belongs there.
1498 tag_state("G10_POST ", parms.BioIter);
1499
1500 if (use_carbon) {
1501 for (int k = 0; k <= N; ++k) {
1502 tic(i,j,k) = amrex::min(tic(i,j,k), Real(3000.0));
1503 tic(i,j,k) = amrex::max(tic(i,j,k), Real(400.0));
1504 }
1505 }
1506
1507 // Update global tracer variables: add increment due to BGC
1508 // processes to tracer array in time index nnew. Index nnew is
1509 // the solution after advection and mixing and has transport
1510 // units (m Tunits), hence the increment is multiplied by Hz.
1511 // Subtract the original values at the top of the routine to
1512 // account only for concentrations affected by BGC processes.
1513 // If Bio were unchanged by BGC processes, the increment would
1514 // be exactly zero. Final tracer values are not bounded >=0 so
1515 // total inventory can be preserved even when advection causes
1516 // tracer concentration to go negative.
1517 const Real rmask = mskr(i,j,0);
1518 for (int k = 0; k <= N; ++k) {
1519 const Real old_no3 = amrex::max(zero, state_old(i,j,k,bio_comp.no3));
1520 const Real old_nh4 = amrex::max(zero, state_old(i,j,k,bio_comp.nh4));
1521 Real old_po4 = zero;
1522 if (use_po4) {
1523 old_po4 = amrex::max(zero, state_old(i,j,k,bio_comp.po4));
1524 }
1525 const Real old_chlo = amrex::max(zero, state_old(i,j,k,bio_comp.chlo));
1526 const Real old_phyt = amrex::max(zero, state_old(i,j,k,bio_comp.phyt));
1527 const Real old_zoop = amrex::max(zero, state_old(i,j,k,bio_comp.zoop));
1528 const Real old_lden = amrex::max(zero, state_old(i,j,k,bio_comp.lden));
1529 const Real old_sden = amrex::max(zero, state_old(i,j,k,bio_comp.sden));
1530 Real old_rden = zero;
1531 if (use_river_don) {
1532 old_rden = amrex::max(zero, state_old(i,j,k,bio_comp.rden));
1533 }
1534 Real old_rdec = zero;
1535 if (use_river_don_c) {
1536 old_rdec = amrex::max(zero, state_old(i,j,k,bio_comp.rdec));
1537 }
1538 Real old_ldec = zero;
1539 Real old_sdec = zero;
1540 Real old_tic = zero;
1541 Real old_talk = zero;
1542 if (use_carbon) {
1543 old_ldec = amrex::max(zero, state_old(i,j,k,bio_comp.ldec));
1544 old_sdec = amrex::max(zero, state_old(i,j,k,bio_comp.sdec));
1545 old_tic = amrex::max(zero, state_old(i,j,k,bio_comp.tic));
1546 old_tic = amrex::min(old_tic, Real(3000.0));
1547 old_tic = amrex::max(old_tic, Real(400.0));
1548 old_talk = amrex::max(zero, state_old(i,j,k,bio_comp.talk));
1549 }
1550 Real old_oxyg = zero;
1551 if (use_oxygen) {
1552 old_oxyg = amrex::max(zero, state_old(i,j,k,bio_comp.oxyg));
1553 }
1554 Real old_odu = zero;
1555 if (use_odu) {
1556 old_odu = amrex::max(zero, state_old(i,j,k,bio_comp.odu));
1557 }
1558
1559 state_new(i,j,k,bio_comp.no3) += (no3(i,j,k) - old_no3) * rmask * Hz(i,j,k);
1560 state_new(i,j,k,bio_comp.nh4) += (nh4(i,j,k) - old_nh4) * rmask * Hz(i,j,k);
1561 if (use_po4) {
1562 state_new(i,j,k,bio_comp.po4) += (po4(i,j,k) - old_po4) * rmask * Hz(i,j,k);
1563 }
1564 state_new(i,j,k,bio_comp.chlo) += (chlo(i,j,k) - old_chlo) * rmask * Hz(i,j,k);
1565 state_new(i,j,k,bio_comp.phyt) += (phyt(i,j,k) - old_phyt) * rmask * Hz(i,j,k);
1566 state_new(i,j,k,bio_comp.zoop) += (zoop(i,j,k) - old_zoop) * rmask * Hz(i,j,k);
1567 state_new(i,j,k,bio_comp.lden) += (lden(i,j,k) - old_lden) * rmask * Hz(i,j,k);
1568 state_new(i,j,k,bio_comp.sden) += (sden(i,j,k) - old_sden) * rmask * Hz(i,j,k);
1569 if (use_river_don) {
1570 state_new(i,j,k,bio_comp.rden) += (rden(i,j,k) - old_rden) * rmask * Hz(i,j,k);
1571 }
1572 if (use_carbon) {
1573 state_new(i,j,k,bio_comp.ldec) += (ldec(i,j,k) - old_ldec) * rmask * Hz(i,j,k);
1574 state_new(i,j,k,bio_comp.sdec) += (sdec(i,j,k) - old_sdec) * rmask * Hz(i,j,k);
1575 state_new(i,j,k,bio_comp.tic) += (tic(i,j,k) - old_tic) * rmask * Hz(i,j,k);
1576 state_new(i,j,k,bio_comp.talk) += (talk(i,j,k) - old_talk) * rmask * Hz(i,j,k);
1577 if (use_river_don) {
1578 state_new(i,j,k,bio_comp.rdec) += (rdec(i,j,k) - old_rdec) * rmask * Hz(i,j,k);
1579 }
1580 }
1581 if (use_oxygen) {
1582 state_new(i,j,k,bio_comp.oxyg) += (oxyg(i,j,k) - old_oxyg) * rmask * Hz(i,j,k);
1583 }
1584 if (use_odu) {
1585 state_new(i,j,k,bio_comp.odu) += (odu(i,j,k) - old_odu) * rmask * Hz(i,j,k);
1586 }
1587 }
1588
1589#ifdef REMORA_USE_BIOLOGY_DIAG
1590 // G11_UPDATE reads the updated global tracer array rather than
1591 // the Bio scratch, matching fennel_tag_t on the Fortran path.
1592 // Guarded separately from tag_state because it reads state_new
1593 // rather than the Bio scratch and so cannot go through it.
1594 if (dbg_level > 0 &&
1595 (dbg_level == 2 || (i == dbg_i && j == dbg_j))) {
1596 const int it = parms.BioIter;
1597 for (int k = 0; k <= N; ++k) {
1598 remora_fennel_tag("G11_UPDATE", it, i, j, k, "NO3 ",
1599 state_new(i,j,k,bio_comp.no3));
1600 remora_fennel_tag("G11_UPDATE", it, i, j, k, "NH4 ",
1601 state_new(i,j,k,bio_comp.nh4));
1602 remora_fennel_tag("G11_UPDATE", it, i, j, k, "CHLO",
1603 state_new(i,j,k,bio_comp.chlo));
1604 remora_fennel_tag("G11_UPDATE", it, i, j, k, "PHYT",
1605 state_new(i,j,k,bio_comp.phyt));
1606 remora_fennel_tag("G11_UPDATE", it, i, j, k, "ZOOP",
1607 state_new(i,j,k,bio_comp.zoop));
1608 remora_fennel_tag("G11_UPDATE", it, i, j, k, "LDEN",
1609 state_new(i,j,k,bio_comp.lden));
1610 remora_fennel_tag("G11_UPDATE", it, i, j, k, "SDEN",
1611 state_new(i,j,k,bio_comp.sden));
1612 if (use_river_don) {
1613 remora_fennel_tag("G11_UPDATE", it, i, j, k, "RDEN",
1614 state_new(i,j,k,bio_comp.rden));
1615 }
1616 if (use_po4) {
1617 remora_fennel_tag("G11_UPDATE", it, i, j, k, "PO4 ",
1618 state_new(i,j,k,bio_comp.po4));
1619 }
1620 if (use_carbon) {
1621 remora_fennel_tag("G11_UPDATE", it, i, j, k, "LDEC",
1622 state_new(i,j,k,bio_comp.ldec));
1623 remora_fennel_tag("G11_UPDATE", it, i, j, k, "SDEC",
1624 state_new(i,j,k,bio_comp.sdec));
1625 remora_fennel_tag("G11_UPDATE", it, i, j, k, "TIC ",
1626 state_new(i,j,k,bio_comp.tic));
1627 remora_fennel_tag("G11_UPDATE", it, i, j, k, "TALK",
1628 state_new(i,j,k,bio_comp.talk));
1629 if (use_river_don) {
1630 remora_fennel_tag("G11_UPDATE", it, i, j, k, "RDEC",
1631 state_new(i,j,k,bio_comp.rdec));
1632 }
1633 }
1634 if (use_oxygen) {
1635 remora_fennel_tag("G11_UPDATE", it, i, j, k, "OXYG",
1636 state_new(i,j,k,bio_comp.oxyg));
1637 }
1638 if (use_odu) {
1639 remora_fennel_tag("G11_UPDATE", it, i, j, k, "ODU ",
1640 state_new(i,j,k,bio_comp.odu));
1641 }
1642 }
1643 }
1644#endif
1645 });
1646 }
1647}
REMORA_Constants.H
two
constexpr amrex::Real two
Definition REMORA_Constants.H:8
one
constexpr amrex::Real one
Definition REMORA_Constants.H:7
zero
constexpr amrex::Real zero
Definition REMORA_Constants.H:6
REMORA_DataStruct.H
REMORA_DateClock.H
remora_caldate
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void remora_caldate(double time_ref, double current_time_days, int &year, int &month, int &day, double &yday) noexcept
Model time in days to calendar date. ROMS caldate.
Definition REMORA_DateClock.H:336
REMORA_FennelBiology.H
Temp_comp
#define Temp_comp
Definition REMORA_IndexDefines.H:10
Salt_comp
#define Salt_comp
Definition REMORA_IndexDefines.H:11
cff1
Real cff1
Definition REMORA_InitAnalyticBiology_BioToy.H:48
ldec
const int ldec
Definition REMORA_InitAnalyticBiology_BioToy.H:37
zoop
const int zoop
Definition REMORA_InitAnalyticBiology_BioToy.H:32
po4
const int po4
Definition REMORA_InitAnalyticBiology_BioToy.H:36
tic
const int tic
Definition REMORA_InitAnalyticBiology_BioToy.H:39
rden
const int rden
Definition REMORA_InitAnalyticBiology_BioToy.H:35
lden
const int lden
Definition REMORA_InitAnalyticBiology_BioToy.H:33
rdec
const int rdec
Definition REMORA_InitAnalyticBiology_BioToy.H:41
talk
const int talk
Definition REMORA_InitAnalyticBiology_BioToy.H:40
odu
const int odu
Definition REMORA_InitAnalyticBiology_BioToy.H:43
nh4
const int nh4
Definition REMORA_InitAnalyticBiology_BioToy.H:29
phyt
const int phyt
Definition REMORA_InitAnalyticBiology_BioToy.H:31
sdec
const int sdec
Definition REMORA_InitAnalyticBiology_BioToy.H:38
cff2
Real cff2
Definition REMORA_InitAnalyticBiology_BioToy.H:49
chlo
const int chlo
Definition REMORA_InitAnalyticBiology_BioToy.H:30
no3
const int no3
Definition REMORA_InitAnalyticBiology_BioToy.H:28
oxyg
const int oxyg
Definition REMORA_InitAnalyticBiology_BioToy.H:42
sden
const int sden
Definition REMORA_InitAnalyticBiology_BioToy.H:34
rho0
Real rho0
Definition REMORA_InitAnalyticSMFlux_ParticleAdvectionFlat.H:5
REMORABiology::FennelBiology::advance
void advance(const BiologyAdvanceContext &context) const override
Definition REMORA_FennelBiology.cpp:450
REMORABiology::FennelBiology::configure
void configure(const std::string &remora_prefix) override
Definition REMORA_FennelBiology.cpp:393
REMORABiology::Fennel::components
Components components(FennelParameters const &parameters, int base) noexcept
Map the Fennel tracers onto cons components starting at base.
Definition REMORA_FennelBiology.H:229
REMORABiology::PCO2AirType::secular
@ secular
secular trend plus harmonics, referenced to 1951
REMORABiology::PCO2AirType::constant
@ constant
remora.fennel.pCO2air, unchanging
REMORABiology::PCO2AirType::data
@ data
annual climatology of Laurent et al. (2017)
REMORABiology::BiologyAdvanceContext
Host views at the existing biology source-update stage.
Definition REMORA_Biology.H:64
REMORABiology::BiologyAdvanceContext::old_state
const amrex::MultiFab & old_state
Definition REMORA_Biology.H:65
REMORABiology::BiologyAdvanceContext::shortwave
const amrex::MultiFab * shortwave
Definition REMORA_Biology.H:76
REMORABiology::BiologyDeclaration::plot_aliases
amrex::Vector< std::string > plot_aliases
Definition REMORA_Biology.H:25
REMORABiology::BiologyDeclaration::tracer_names
amrex::Vector< std::string > tracer_names
Definition REMORA_Biology.H:24
REMORABiology::FennelParameters::init_params
void init_params(const std::string &remora_prefix)
Definition REMORA_FennelBiology.cpp:292