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