REMORA
Regional Modeling of Oceans Refined Adaptively
Loading...
Searching...
No Matches
REMORA_advance_3d.cpp
Go to the documentation of this file.
1#include <REMORA.H>
2
3using namespace amrex;
4
5/** Corresponds to step3d_uv.F and step3d_t.F in ROMS
6 *
7 * @param[in ] lev level of refinement (coarsest level is 0)
8 * @param[inout] mf_cons scalar variables (temp, salt, scalar, etc)
9 * @param[inout] mf_u u-velocity
10 * @param[inout] mf_v v-velocity
11 * @param[inout] mf_sstore intermediate-step scalar variables
12 * @param[inout] mf_ru RHS of total u-velocity
13 * @param[inout] mf_rv RHS of total v-velocity
14 * @param[in ] mf_DU_avg1 time-averaged u-flux for 2D equations
15 * @param[in ] mf_DU_avg2 time-averaged u-flux for 3D equation coupling
16 * @param[in ] mf_DV_avg1 time-averaged v-flux for 2D equations
17 * @param[in ] mf_DV_avg2 time-averaged v-flux for 3D equation coupling
18 * @param[inout] mf_ubar vertically integrated u-momentum
19 * @param[inout] mf_vbar vertically integrated v-momentum
20 * @param[inout] mf_Akv vertical viscosity coefficient
21 * @param[inout] mf_Akt vertical diffusivity coefficient
22 * @param[inout] mf_Hz vertical height of cells
23 * @param[inout] mf_Huon u-volume flux
24 * @param[inout] mf_Huon v-volume flux
25 * @param[inout] mf_z_w vertical coordinates on w points
26 * @param[in ] mf_h bathymetry
27 * @param[in ] mf_pm 1 / dx
28 * @param[in ] mf_pn 1 / dy
29 * @param[in ] mf_msku land-sea mask at u-points
30 * @param[in ] mf_mskv land-sea mask at v-points
31 * @param[in ] N number of vertical levels
32 * @param[in ] dt_lev time step at this refinement level
33 */
34void
35REMORA::advance_3d (int lev, MultiFab& mf_cons,
36 MultiFab& mf_u , MultiFab& mf_v ,
37 MultiFab* mf_sstore,
38 MultiFab* mf_ru , MultiFab* mf_rv,
39 std::unique_ptr<MultiFab>& mf_DU_avg1,
40 std::unique_ptr<MultiFab>& mf_DU_avg2,
41 std::unique_ptr<MultiFab>& mf_DV_avg1,
42 std::unique_ptr<MultiFab>& mf_DV_avg2,
43 std::unique_ptr<MultiFab>& mf_ubar,
44 std::unique_ptr<MultiFab>& mf_vbar,
45 std::unique_ptr<MultiFab>& mf_Akv,
46 std::unique_ptr<MultiFab>& mf_Akt,
47 std::unique_ptr<MultiFab>& mf_Hz,
48 std::unique_ptr<MultiFab>& mf_Huon,
49 std::unique_ptr<MultiFab>& mf_Hvom,
50 std::unique_ptr<MultiFab>& mf_z_w,
51 MultiFab const* mf_h,
52 MultiFab const* mf_pm,
53 MultiFab const* mf_pn,
54 MultiFab const* mf_mskr,
55 MultiFab const* mf_msku,
56 MultiFab const* mf_mskv,
57 const int N, Real dt_lev)
58{
59 BL_PROFILE("REMORA::advance_3d()");
60 const int nrhs = 0;
61 int nnew = 0;
62
63 int iic = istep[lev];
64 int ntfirst = 0;
65
66 // Because zeta may have changed
67 stretch_transform(lev);
68
69 // These temporaries used to be made in advance_3d_ml and passed in;
70 // now we make them here
71
72 const BoxArray& ba = mf_cons.boxArray();
73 const DistributionMapping& dm = mf_cons.DistributionMap();
74
75 //Only used locally, probably should be rearranged into FArrayBox declaration
76 MultiFab mf_AK (convert(ba,IntVect(0,0,1)),dm,1,IntVect(NGROW,NGROW,0)); //2d missing j coordinate
77 MultiFab mf_DC (ba,dm,1,IntVect(NGROW,NGROW,NGROW-1)); //2d missing j coordinate
78 MultiFab mf_Hzk(ba,dm,1,IntVect(NGROW,NGROW,NGROW-1)); //2d missing j coordinate
79
80 for ( MFIter mfi(mf_cons, TilingIfNotGPU()); mfi.isValid(); ++mfi )
81 {
82 Array4<Real > const& u = mf_u.array(mfi);
83 Array4<Real > const& v = mf_v.array(mfi);
84
85 Array4<Real > const& ru = mf_ru->array(mfi);
86 Array4<Real > const& rv = mf_rv->array(mfi);
87
88 Array4<Real > const& AK = mf_AK.array(mfi);
89 Array4<Real > const& DC = mf_DC.array(mfi);
90
91 Array4<Real > const& Hzk = mf_Hzk.array(mfi);
92 Array4<Real > const& Akv = mf_Akv->array(mfi);
93
94 Array4<Real const> const& Hz = mf_Hz->const_array(mfi);
95
96 Array4<Real const> const& DU_avg1 = mf_DU_avg1->const_array(mfi);
97 Array4<Real const> const& DV_avg1 = mf_DV_avg1->const_array(mfi);
98
99 Array4<Real const> const& pm = mf_pm->const_array(mfi);
100 Array4<Real const> const& pn = mf_pn->const_array(mfi);
101
102 Array4<Real const> const& msku = mf_msku->const_array(mfi);
103 Array4<Real const> const& mskv = mf_mskv->const_array(mfi);
104
105 Box bx = mfi.tilebox();
106 Box gbx2 = mfi.growntilebox(IntVect(NGROW,NGROW,0));
107 Box gbx21 = mfi.growntilebox(IntVect(NGROW,NGROW,NGROW-1));
108
109 Box xbx = mfi.nodaltilebox(0);
110 Box ybx = mfi.nodaltilebox(1);
111
112 Box gbx2D = gbx2;
113 gbx2D.makeSlab(2,0);
114
115 Box tbxp1 = bx;
116 Box tbxp11 = bx;
117 Box tbxp2 = bx;
118 tbxp1.grow(IntVect(NGROW-1,NGROW-1,0));
119 tbxp2.grow(IntVect(NGROW,NGROW,0));
120 tbxp11.grow(IntVect(NGROW-1,NGROW-1,NGROW-1));
121
122 FArrayBox fab_FC(convert(gbx2,IntVect(0,0,1)),1,amrex::The_Async_Arena());
123 FArrayBox fab_BC(convert(gbx2,IntVect(0,0,1)),1,amrex::The_Async_Arena());
124 FArrayBox fab_CF(gbx21,1,amrex::The_Async_Arena());
125 FArrayBox fab_W(tbxp2,1,amrex::The_Async_Arena());
126
127 auto FC = fab_FC.array();
128 auto BC = fab_BC.array();
129 auto CF = fab_CF.array();
130
131 Real cff;
132 if (iic==ntfirst) {
133 cff=Real(0.25)*dt_lev;
134 } else if (iic==ntfirst+1) {
135 cff=Real(0.25)*dt_lev*Real(3.0)/Real(2.0);
136 } else {
137 cff=Real(0.25)*dt_lev*Real(23.0)/Real(12.0);
138 }
139
140 ParallelFor(xbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
141 {
142 u(i,j,k) += cff * (pm(i,j,0)+pm(i-1,j,0)) * (pn(i,j,0)+pn(i-1,j,0)) * ru(i,j,k,nrhs);
143 u(i,j,k) *= two / (Hz(i-1,j,k) + Hz(i,j,k));
144 });
145
146 ParallelFor(ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
147 {
148 v(i,j,k) += cff * (pm(i,j,0)+pm(i,j-1,0)) * (pn(i,j,0)+pn(i,j-1,0)) * rv(i,j,k,nrhs);
149 v(i,j,k) *= two / (Hz(i,j-1,k) + Hz(i,j,k));
150 });
151
152 // NOTE: DC is only used as scratch in vert_visc_3d -- no need to pass or return a value
153 // NOTE: may not actually need to set these to zero
154
155 // Reset to zero on the box on which they'll be used
156 mf_DC[mfi].template setVal<RunOn::Device>(zero,xbx);
157 fab_CF.template setVal<RunOn::Device>(zero,xbx);
158
159 vert_visc_3d(xbx,1,0,u,Hz,Hzk,AK,Akv,BC,DC,FC,CF,nnew,N,dt_lev);
160
161 // Reset to zero on the box on which they'll be used
162 mf_DC[mfi].template setVal<RunOn::Device>(zero,ybx);
163 fab_CF.template setVal<RunOn::Device>(zero,ybx);
164
165 vert_visc_3d(ybx,0,1,v,Hz,Hzk,AK,Akv,BC,DC,FC,CF,nnew,N,dt_lev);
166
167 // Reset to zero on the box on which they'll be used
168 mf_DC[mfi].template setVal<RunOn::Device>(zero,xbx);
169 fab_CF.template setVal<RunOn::Device>(zero,xbx);
170
171 vert_mean_3d(xbx,1,0,u,Hz,DU_avg1,DC,CF,pn,msku,nnew,N);
172
173 // Reset to zero on the box on which they'll be used
174 mf_DC[mfi].template setVal<RunOn::Device>(zero,ybx);
175 fab_CF.template setVal<RunOn::Device>(zero,ybx);
176
177 vert_mean_3d(ybx,0,1,v,Hz,DV_avg1,DC,CF,pm,mskv,nnew,N);
178 }
179
180 // Apply physical boundary conditions to u and v
181 (*physbcs[lev])(mf_u,*mf_msku,0,1,mf_u.nGrowVect(),t_old[lev],xvel_bc(),0,*xvel_old[lev]);
182 (*physbcs[lev])(mf_v,*mf_mskv,0,1,mf_v.nGrowVect(),t_old[lev],yvel_bc(),0,*yvel_old[lev]);
183
184#ifdef REMORA_USE_NETCDF
185 // Fill the data which is stored in the boundary data read from netcdf files
186 if (solverChoice.boundary_from_netcdf)
187 {
188 fill_from_bdyfiles(lev, mf_u,*mf_msku,t_old[lev],xvel_bc(),BdyVars::u,0,0,*xvel_old[lev],dt_lev);
189 fill_from_bdyfiles(lev, mf_v,*mf_mskv,t_old[lev],yvel_bc(),BdyVars::v,0,0,*yvel_old[lev],dt_lev);
190 }
191
192 if (solverChoice.do_rivers) {
193 river_source_transport->update_interpolated_to_time(t_old[lev]);
194 }
195#endif
196 for ( MFIter mfi(mf_cons, TilingIfNotGPU()); mfi.isValid(); ++mfi )
197 {
198 Array4<Real > const& u = mf_u.array(mfi);
199 Array4<Real > const& v = mf_v.array(mfi);
200
201 Array4<Real > const& DC = mf_DC.array(mfi);
202
203 Array4<Real const> const& Hz = mf_Hz->const_array(mfi);
204
205 Array4<Real const> const& DU_avg1 = mf_DU_avg1->const_array(mfi);
206 Array4<Real const> const& DV_avg1 = mf_DV_avg1->const_array(mfi);
207
208 Array4<Real const> const& DU_avg2 = mf_DU_avg2->const_array(mfi);
209 Array4<Real const> const& DV_avg2 = mf_DV_avg2->const_array(mfi);
210
211 Array4<Real> const& ubar = mf_ubar->array(mfi);
212 Array4<Real> const& vbar = mf_vbar->array(mfi);
213
214 Array4<Real> const& Huon = mf_Huon->array(mfi);
215 Array4<Real> const& Hvom = mf_Hvom->array(mfi);
216
217 Array4<Real const> const& pm = mf_pm->const_array(mfi);
218 Array4<Real const> const& pn = mf_pn->const_array(mfi);
219 Array4<Real const> const& msku = mf_msku->const_array(mfi);
220 Array4<Real const> const& mskv = mf_mskv->const_array(mfi);
221
222 Box gbx2 = mfi.growntilebox(IntVect(NGROW,NGROW,0));
223
224 FArrayBox fab_FC(gbx2,1,amrex::The_Async_Arena());
225 auto FC = fab_FC.array();
226
227#ifdef REMORA_USE_NETCDF
228 if (solverChoice.do_rivers) {
229 Box gbx1 = mfi.growntilebox(IntVect(NGROW-1,NGROW-1,0));
230 Array4<Real const> const& z_w = mf_z_w->const_array(mfi);
231 Array4<int const> const& river_pos = vec_river_position[lev]->const_array(mfi);
232 Array4<Real const> const& river_transport = river_source_transport->fab_interp->array();
233 int* river_direction_d = river_direction.data();
234 ParallelFor(gbx1, [=] AMREX_GPU_DEVICE(int i, int j, int k)
235 {
236 int iriver = river_pos(i,j,0);
237 if (iriver >= 0) {
238 if (river_direction_d[iriver] == 0) {
239 Real on_u = two / (pn(i,j,0)+pn(i-1,j,0));
240 Real cff = one / (on_u * Real(0.5) * (z_w(i-1,j,k+1) - z_w(i-1,j,k) + z_w(i,j,k+1) - z_w(i,j,k)));
241 u(i,j,k) = cff * river_transport(iriver,0,k);
242 } else {
243 Real om_v = two / (pm(i,j,0)+pm(i,j-1,0));
244 Real cff = one / (om_v * Real(0.5) * (z_w(i,j-1,k+1) - z_w(i,j-1,k) + z_w(i,j,k+1) - z_w(i,j,k)));
245 v(i,j,k) = cff * river_transport(iriver,0,k);
246 }
247 }
248 });
249 }
250#endif
251
252#if 0
253 // Reset to zero on the box on which they'll be used
254 mf_DC[mfi].template setVal<RunOn::Device>(zero,grow(xbx,IntVect(0,0,1)));
255 fab_CF.template setVal<RunOn::Device>(zero,grow(xbx,IntVect(0,0,1)));
256
257 update_massflux_3d(xbx,1,0,u,ubar,Huon,Hz,pn,DU_avg1,DU_avg2,DC,FC,nnew);
258
259 // Reset to zero on the box on which they'll be used
260 mf_DC[mfi].template setVal<RunOn::Device>(zero,grow(ybx,IntVect(0,0,1)));
261 fab_CF.template setVal<RunOn::Device>(zero,grow(ybx,IntVect(0,0,1)));
262
263 update_massflux_3d(ybx,0,1,v,vbar,Hvom,Hz,pm,DV_avg1,DV_avg2,DC,FC,nnew);
264
265#else
266
267 // Reset to zero on the box on which they'll be used
268 fab_FC.template setVal<RunOn::Device>(zero,gbx2);
269 mf_DC[mfi].template setVal<RunOn::Device>(zero,grow(gbx2,IntVect(0,0,1)));
270 update_massflux_3d(lev,gbx2,1,0,u,ubar,Huon,Hz,pn,DU_avg1,DU_avg2,DC,FC,msku,nnew);
271
272 // Reset to zero on the box on which they'll be used
273 fab_FC.template setVal<RunOn::Device>(zero,gbx2);
274 mf_DC[mfi].template setVal<RunOn::Device>(zero,grow(gbx2,IntVect(0,0,1)));
275 update_massflux_3d(lev,gbx2,0,1,v,vbar,Hvom,Hz,pm,DV_avg1,DV_avg2,DC,FC,mskv,nnew);
276
277#endif
278 }
279
280 // WE BELIEVE THESE VALUES SHOULD ALREADY BE FILLED
281 // mf_Huon->FillBoundary(geom[lev].periodicity());
282 // mf_Hvom->FillBoundary(geom[lev].periodicity());
283
284 // ************************************************************************
285 // This should fill both temp and salt with temp/salt currently in cons_old
286 // ************************************************************************
287
288 MultiFab mf_W(convert(ba,IntVect(0,0,1)),dm,1,IntVect(NGROW+1,NGROW+1,0));
289 mf_W.setVal(zero);
290 for ( MFIter mfi(mf_cons, TilingIfNotGPU()); mfi.isValid(); ++mfi )
291 {
292 Array4<Real> const& Huon = mf_Huon->array(mfi);
293 Array4<Real> const& Hvom = mf_Hvom->array(mfi);
294
295 Array4<Real const> const& z_w = mf_z_w->const_array(mfi);
296 Array4<Real const> const& h = mf_h->const_array(mfi);
297
298 Box bx = mfi.tilebox();
299 Box gbx1 = mfi.growntilebox(IntVect(NGROW-1,NGROW-1,0));
300 Box gbx2 = mfi.growntilebox(IntVect(NGROW,NGROW,0));
301 Box gbx21 = mfi.growntilebox(IntVect(NGROW,NGROW,NGROW-1));
302
303 Box tbxp1 = bx;
304 Box tbxp11 = bx;
305 Box tbxp2 = bx;
306 tbxp1.grow(IntVect(NGROW-1,NGROW-1,0));
307 tbxp2.grow(IntVect(NGROW,NGROW,0));
308 tbxp11.grow(IntVect(NGROW-1,NGROW-1,NGROW-1));
309
310 FArrayBox fab_FC(surroundingNodes(gbx2,2),1,amrex::The_Async_Arena());
311 FArrayBox fab_BC(gbx2,1,amrex::The_Async_Arena());
312 FArrayBox fab_CF(gbx21,1,amrex::The_Async_Arena());
313
314 auto W = mf_W.array(mfi);
315
316 //
317 //------------------------------------------------------------------------
318 // Vertically integrate horizontal mass flux divergence.
319 //------------------------------------------------------------------------
320 //
321 //Should really use gbx3uneven
322 //TODO: go over these boxes and compare to other spots where we do the same thing
323 Box gbx1D = gbx1;
324 gbx1D.makeSlab(2,0);
325
326 ParallelFor(gbx1D, N+1,
327 [=] AMREX_GPU_DEVICE (int i, int j, int , int kk)
328 {
329 // Starting with zero vertical velocity at the bottom, integrate
330 // from the bottom (k=0) to the free-surface (k=N). The w(:,:,N(ng))
331 // contains the vertical velocity at the free-surface, d(zeta)/d(t).
332 // Notice that barotropic mass flux divergence is not used directly.
333 //
334 int k = kk + 1;
335 W(i,j,k) = W(i,j,k-1) - (Huon(i+1,j,k-1)-Huon(i,j,k-1)) - (Hvom(i,j+1,k-1)-Hvom(i,j,k-1));
336 });
337 ParallelFor(gbx1D, [=] AMREX_GPU_DEVICE (int i, int j, int )
338 {
339 W(i,j,N+1)=W(i,j,N+1)/(z_w(i,j,N+1)+h(i,j,0,0)); // wrk_i
340 });
341 ParallelFor(gbx1D, N,
342 [=] AMREX_GPU_DEVICE (int i, int j, int , int kk)
343 {
344 int k = kk + 1;
345 W(i,j,k) = W(i,j,k)- W(i,j,N+1)*(z_w(i,j,k)+h(i,j,0,0));
346 });
347 ParallelFor(gbx1D, [=] AMREX_GPU_DEVICE (int i, int j, int )
348 {
349 W(i,j,N+1) = zero;
350 });
351 }
352
353 const int nstp = (iic) % 2;
354 nnew = 1-nstp;
355 if (solverChoice.vert_mixing_type == VertMixingType::GLS) {
356 gls_corrector(lev, vec_gls[lev].get(), vec_tke[lev].get(), mf_W, vec_Akv[lev].get(),
357 vec_Akt[lev].get(),vec_Akk[lev].get(), vec_Akp[lev].get(), vec_mskr[lev].get(),
358 vec_msku[lev].get(), vec_mskv[lev].get(),
359 nstp, nnew, N, dt_lev);
360 }
361 advance_biology(lev, *cons_old[lev], mf_cons, N, dt_lev);
362 nnew = 0;
363
364 for ( MFIter mfi(mf_cons, TilingIfNotGPU()); mfi.isValid(); ++mfi )
365 {
366 Array4<Real> const& Hz = mf_Hz->array(mfi);
367
368 Array4<Real> const& Huon = mf_Huon->array(mfi);
369 Array4<Real> const& Hvom = mf_Hvom->array(mfi);
370
371 Array4<Real const> const& pm = mf_pm->const_array(mfi);
372 Array4<Real const> const& pn = mf_pn->const_array(mfi);
373 Array4<Real const> const& mskr = mf_mskr->const_array(mfi);
374 Array4<Real const> const& msku = mf_msku->const_array(mfi);
375 Array4<Real const> const& mskv = mf_mskv->const_array(mfi);
376
377 Box bx = mfi.tilebox();
378 Box gbx2 = mfi.growntilebox(IntVect(NGROW,NGROW,0));
379 Box gbx21 = mfi.growntilebox(IntVect(NGROW,NGROW,NGROW-1));
380
381 Box tbxp1 = bx;
382 Box tbxp11 = bx;
383 Box tbxp2 = bx;
384 tbxp1.grow(IntVect(NGROW-1,NGROW-1,0));
385 tbxp2.grow(IntVect(NGROW,NGROW,0));
386 tbxp11.grow(IntVect(NGROW-1,NGROW-1,NGROW-1));
387
388 FArrayBox fab_FC(surroundingNodes(gbx2,2),1,amrex::The_Async_Arena());
389 FArrayBox fab_BC(gbx2,1,amrex::The_Async_Arena());
390 FArrayBox fab_CF(gbx21,1,amrex::The_Async_Arena());
391
392 auto FC = fab_FC.array();
393 auto W = mf_W.array(mfi);
394 //
395 //-----------------------------------------------------------------------
396 // rhs_t_3d
397 //-----------------------------------------------------------------------
398 //
399 for (int i_comp=0; i_comp < ncons; i_comp++)
400 {
401#ifdef REMORA_USE_NETCDF
402 FArrayBox* fab_river_source;
403 if (solverChoice.do_rivers_cons[i_comp]) {
404 fab_river_source = river_source_cons[i_comp]->fab_interp;
405 }
406 const Array4<const int>& river_pos = (solverChoice.do_rivers) ? vec_river_position[lev]->const_array(mfi) : Array4<const int>();
407 const Array4<const Real>& river_source = (solverChoice.do_rivers_cons[i_comp]) ? fab_river_source->const_array() : Array4<const Real>();
408#else
409 const Array4<const int >& river_pos = Array4<const int>();
410 const Array4<const Real>& river_source = Array4<const Real>();
411#endif
412 Array4<Real> const& sstore = mf_sstore->array(mfi, i_comp);
413 rhs_t_3d(lev,bx, mf_cons.array(mfi,i_comp), sstore, Huon, Hvom,
414 Hz, pn, pm, W, FC, mskr, msku, mskv, river_pos, river_source, nrhs, nnew, N,dt_lev);
415 }
416 } // mfi
417
418 FillPatch(lev, t_old[lev], mf_cons, cons_new, BCVars::cons_bc, BdyVars::t,0,true,false,0,0,dt_lev,*cons_old[lev]);
419
420 for ( MFIter mfi(mf_cons, TilingIfNotGPU()); mfi.isValid(); ++mfi )
421 {
422 Array4<Real> const& AK = mf_AK.array(mfi);
423 Array4<Real> const& DC = mf_DC.array(mfi);
424
425 Array4<Real> const& Hzk = mf_Hzk.array(mfi);
426 Array4<Real const> const& Hz = mf_Hz->const_array(mfi);
427
428 Box bx = mfi.tilebox();
429
430 // Copy the tilebox
431 Box tbxp1 = bx;
432 Box tbxp11 = bx;
433 Box tbxp2 = bx;
434 Box tbxp21 = bx;
435 //make only gbx be grown to match multifabs
436 tbxp21.grow(IntVect(NGROW,NGROW,NGROW-1));
437 tbxp2.grow(IntVect(NGROW,NGROW,0));
438 tbxp1.grow(IntVect(NGROW-1,NGROW-1,0));
439 tbxp11.grow(IntVect(NGROW-1,NGROW-1,NGROW-1));
440
441 FArrayBox fab_FC(tbxp2,1,amrex::The_Async_Arena());
442 FArrayBox fab_BC(tbxp2,1,amrex::The_Async_Arena());
443 FArrayBox fab_CF(tbxp21,1,amrex::The_Async_Arena());
444 FArrayBox fab_W(tbxp2,1,amrex::The_Async_Arena());
445
446 auto FC = fab_FC.array();
447 auto BC = fab_BC.array();
448 auto CF = fab_CF.array();
449
450 for (int i_comp=0; i_comp < ncons; i_comp++) {
451 vert_visc_3d(bx,0,0,mf_cons.array(mfi,i_comp),Hz,Hzk,
452 AK,mf_Akt->array(mfi,akt_comp(i_comp)),BC,DC,FC,CF,nnew,N,dt_lev);
453 }
454 } // MFiter
455 FillPatch(lev, t_old[lev], *cons_new[lev], cons_new, BCVars::cons_bc, BdyVars::t,0,true,false,0,0,dt_lev,*cons_old[lev]);
456
457#ifdef REMORA_USE_NETCDF
458 // Nudge every tracer that was given a climatology toward it. temp and salt are just
459 // components 0 and 1 of this loop, so their behavior is unchanged.
460 for (int icomp = 0; icomp < ncons; ++icomp) {
461 if (!solverChoice.do_cons_clim_nudg[icomp]) { continue; }
462
463 cons_clim_data_from_file[icomp]->update_interpolated_to_time(t_old[lev], lev, cons_new[lev], geom, ref_ratio);
464
465 for ( MFIter mfi(mf_cons, TilingIfNotGPU()); mfi.isValid(); ++mfi )
466 {
467 Box bx = mfi.growntilebox(IntVect(1,1,0));
468 Array4<const Real> const& cons_nudg_coeff = vec_nudg_coeff[BdyVars::cons(icomp)][lev]->const_array(mfi);
469 Array4<const Real> const& cons_clim = cons_clim_data_from_file[icomp]->get_interpolated_mf(lev)->const_array(mfi);
470 Array4< Real> const& cons = cons_new[lev]->array(mfi, icomp);
471 Array4<const Real> const& Hz = vec_Hz[lev]->const_array(mfi);
472 Array4<const Real> const& pm = vec_pm[lev]->const_array(mfi);
473 Array4<const Real> const& pn = vec_pn[lev]->const_array(mfi);
474
475 apply_clim_nudg(bx, 0, 0, cons, cons, cons_clim, cons_nudg_coeff, Hz, pm, pn, dt_lev);
476 }
477 }
478#endif
479}
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
VertMixingType::GLS
@ GLS
NGROW
#define NGROW
Definition REMORA_IndexDefines.H:30
akt_comp
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE int akt_comp(int icomp) noexcept
Definition REMORA_IndexDefines.H:28
setVal
mf_h setVal(geomdata.ProbHi(2))
REMORA::ncons
int ncons
Number of conserved scalars in the state (temperature + salt + passive scalars + biology tracers)
Definition REMORA.H:1644
REMORA::update_massflux_3d
void update_massflux_3d(int lev, const amrex::Box &bx, const int ioff, const int joff, const amrex::Array4< amrex::Real > &phi, const amrex::Array4< amrex::Real > &phibar, const amrex::Array4< amrex::Real > &Hphi, const amrex::Array4< amrex::Real const > &Hz, const amrex::Array4< amrex::Real const > &pm_or_pn, const amrex::Array4< amrex::Real const > &Dphi1, const amrex::Array4< amrex::Real const > &Dphi2, const amrex::Array4< amrex::Real > &DC, const amrex::Array4< amrex::Real > &FC, const amrex::Array4< amrex::Real const > &msk, const int nnew)
Correct mass flux.
Definition REMORA_update_massflux_3d.cpp:23
REMORA::xvel_bc
int xvel_bc() const noexcept
Definition REMORA.H:1313
REMORA::vec_pm
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_pm
horizontal scaling factor: 1 / dx (2D)
Definition REMORA.H:568
REMORA::cons_new
amrex::Vector< amrex::MultiFab * > cons_new
multilevel data container for current step's scalar data: temperature, salinity, passive tracer
Definition REMORA.H:386
REMORA::stretch_transform
void stretch_transform(int lev)
Calculate vertical stretched coordinates.
Definition REMORA_DepthStretchTransform.H:78
REMORA::vec_mskr
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_mskr
land/sea mask at cell centers (2D)
Definition REMORA.H:557
REMORA::advance_3d
void advance_3d(int lev, amrex::MultiFab &mf_cons, amrex::MultiFab &mf_u, amrex::MultiFab &mf_v, amrex::MultiFab *mf_sstore, amrex::MultiFab *mf_ru, amrex::MultiFab *mf_rv, std::unique_ptr< amrex::MultiFab > &mf_DU_avg1, std::unique_ptr< amrex::MultiFab > &mf_DU_avg2, std::unique_ptr< amrex::MultiFab > &mf_DV_avg1, std::unique_ptr< amrex::MultiFab > &mf_DV_avg2, std::unique_ptr< amrex::MultiFab > &mf_ubar, std::unique_ptr< amrex::MultiFab > &mf_vbar, std::unique_ptr< amrex::MultiFab > &mf_Akv, std::unique_ptr< amrex::MultiFab > &mf_Akt, std::unique_ptr< amrex::MultiFab > &mf_Hz, std::unique_ptr< amrex::MultiFab > &mf_Huon, std::unique_ptr< amrex::MultiFab > &mf_Hvom, std::unique_ptr< amrex::MultiFab > &mf_z_w, amrex::MultiFab const *mf_h, amrex::MultiFab const *mf_pm, amrex::MultiFab const *mf_pn, amrex::MultiFab const *mf_mskr, amrex::MultiFab const *mf_msku, amrex::MultiFab const *mf_mskv, const int N, const amrex::Real dt_lev)
Advance the 3D variables.
Definition REMORA_advance_3d.cpp:35
REMORA::vec_tke
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_tke
Turbulent kinetic energy.
Definition REMORA.H:632
REMORA::vec_gls
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_gls
Turbulent generic length scale.
Definition REMORA.H:634
REMORA::yvel_bc
int yvel_bc() const noexcept
Definition REMORA.H:1314
REMORA::fill_from_bdyfiles
void fill_from_bdyfiles(int lev, amrex::MultiFab &mf_to_fill, const amrex::MultiFab &mf_mask, const amrex::Real time, const int bccomp, const int bdy_var_type, const int icomp_to_fill, const int icomp_calc=0, const amrex::MultiFab &mf_calc=amrex::MultiFab(), const amrex::Real=zero)
Fill boundary data from netcdf file.
Definition REMORA_BoundaryConditions_netcdf.cpp:20
REMORA::vert_mean_3d
void vert_mean_3d(const amrex::Box &bx, const int ioff, const int joff, const amrex::Array4< amrex::Real > &phi, const amrex::Array4< amrex::Real const > &Hz, const amrex::Array4< amrex::Real const > &Dphi_avg1, const amrex::Array4< amrex::Real > &DC, const amrex::Array4< amrex::Real > &CF, const amrex::Array4< amrex::Real const > &pm_or_pn, const amrex::Array4< amrex::Real const > &msk, const int nnew, const int N)
Adjust 3D momentum variables based on vertical mean momentum.
Definition REMORA_vert_mean_3d.cpp:19
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_Akt
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Akt
Vertical diffusion coefficient (3D)
Definition REMORA.H:432
REMORA::vec_msku
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_msku
land/sea mask at x-faces (2D)
Definition REMORA.H:559
REMORA::xvel_old
amrex::Vector< amrex::MultiFab * > xvel_old
multilevel data container for last step's x velocities (u in ROMS)
Definition REMORA.H:379
REMORA::vec_mskv
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_mskv
land/sea mask at y-faces (2D)
Definition REMORA.H:561
REMORA::physbcs
amrex::Vector< std::unique_ptr< REMORAPhysBCFunct > > physbcs
Vector (over level) of functors to apply physical boundary conditions.
Definition REMORA.H:1566
REMORA::istep
amrex::Vector< int > istep
which step?
Definition REMORA.H:1552
REMORA::yvel_old
amrex::Vector< amrex::MultiFab * > yvel_old
multilevel data container for last step's y velocities (v in ROMS)
Definition REMORA.H:381
REMORA::rhs_t_3d
void rhs_t_3d(int lev, const amrex::Box &bx, const amrex::Array4< amrex::Real > &t, const amrex::Array4< amrex::Real const > &tempstore, const amrex::Array4< amrex::Real const > &Huon, const amrex::Array4< amrex::Real const > &Hvom, const amrex::Array4< amrex::Real const > &Hz, const amrex::Array4< amrex::Real const > &pn, const amrex::Array4< amrex::Real const > &pm, const amrex::Array4< amrex::Real const > &W, const amrex::Array4< amrex::Real > &FC, const amrex::Array4< amrex::Real const > &mskr, const amrex::Array4< amrex::Real const > &msku, const amrex::Array4< amrex::Real const > &mskv, const amrex::Array4< int const > &river_pos, const amrex::Array4< amrex::Real const > &river_source, int nrhs, int nnew, int N, const amrex::Real dt_lev)
RHS terms for tracer.
Definition REMORA_rhs_t_3d.cpp:28
REMORA::apply_clim_nudg
void apply_clim_nudg(const amrex::Box &bx, int ioff, int joff, const amrex::Array4< amrex::Real > &var, const amrex::Array4< amrex::Real const > &var_old, const amrex::Array4< amrex::Real const > &var_clim, const amrex::Array4< amrex::Real const > &clim_coeff, const amrex::Array4< amrex::Real const > &Hz, const amrex::Array4< amrex::Real const > &pm, const amrex::Array4< amrex::Real const > &pn, const amrex::Real dt_lev=zero)
Apply climatology nudging.
Definition REMORA_apply_clim_nudg.cpp:19
REMORA::solverChoice
static SolverChoice solverChoice
Container for algorithmic choices.
Definition REMORA.H:1717
REMORA::vec_river_position
amrex::Vector< std::unique_ptr< amrex::iMultiFab > > vec_river_position
iMultiFab for river positions; contents are indices of rivers
Definition REMORA.H:1494
REMORA::vec_Akk
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Akk
Turbulent kinetic energy vertical diffusion coefficient.
Definition REMORA.H:638
REMORA::vert_visc_3d
void vert_visc_3d(const amrex::Box &bx, const int ioff, const int joff, const amrex::Array4< amrex::Real > &phi, const amrex::Array4< amrex::Real const > &Hz, const amrex::Array4< amrex::Real > &Hzk, const amrex::Array4< amrex::Real > &AK, const amrex::Array4< amrex::Real const > &Akv, const amrex::Array4< amrex::Real > &BC, const amrex::Array4< amrex::Real > &DC, const amrex::Array4< amrex::Real > &FC, const amrex::Array4< amrex::Real > &CF, const int nnew, const int N, const amrex::Real dt_lev)
Calculate effects of vertical viscosity or diffusivity.
Definition REMORA_vert_visc_3d.cpp:23
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::river_direction
amrex::Gpu::DeviceVector< int > river_direction
Vector over rivers of river direction: 0: u-face; 1: v-face; 2: w-face.
Definition REMORA.H:1496
REMORA::cons_old
amrex::Vector< amrex::MultiFab * > cons_old
multilevel data container for last step's scalar data: temperature, salinity, passive tracer
Definition REMORA.H:377
REMORA::FillPatch
void FillPatch(int lev, amrex::Real time, amrex::MultiFab &mf_to_be_filled, amrex::Vector< amrex::MultiFab * > const &mfs, const int bccomp, const int bdy_var_type=BdyVars::null, const int icomp=0, const bool fill_all=true, const bool fill_set=false, const int n_not_fill=0, const int icomp_calc=0, const amrex::Real dt=zero, const amrex::MultiFab &mf_calc=amrex::MultiFab())
Fill a new MultiFab by copying in phi from valid region and filling ghost cells.
Definition REMORA_FillPatch.cpp:30
REMORA::cons_clim_data_from_file
amrex::Vector< std::unique_ptr< NCTimeSeries > > cons_clim_data_from_file
Vector over cons components of climatology data read from file.
Definition REMORA.H:1480
REMORA::river_source_cons
amrex::Vector< std::unique_ptr< NCTimeSeriesRiver > > river_source_cons
Vector of data containers for scalar data in rivers.
Definition REMORA.H:1483
REMORA::river_source_transport
std::unique_ptr< NCTimeSeriesRiver > river_source_transport
Data container for momentum transport in rivers.
Definition REMORA.H:1485
REMORA::vec_nudg_coeff
amrex::Vector< amrex::Vector< std::unique_ptr< amrex::MultiFab > > > vec_nudg_coeff
Climatology nudging coefficients.
Definition REMORA.H:643
REMORA::vec_pn
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_pn
horizontal scaling factor: 1 / dy (2D)
Definition REMORA.H:570
REMORA::vec_Akv
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Akv
Vertical viscosity coefficient (3D)
Definition REMORA.H:430
REMORA::t_old
amrex::Vector< amrex::Real > t_old
old time at each level
Definition REMORA.H:1558
REMORA::vec_Akp
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Akp
Turbulent length scale vertical diffusion coefficient.
Definition REMORA.H:640
REMORA::gls_corrector
void gls_corrector(int lev, amrex::MultiFab *mf_gls, amrex::MultiFab *mf_tke, amrex::MultiFab &mf_W, amrex::MultiFab *mf_Akv, amrex::MultiFab *mf_Akt, amrex::MultiFab *mf_Akk, amrex::MultiFab *mf_Akp, amrex::MultiFab *mf_mskr, amrex::MultiFab *mf_msku, amrex::MultiFab *mf_mskv, const int nstp, const int nnew, const int N, const amrex::Real dt_lev)
Corrector step for GLS calculation.
Definition REMORA_gls.cpp:247
BCVars::cons_bc
static constexpr int cons_bc
Definition REMORA_IndexDefines.H:36
BdyVars::t
static constexpr int t
cons component Temp_comp
Definition REMORA_IndexDefines.H:73
BdyVars::u
static constexpr int u
Definition REMORA_IndexDefines.H:71
BdyVars::v
static constexpr int v
Definition REMORA_IndexDefines.H:72
BdyVars::cons
int cons(int icomp) noexcept
Definition REMORA_IndexDefines.H:76
SolverChoice::do_rivers_cons
amrex::Vector< int > do_rivers_cons
Definition REMORA_DataStruct.H:900
SolverChoice::vert_mixing_type
VertMixingType vert_mixing_type
Definition REMORA_DataStruct.H:936
SolverChoice::do_cons_clim_nudg
amrex::Vector< int > do_cons_clim_nudg
Definition REMORA_DataStruct.H:1082