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=0.25_rt*dt_lev;
134 } else if (iic==ntfirst+1) {
135 cff=0.25_rt*dt_lev*3.0_rt/2.0_rt;
136 } else {
137 cff=0.25_rt*dt_lev*23.0_rt/12.0_rt;
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) *= 2.0_rt / (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) *= 2.0_rt / (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>(0.,xbx);
157 fab_CF.template setVal<RunOn::Device>(0.,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>(0.,ybx);
163 fab_CF.template setVal<RunOn::Device>(0.,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>(0.,xbx);
169 fab_CF.template setVal<RunOn::Device>(0.,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>(0.,ybx);
175 fab_CF.template setVal<RunOn::Device>(0.,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],BCVars::xvel_bc,0,*xvel_old[lev]);
182 (*physbcs[lev])(mf_v,*mf_mskv,0,1,mf_v.nGrowVect(),t_old[lev],BCVars::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) && (lev==0) )
187 {
188 fill_from_bdyfiles(mf_u,*mf_msku,t_old[lev],BCVars::xvel_bc,BdyVars::u,0,0,*xvel_old[lev],dt_lev);
189 fill_from_bdyfiles(mf_v,*mf_mskv,t_old[lev],BCVars::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 = 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(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>(0.,gbx2);
274 mf_DC[mfi].template setVal<RunOn::Device>(0.,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#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 ParallelFor(gbx1D, N+1,
326 [=] AMREX_GPU_DEVICE (int i, int j, int , int kk)
327 {
328 // Starting with zero vertical velocity at the bottom, integrate
329 // from the bottom (k=0) to the free-surface (k=N). The w(:,:,N(ng))
330 // contains the vertical velocity at the free-surface, d(zeta)/d(t).
331 // Notice that barotropic mass flux divergence is not used directly.
332 //
333 int k = kk + 1;
334 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));
335 });
336 ParallelFor(gbx1D, [=] AMREX_GPU_DEVICE (int i, int j, int )
337 {
338 W(i,j,N+1)=W(i,j,N+1)/(z_w(i,j,N+1)+h(i,j,0,0)); // wrk_i
339 });
340 ParallelFor(gbx1D, N,
341 [=] AMREX_GPU_DEVICE (int i, int j, int , int kk)
342 {
343 int k = kk + 1;
344 W(i,j,k) = W(i,j,k)- W(i,j,N+1)*(z_w(i,j,k)+h(i,j,0,0));
345 });
346 ParallelFor(gbx1D, [=] AMREX_GPU_DEVICE (int i, int j, int )
347 {
348 W(i,j,N+1) = 0.0_rt;
349 });
350 }
351
352 const int nstp = (iic) % 2;
353 nnew = 1-nstp;
354 if (solverChoice.vert_mixing_type == VertMixingType::GLS) {
355 gls_corrector(lev, vec_gls[lev].get(), vec_tke[lev].get(), mf_W, vec_Akv[lev].get(),
356 vec_Akt[lev].get(),vec_Akk[lev].get(), vec_Akp[lev].get(), vec_mskr[lev].get(),
357 vec_msku[lev].get(), vec_mskv[lev].get(),
358 nstp, nnew, N, dt_lev);
359 }
360 nnew = 0;
361
362 for ( MFIter mfi(mf_cons, TilingIfNotGPU()); mfi.isValid(); ++mfi )
363 {
364 Array4<Real> const& Hz = mf_Hz->array(mfi);
365
366 Array4<Real> const& Huon = mf_Huon->array(mfi);
367 Array4<Real> const& Hvom = mf_Hvom->array(mfi);
368
369 Array4<Real const> const& pm = mf_pm->const_array(mfi);
370 Array4<Real const> const& pn = mf_pn->const_array(mfi);
371 Array4<Real const> const& mskr = mf_mskr->const_array(mfi);
372 Array4<Real const> const& msku = mf_msku->const_array(mfi);
373 Array4<Real const> const& mskv = mf_mskv->const_array(mfi);
374
375 Box bx = mfi.tilebox();
376 Box gbx2 = mfi.growntilebox(IntVect(NGROW,NGROW,0));
377 Box gbx21 = mfi.growntilebox(IntVect(NGROW,NGROW,NGROW-1));
378
379 Box tbxp1 = bx;
380 Box tbxp11 = bx;
381 Box tbxp2 = bx;
382 tbxp1.grow(IntVect(NGROW-1,NGROW-1,0));
383 tbxp2.grow(IntVect(NGROW,NGROW,0));
384 tbxp11.grow(IntVect(NGROW-1,NGROW-1,NGROW-1));
385
386 FArrayBox fab_FC(surroundingNodes(gbx2,2),1,amrex::The_Async_Arena());
387 FArrayBox fab_BC(gbx2,1,amrex::The_Async_Arena());
388 FArrayBox fab_CF(gbx21,1,amrex::The_Async_Arena());
389
390 auto FC = fab_FC.array();
391 auto W = mf_W.array(mfi);
392 //
393 //-----------------------------------------------------------------------
394 // rhs_t_3d
395 //-----------------------------------------------------------------------
396 //
397 for (int i_comp=0; i_comp < NCONS; i_comp++)
398 {
399#ifdef REMORA_USE_NETCDF
400 FArrayBox* fab_river_source;
401 if (solverChoice.do_rivers_cons[i_comp]) {
402 fab_river_source = river_source_cons[i_comp]->fab_interp;
403 }
404 const Array4<const int>& river_pos = (solverChoice.do_rivers) ? vec_river_position[lev]->const_array(mfi) : Array4<const int>();
405 const Array4<const Real>& river_source = (solverChoice.do_rivers_cons[i_comp]) ? fab_river_source->const_array() : Array4<const Real>();
406#else
407 const Array4<const int >& river_pos = Array4<const int>();
408 const Array4<const Real>& river_source = Array4<const Real>();
409#endif
410 Array4<Real> const& sstore = mf_sstore->array(mfi, i_comp);
411 rhs_t_3d(lev,bx, mf_cons.array(mfi,i_comp), sstore, Huon, Hvom,
412 Hz, pn, pm, W, FC, mskr, msku, mskv, river_pos, river_source, nrhs, nnew, N,dt_lev);
413 }
414 } // mfi
415
416 FillPatch(lev, t_old[lev], mf_cons, cons_new, BCVars::cons_bc, BdyVars::t,0,true,false,0,0,dt_lev,*cons_old[lev]);
417
418 for ( MFIter mfi(mf_cons, TilingIfNotGPU()); mfi.isValid(); ++mfi )
419 {
420 Array4<Real> const& AK = mf_AK.array(mfi);
421 Array4<Real> const& DC = mf_DC.array(mfi);
422
423 Array4<Real> const& Hzk = mf_Hzk.array(mfi);
424 Array4<Real const> const& Hz = mf_Hz->const_array(mfi);
425
426 Box bx = mfi.tilebox();
427
428 // Copy the tilebox
429 Box tbxp1 = bx;
430 Box tbxp11 = bx;
431 Box tbxp2 = bx;
432 Box tbxp21 = bx;
433 //make only gbx be grown to match multifabs
434 tbxp21.grow(IntVect(NGROW,NGROW,NGROW-1));
435 tbxp2.grow(IntVect(NGROW,NGROW,0));
436 tbxp1.grow(IntVect(NGROW-1,NGROW-1,0));
437 tbxp11.grow(IntVect(NGROW-1,NGROW-1,NGROW-1));
438
439 FArrayBox fab_FC(tbxp2,1,amrex::The_Async_Arena());
440 FArrayBox fab_BC(tbxp2,1,amrex::The_Async_Arena());
441 FArrayBox fab_CF(tbxp21,1,amrex::The_Async_Arena());
442 FArrayBox fab_W(tbxp2,1,amrex::The_Async_Arena());
443
444 auto FC = fab_FC.array();
445 auto BC = fab_BC.array();
446 auto CF = fab_CF.array();
447
448 for (int i_comp=0; i_comp < NCONS; i_comp++) {
449 vert_visc_3d(bx,0,0,mf_cons.array(mfi,i_comp),Hz,Hzk,
450 AK,mf_Akt->array(mfi,i_comp),BC,DC,FC,CF,nnew,N,dt_lev);
451 }
452 } // MFiter
453 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]);
454
455#ifdef REMORA_USE_NETCDF
456 if (solverChoice.do_temp_clim_nudg) {
457 temp_clim_data_from_file->update_interpolated_to_time(t_old[lev]);
458
459 for ( MFIter mfi(mf_cons, TilingIfNotGPU()); mfi.isValid(); ++mfi )
460 {
461 Box bx = mfi.growntilebox(IntVect(1,1,0));
462 Array4<const Real> const& temp_nudg_coeff = vec_nudg_coeff[BdyVars::t][lev]->const_array(mfi);
463 Array4<const Real> const& temp_clim = temp_clim_data_from_file->mf_interpolated->const_array(mfi);
464 Array4< Real> const& temp = cons_new[lev]->array(mfi, Temp_comp);
465 Array4<const Real> const& Hz = vec_Hz[lev]->const_array(mfi);
466 Array4<const Real> const& pm = vec_pm[lev]->const_array(mfi);
467 Array4<const Real> const& pn = vec_pn[lev]->const_array(mfi);
468
469 apply_clim_nudg(bx, 0, 0, temp, temp, temp_clim, temp_nudg_coeff, Hz, pm, pn, dt_lev);
470 }
471 }
472 if (solverChoice.do_salt_clim_nudg) {
473 salt_clim_data_from_file->update_interpolated_to_time(t_old[lev]);
474
475 for ( MFIter mfi(mf_cons, TilingIfNotGPU()); mfi.isValid(); ++mfi )
476 {
477 Box bx = mfi.growntilebox(IntVect(1,1,0));
478 Array4<const Real> const& salt_nudg_coeff = vec_nudg_coeff[BdyVars::s][lev]->const_array(mfi);
479 Array4<const Real> const& salt_clim = salt_clim_data_from_file->mf_interpolated->const_array(mfi);
480 Array4< Real> const& salt = cons_new[lev]->array(mfi, Salt_comp);
481 Array4<const Real> const& Hz = vec_Hz[lev]->const_array(mfi);
482 Array4<const Real> const& pm = vec_pm[lev]->const_array(mfi);
483 Array4<const Real> const& pn = vec_pn[lev]->const_array(mfi);
484
485 apply_clim_nudg(bx, 0, 0, salt, salt, salt_clim, salt_nudg_coeff, Hz, pm, pn, dt_lev);
486 }
487 }
488#endif
489}
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::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::river_source_cons
amrex::Vector< NCTimeSeriesRiver * > river_source_cons
Vector of data containers for scalar data in rivers.
Definition REMORA.H:1087
REMORA::vec_pm
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_pm
horizontal scaling factor: 1 / dx (2D)
Definition REMORA.H:366
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:223
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:357
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:414
REMORA::vec_gls
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_gls
Turbulent generic length scale.
Definition REMORA.H:416
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:238
REMORA::vec_Akt
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Akt
Vertical diffusion coefficient (3D)
Definition REMORA.H:258
REMORA::vec_msku
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_msku
land/sea mask at x-faces (2D)
Definition REMORA.H:359
REMORA::xvel_old
amrex::Vector< amrex::MultiFab * > xvel_old
multilevel data container for last step's x velocities (u in ROMS)
Definition REMORA.H:216
REMORA::vec_mskv
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_mskv
land/sea mask at y-faces (2D)
Definition REMORA.H:361
REMORA::physbcs
amrex::Vector< std::unique_ptr< REMORAPhysBCFunct > > physbcs
Vector (over level) of functors to apply physical boundary conditions.
Definition REMORA.H:1166
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:1152
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::yvel_old
amrex::Vector< amrex::MultiFab * > yvel_old
multilevel data container for last step's y velocities (v in ROMS)
Definition REMORA.H:218
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::river_source_transport
NCTimeSeriesRiver * river_source_transport
Data container for momentum transport in rivers.
Definition REMORA.H:1089
REMORA::temp_clim_data_from_file
NCTimeSeries * temp_clim_data_from_file
Data container for temperature climatology data read from file.
Definition REMORA.H:1082
REMORA::solverChoice
static SolverChoice solverChoice
Container for algorithmic choices.
Definition REMORA.H:1279
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:1094
REMORA::vec_Akk
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Akk
Turbulent kinetic energy vertical diffusion coefficient.
Definition REMORA.H:420
REMORA::salt_clim_data_from_file
NCTimeSeries * salt_clim_data_from_file
Data container for salinity climatology data read from file.
Definition REMORA.H:1084
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:1096
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:214
REMORA::vec_nudg_coeff
amrex::Vector< amrex::Vector< std::unique_ptr< amrex::MultiFab > > > vec_nudg_coeff
Climatology nudging coefficients.
Definition REMORA.H:425
REMORA::vec_pn
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_pn
horizontal scaling factor: 1 / dy (2D)
Definition REMORA.H:368
REMORA::vec_Akv
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Akv
Vertical viscosity coefficient (3D)
Definition REMORA.H:256
REMORA::t_old
amrex::Vector< amrex::Real > t_old
old time at each level
Definition REMORA.H:1158
REMORA::vec_Akp
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Akp
Turbulent length scale vertical diffusion coefficient.
Definition REMORA.H:422
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
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