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 const int nrhs = 0;
60 int nnew = 0;
61
62 int iic = istep[lev];
63 int ntfirst = 0;
64
65 // Because zeta may have changed
66 stretch_transform(lev);
67
68 // These temporaries used to be made in advance_3d_ml and passed in;
69 // now we make them here
70
71 const BoxArray& ba = mf_cons.boxArray();
72 const DistributionMapping& dm = mf_cons.DistributionMap();
73
74 //Only used locally, probably should be rearranged into FArrayBox declaration
75 MultiFab mf_AK (convert(ba,IntVect(0,0,1)),dm,1,IntVect(NGROW,NGROW,0)); //2d missing j coordinate
76 MultiFab mf_DC (ba,dm,1,IntVect(NGROW,NGROW,NGROW-1)); //2d missing j coordinate
77 MultiFab mf_Hzk(ba,dm,1,IntVect(NGROW,NGROW,NGROW-1)); //2d missing j coordinate
78
79 for ( MFIter mfi(mf_cons, TilingIfNotGPU()); mfi.isValid(); ++mfi )
80 {
81 Array4<Real > const& u = mf_u.array(mfi);
82 Array4<Real > const& v = mf_v.array(mfi);
83
84 Array4<Real > const& ru = mf_ru->array(mfi);
85 Array4<Real > const& rv = mf_rv->array(mfi);
86
87 Array4<Real > const& AK = mf_AK.array(mfi);
88 Array4<Real > const& DC = mf_DC.array(mfi);
89
90 Array4<Real > const& Hzk = mf_Hzk.array(mfi);
91 Array4<Real > const& Akv = mf_Akv->array(mfi);
92
93 Array4<Real const> const& Hz = mf_Hz->const_array(mfi);
94
95 Array4<Real const> const& DU_avg1 = mf_DU_avg1->const_array(mfi);
96 Array4<Real const> const& DV_avg1 = mf_DV_avg1->const_array(mfi);
97
98 Array4<Real const> const& pm = mf_pm->const_array(mfi);
99 Array4<Real const> const& pn = mf_pn->const_array(mfi);
100
101 Array4<Real const> const& msku = mf_msku->const_array(mfi);
102 Array4<Real const> const& mskv = mf_mskv->const_array(mfi);
103
104 Box bx = mfi.tilebox();
105 Box gbx2 = mfi.growntilebox(IntVect(NGROW,NGROW,0));
106 Box gbx21 = mfi.growntilebox(IntVect(NGROW,NGROW,NGROW-1));
107
108 Box xbx = mfi.nodaltilebox(0);
109 Box ybx = mfi.nodaltilebox(1);
110
111 Box gbx2D = gbx2;
112 gbx2D.makeSlab(2,0);
113
114 Box tbxp1 = bx;
115 Box tbxp11 = bx;
116 Box tbxp2 = bx;
117 tbxp1.grow(IntVect(NGROW-1,NGROW-1,0));
118 tbxp2.grow(IntVect(NGROW,NGROW,0));
119 tbxp11.grow(IntVect(NGROW-1,NGROW-1,NGROW-1));
120
121 FArrayBox fab_FC(convert(gbx2,IntVect(0,0,1)),1,amrex::The_Async_Arena());
122 FArrayBox fab_BC(convert(gbx2,IntVect(0,0,1)),1,amrex::The_Async_Arena());
123 FArrayBox fab_CF(gbx21,1,amrex::The_Async_Arena());
124 FArrayBox fab_W(tbxp2,1,amrex::The_Async_Arena());
125
126 auto FC = fab_FC.array();
127 auto BC = fab_BC.array();
128 auto CF = fab_CF.array();
129
130 Real cff;
131 if (iic==ntfirst) {
132 cff=0.25_rt*dt_lev;
133 } else if (iic==ntfirst+1) {
134 cff=0.25_rt*dt_lev*3.0_rt/2.0_rt;
135 } else {
136 cff=0.25_rt*dt_lev*23.0_rt/12.0_rt;
137 }
138
139 ParallelFor(xbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
140 {
141 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);
142 u(i,j,k) *= 2.0_rt / (Hz(i-1,j,k) + Hz(i,j,k));
143 });
144
145 ParallelFor(ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
146 {
147 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);
148 v(i,j,k) *= 2.0_rt / (Hz(i,j-1,k) + Hz(i,j,k));
149 });
150
151 // NOTE: DC is only used as scratch in vert_visc_3d -- no need to pass or return a value
152 // NOTE: may not actually need to set these to zero
153
154 // Reset to zero on the box on which they'll be used
155 mf_DC[mfi].template setVal<RunOn::Device>(0.,xbx);
156 fab_CF.template setVal<RunOn::Device>(0.,xbx);
157
158 vert_visc_3d(xbx,1,0,u,Hz,Hzk,AK,Akv,BC,DC,FC,CF,nnew,N,dt_lev);
159
160 // Reset to zero on the box on which they'll be used
161 mf_DC[mfi].template setVal<RunOn::Device>(0.,ybx);
162 fab_CF.template setVal<RunOn::Device>(0.,ybx);
163
164 vert_visc_3d(ybx,0,1,v,Hz,Hzk,AK,Akv,BC,DC,FC,CF,nnew,N,dt_lev);
165
166 // Reset to zero on the box on which they'll be used
167 mf_DC[mfi].template setVal<RunOn::Device>(0.,xbx);
168 fab_CF.template setVal<RunOn::Device>(0.,xbx);
169
170 vert_mean_3d(xbx,1,0,u,Hz,DU_avg1,DC,CF,pn,msku,nnew,N);
171
172 // Reset to zero on the box on which they'll be used
173 mf_DC[mfi].template setVal<RunOn::Device>(0.,ybx);
174 fab_CF.template setVal<RunOn::Device>(0.,ybx);
175
176 vert_mean_3d(ybx,0,1,v,Hz,DV_avg1,DC,CF,pm,mskv,nnew,N);
177 }
178
179 // Apply physical boundary conditions to u and v
180 (*physbcs[lev])(mf_u,*mf_msku,0,1,mf_u.nGrowVect(),t_old[lev],BCVars::xvel_bc,0,*xvel_old[lev]);
181 (*physbcs[lev])(mf_v,*mf_mskv,0,1,mf_v.nGrowVect(),t_old[lev],BCVars::yvel_bc,0,*yvel_old[lev]);
182
183#ifdef REMORA_USE_NETCDF
184 // Fill the data which is stored in the boundary data read from netcdf files
185 if ( (solverChoice.boundary_from_netcdf) && (lev==0) )
186 {
187 fill_from_bdyfiles(mf_u,*mf_msku,t_old[lev],BCVars::xvel_bc,BdyVars::u,0,0,*xvel_old[lev],dt_lev);
188 fill_from_bdyfiles(mf_v,*mf_mskv,t_old[lev],BCVars::yvel_bc,BdyVars::v,0,0,*yvel_old[lev],dt_lev);
189 }
190
191 if (solverChoice.do_rivers) {
192 river_source_transport->update_interpolated_to_time(t_old[lev]);
193 }
194#endif
195 for ( MFIter mfi(mf_cons, TilingIfNotGPU()); mfi.isValid(); ++mfi )
196 {
197 Array4<Real > const& u = mf_u.array(mfi);
198 Array4<Real > const& v = mf_v.array(mfi);
199
200 Array4<Real > const& DC = mf_DC.array(mfi);
201
202 Array4<Real const> const& Hz = mf_Hz->const_array(mfi);
203
204 Array4<Real const> const& DU_avg1 = mf_DU_avg1->const_array(mfi);
205 Array4<Real const> const& DV_avg1 = mf_DV_avg1->const_array(mfi);
206
207 Array4<Real const> const& DU_avg2 = mf_DU_avg2->const_array(mfi);
208 Array4<Real const> const& DV_avg2 = mf_DV_avg2->const_array(mfi);
209
210 Array4<Real> const& ubar = mf_ubar->array(mfi);
211 Array4<Real> const& vbar = mf_vbar->array(mfi);
212
213 Array4<Real> const& Huon = mf_Huon->array(mfi);
214 Array4<Real> const& Hvom = mf_Hvom->array(mfi);
215
216 Array4<Real const> const& pm = mf_pm->const_array(mfi);
217 Array4<Real const> const& pn = mf_pn->const_array(mfi);
218 Array4<Real const> const& msku = mf_msku->const_array(mfi);
219 Array4<Real const> const& mskv = mf_mskv->const_array(mfi);
220
221 Array4<Real const> const& z_w = mf_z_w->const_array(mfi);
222
223 Box gbx1 = mfi.growntilebox(IntVect(NGROW-1,NGROW-1,0));
224 Box gbx2 = mfi.growntilebox(IntVect(NGROW,NGROW,0));
225
226 FArrayBox fab_FC(gbx2,1,amrex::The_Async_Arena());
227 auto FC = fab_FC.array();
228
229#ifdef REMORA_USE_NETCDF
230 if (solverChoice.do_rivers) {
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 = 2.0_rt / (pn(i,j,0)+pn(i-1,j,0));
240 Real cff = 1.0_rt / (on_u * 0.5_rt * (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 = 2.0_rt / (pm(i,j,0)+pm(i,j-1,0));
244 Real cff = 1.0_rt / (om_v * 0.5_rt * (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>(0.,grow(xbx,IntVect(0,0,1)));
255 fab_CF.template setVal<RunOn::Device>(0.,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>(0.,grow(ybx,IntVect(0,0,1)));
261 fab_CF.template setVal<RunOn::Device>(0.,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>(0.,gbx2);
269 mf_DC[mfi].template setVal<RunOn::Device>(0.,grow(gbx2,IntVect(0,0,1)));
270 update_massflux_3d(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>(0.,gbx2);
274 mf_DC[mfi].template setVal<RunOn::Device>(0.,grow(gbx2,IntVect(0,0,1)));
275 update_massflux_3d(gbx2,0,1,v,vbar,Hvom,Hz,pm,DV_avg1,DV_avg2,DC,FC,mskv,nnew);
276#endif
277 }
278
279 // WE BELIEVE THESE VALUES SHOULD ALREADY BE FILLED
280 // mf_Huon->FillBoundary(geom[lev].periodicity());
281 // mf_Hvom->FillBoundary(geom[lev].periodicity());
282
283 // ************************************************************************
284 // This should fill both temp and salt with temp/salt currently in cons_old
285 // ************************************************************************
286
287 MultiFab mf_W(convert(ba,IntVect(0,0,1)),dm,1,IntVect(NGROW+1,NGROW+1,0));
288 mf_W.setVal(0.0_rt);
289 for ( MFIter mfi(mf_cons, TilingIfNotGPU()); mfi.isValid(); ++mfi )
290 {
291 Array4<Real> const& Huon = mf_Huon->array(mfi);
292 Array4<Real> const& Hvom = mf_Hvom->array(mfi);
293
294 Array4<Real const> const& z_w = mf_z_w->const_array(mfi);
295 Array4<Real const> const& h = mf_h->const_array(mfi);
296
297 Box bx = mfi.tilebox();
298 Box gbx1 = mfi.growntilebox(IntVect(NGROW-1,NGROW-1,0));
299 Box gbx2 = mfi.growntilebox(IntVect(NGROW,NGROW,0));
300 Box gbx21 = mfi.growntilebox(IntVect(NGROW,NGROW,NGROW-1));
301
302 Box tbxp1 = bx;
303 Box tbxp11 = bx;
304 Box tbxp2 = bx;
305 tbxp1.grow(IntVect(NGROW-1,NGROW-1,0));
306 tbxp2.grow(IntVect(NGROW,NGROW,0));
307 tbxp11.grow(IntVect(NGROW-1,NGROW-1,NGROW-1));
308
309 FArrayBox fab_FC(surroundingNodes(gbx2,2),1,amrex::The_Async_Arena());
310 FArrayBox fab_BC(gbx2,1,amrex::The_Async_Arena());
311 FArrayBox fab_CF(gbx21,1,amrex::The_Async_Arena());
312
313 auto W = mf_W.array(mfi);
314
315 //
316 //------------------------------------------------------------------------
317 // Vertically integrate horizontal mass flux divergence.
318 //------------------------------------------------------------------------
319 //
320 //Should really use gbx3uneven
321 //TODO: go over these boxes and compare to other spots where we do the same thing
322 Box gbx1D = gbx1;
323 gbx1D.makeSlab(2,0);
324
325 // Starting with zero vertical velocity at the bottom, integrate
326 // from the bottom (k=0) to the free-surface (k=N). The w(:,:,N(ng))
327 // contains the vertical velocity at the free-surface, d(zeta)/d(t).
328 // Notice that barotropic mass flux divergence is not used directly.
329 ParallelFor(gbx1D, [=] AMREX_GPU_DEVICE (int i, int j, int )
330 {
331 W(i,j,0) = 0.0_rt;
332 for (int k=1; k<=N+1; k++) {
333 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));
334 }
335 });
336
337 // Starting with zero vertical velocity at the bottom, integrate
338 // from the bottom (k=0) to the free-surface (k=N). The w(:,:,N(ng))
339 // contains the vertical velocity at the free-surface, d(zeta)/d(t).
340 // Notice that barotropic mass flux divergence is not used directly.
341 //
342 ParallelFor(gbx1D, [=] AMREX_GPU_DEVICE (int i, int j, int )
343 {
344 Real wrk_ij = W(i,j,N+1) / (z_w(i,j,N+1)+h(i,j,0,0));
345
346 for (int k=1; k<=N; k++) {
347 W(i,j,k) -= wrk_ij * (z_w(i,j,k)+h(i,j,0,0));
348 }
349 W(i,j,N+1) = 0.0_rt;
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 nnew = 0;
362
363 for ( MFIter mfi(mf_cons, TilingIfNotGPU()); mfi.isValid(); ++mfi )
364 {
365 Array4<Real> const& Hz = mf_Hz->array(mfi);
366
367 Array4<Real> const& Huon = mf_Huon->array(mfi);
368 Array4<Real> const& Hvom = mf_Hvom->array(mfi);
369
370 Array4<Real const> const& pm = mf_pm->const_array(mfi);
371 Array4<Real const> const& pn = mf_pn->const_array(mfi);
372 Array4<Real const> const& mskr = mf_mskr->const_array(mfi);
373 Array4<Real const> const& msku = mf_msku->const_array(mfi);
374 Array4<Real const> const& mskv = mf_mskv->const_array(mfi);
375
376 Box bx = mfi.tilebox();
377 Box gbx2 = mfi.growntilebox(IntVect(NGROW,NGROW,0));
378 Box gbx21 = mfi.growntilebox(IntVect(NGROW,NGROW,NGROW-1));
379
380 Box tbxp1 = bx;
381 Box tbxp11 = bx;
382 Box tbxp2 = bx;
383 tbxp1.grow(IntVect(NGROW-1,NGROW-1,0));
384 tbxp2.grow(IntVect(NGROW,NGROW,0));
385 tbxp11.grow(IntVect(NGROW-1,NGROW-1,NGROW-1));
386
387 FArrayBox fab_FC(surroundingNodes(gbx2,2),1,amrex::The_Async_Arena());
388 FArrayBox fab_BC(gbx2,1,amrex::The_Async_Arena());
389 FArrayBox fab_CF(gbx21,1,amrex::The_Async_Arena());
390
391 auto FC = fab_FC.array();
392 auto W = mf_W.array(mfi);
393 //
394 //-----------------------------------------------------------------------
395 // rhs_t_3d
396 //-----------------------------------------------------------------------
397 //
398 for (int i_comp=0; i_comp < NCONS; i_comp++)
399 {
400#ifdef REMORA_USE_NETCDF
401 FArrayBox* fab_river_source;
402 if (solverChoice.do_rivers_cons[i_comp]) {
403 fab_river_source = river_source_cons[i_comp]->fab_interp;
404 }
405 const Array4<const int>& river_pos = (solverChoice.do_rivers) ? vec_river_position[lev]->const_array(mfi) : Array4<const int>();
406 const Array4<const Real>& river_source = (solverChoice.do_rivers_cons[i_comp]) ? fab_river_source->const_array() : Array4<const Real>();
407#else
408 const Array4<const int >& river_pos = Array4<const int>();
409 const Array4<const Real>& river_source = Array4<const Real>();
410#endif
411 Array4<Real> const& sstore = mf_sstore->array(mfi, i_comp);
412 rhs_t_3d(bx, mf_cons.array(mfi,i_comp), sstore, Huon, Hvom,
413 Hz, pn, pm, W, FC, mskr, msku, mskv, river_pos, river_source, nrhs, nnew, N,dt_lev);
414 }
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,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
458#ifdef REMORA_USE_NETCDF
459 if (solverChoice.do_temp_clim_nudg) {
460 temp_clim_data_from_file->update_interpolated_to_time(t_old[lev]);
461
462 for ( MFIter mfi(mf_cons, TilingIfNotGPU()); mfi.isValid(); ++mfi )
463 {
464 Box bx = mfi.growntilebox(IntVect(1,1,0));
465 Array4<const Real> const& temp_nudg_coeff = vec_nudg_coeff[BdyVars::t][lev]->const_array(mfi);
466 Array4<const Real> const& temp_clim = temp_clim_data_from_file->mf_interpolated->const_array(mfi);
467 Array4< Real> const& temp = cons_new[lev]->array(mfi, Temp_comp);
468 Array4<const Real> const& Hz = vec_Hz[lev]->const_array(mfi);
469 Array4<const Real> const& pm = vec_pm[lev]->const_array(mfi);
470 Array4<const Real> const& pn = vec_pn[lev]->const_array(mfi);
471
472 apply_clim_nudg(bx, 0, 0, temp, temp, temp_clim, temp_nudg_coeff, Hz, pm, pn, dt_lev);
473 }
474 }
475 if (solverChoice.do_salt_clim_nudg) {
476 salt_clim_data_from_file->update_interpolated_to_time(t_old[lev]);
477
478 for ( MFIter mfi(mf_cons, TilingIfNotGPU()); mfi.isValid(); ++mfi )
479 {
480 Box bx = mfi.growntilebox(IntVect(1,1,0));
481 Array4<const Real> const& salt_nudg_coeff = vec_nudg_coeff[BdyVars::s][lev]->const_array(mfi);
482 Array4<const Real> const& salt_clim = salt_clim_data_from_file->mf_interpolated->const_array(mfi);
483 Array4< Real> const& salt = cons_new[lev]->array(mfi, Salt_comp);
484 Array4<const Real> const& Hz = vec_Hz[lev]->const_array(mfi);
485 Array4<const Real> const& pm = vec_pm[lev]->const_array(mfi);
486 Array4<const Real> const& pn = vec_pn[lev]->const_array(mfi);
487
488 apply_clim_nudg(bx, 0, 0, salt, salt, salt_clim, salt_nudg_coeff, Hz, pm, pn, dt_lev);
489 }
490 }
491#endif
492}
VertMixingType::GLS
@ GLS
NGROW
#define NGROW
Definition REMORA_IndexDefines.H:13
Temp_comp
#define Temp_comp
Definition REMORA_IndexDefines.H:8
Salt_comp
#define Salt_comp
Definition REMORA_IndexDefines.H:9
NCONS
#define NCONS
Definition REMORA_IndexDefines.H:11
NCTimeSeriesRiver::update_interpolated_to_time
void update_interpolated_to_time(amrex::Real time)
Calculate.
Definition REMORA_NCTimeSeriesRiver.cpp:109
NCTimeSeriesRiver::fab_interp
amrex::FArrayBox * fab_interp
Container for interpolated data; Only used if save_interpolated == true.
Definition REMORA_NCTimeSeriesRiver.H:26
NCTimeSeries::mf_interpolated
amrex::MultiFab * mf_interpolated
Definition REMORA_NCTimeSeries.H:35
NCTimeSeries::update_interpolated_to_time
void update_interpolated_to_time(amrex::Real time)
Calculate interpolated values at time, reading in data as necessary.
Definition REMORA_NCTimeSeries.cpp:96
REMORA::river_source_cons
amrex::Vector< NCTimeSeriesRiver * > river_source_cons
Vector of data containers for scalar data in rivers.
Definition REMORA.H:1045
REMORA::vec_pm
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_pm
horizontal scaling factor: 1 / dx (2D)
Definition REMORA.H:358
REMORA::cons_new
amrex::Vector< amrex::MultiFab * > cons_new
multilevel data container for current step's scalar data: temperature, salinity, passive scalar
Definition REMORA.H:220
REMORA::stretch_transform
void stretch_transform(int lev)
Calculate vertical stretched coordinates.
Definition REMORA_DepthStretchTransform.H:15
REMORA::vec_mskr
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_mskr
land/sea mask at cell centers (2D)
Definition REMORA.H:349
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::rhs_t_3d
void rhs_t_3d(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:27
REMORA::vec_tke
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_tke
Turbulent kinetic energy.
Definition REMORA.H:406
REMORA::vec_gls
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_gls
Turbulent generic length scale.
Definition REMORA.H:408
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:235
REMORA::vec_Akt
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Akt
Vertical diffusion coefficient (3D)
Definition REMORA.H:255
REMORA::vec_msku
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_msku
land/sea mask at x-faces (2D)
Definition REMORA.H:351
REMORA::xvel_old
amrex::Vector< amrex::MultiFab * > xvel_old
multilevel data container for last step's x velocities (u in ROMS)
Definition REMORA.H:213
REMORA::vec_mskv
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_mskv
land/sea mask at y-faces (2D)
Definition REMORA.H:353
REMORA::physbcs
amrex::Vector< std::unique_ptr< REMORAPhysBCFunct > > physbcs
Vector (over level) of functors to apply physical boundary conditions.
Definition REMORA.H:1124
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=true, const int n_not_fill=0, const int icomp_calc=0, const amrex::Real dt=amrex::Real(0.0), 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::istep
amrex::Vector< int > istep
which step?
Definition REMORA.H:1110
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=amrex::Real(0.0))
Apply climatology nudging.
Definition REMORA_apply_clim_nudg.cpp:19
REMORA::fill_from_bdyfiles
void fill_from_bdyfiles(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=amrex::Real(0.0))
Fill boundary data from netcdf file.
Definition REMORA_BoundaryConditions_netcdf.cpp:19
REMORA::update_massflux_3d
void update_massflux_3d(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:22
REMORA::yvel_old
amrex::Vector< amrex::MultiFab * > yvel_old
multilevel data container for last step's y velocities (v in ROMS)
Definition REMORA.H:215
REMORA::river_source_transport
NCTimeSeriesRiver * river_source_transport
Data container for momentum transport in rivers.
Definition REMORA.H:1047
REMORA::temp_clim_data_from_file
NCTimeSeries * temp_clim_data_from_file
Data container for temperature climatology data read from file.
Definition REMORA.H:1040
REMORA::solverChoice
static SolverChoice solverChoice
Container for algorithmic choices.
Definition REMORA.H:1237
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:1052
REMORA::vec_Akk
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Akk
Turbulent kinetic energy vertical diffusion coefficient.
Definition REMORA.H:412
REMORA::salt_clim_data_from_file
NCTimeSeries * salt_clim_data_from_file
Data container for salinity climatology data read from file.
Definition REMORA.H:1042
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::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:1054
REMORA::cons_old
amrex::Vector< amrex::MultiFab * > cons_old
multilevel data container for last step's scalar data: temperature, salinity, passive scalar
Definition REMORA.H:211
REMORA::vec_nudg_coeff
amrex::Vector< amrex::Vector< std::unique_ptr< amrex::MultiFab > > > vec_nudg_coeff
Climatology nudging coefficients.
Definition REMORA.H:417
REMORA::vec_pn
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_pn
horizontal scaling factor: 1 / dy (2D)
Definition REMORA.H:360
REMORA::vec_Akv
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Akv
Vertical viscosity coefficient (3D)
Definition REMORA.H:253
REMORA::t_old
amrex::Vector< amrex::Real > t_old
old time at each level
Definition REMORA.H:1116
REMORA::vec_Akp
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Akp
Turbulent length scale vertical diffusion coefficient.
Definition REMORA.H:414
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:246
SolverChoice::do_rivers_cons
amrex::Vector< int > do_rivers_cons
Definition REMORA_DataStruct.H:560
SolverChoice::vert_mixing_type
VertMixingType vert_mixing_type
Definition REMORA_DataStruct.H:599