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