REMORA
Regional Modeling of Oceans Refined Adaptively
Loading...
Searching...
No Matches
REMORA_init_from_netcdf.cpp
Go to the documentation of this file.
1/**
2 * \file REMORA_init_from_netcdf.cpp
3 */
4
5#include <REMORA.H>
6#include <REMORA_Constants.H>
8#include <REMORA_DataStruct.H>
9
10using namespace amrex;
11
12#ifdef REMORA_USE_NETCDF
13
14/** \brief helper function for reading in initial state data from netcdf */
15void
16read_data_from_netcdf (int /*lev*/, const Box& domain, const std::string& fname,
17 FArrayBox& NC_temp_fab, FArrayBox& NC_salt_fab,
18 FArrayBox& NC_xvel_fab, FArrayBox& NC_yvel_fab,
19 FArrayBox& NC_ubar_fab, FArrayBox& NC_vbar_fab);
20
21/** \brief helper function for reading in land-sea masks from netcdf */
22void
23read_masks_from_netcdf (int /*lev*/, const Box& domain, const std::string& fname,
24 FArrayBox& NC_mskr_fab, FArrayBox& NC_msku_fab,
25 FArrayBox& NC_mskv_fab);
26
27/** \brief helper function for reading boundary data from netcdf */
28Real
29read_bdry_from_netcdf (const Box& domain, const std::string& fname,
30 Vector<Vector<FArrayBox>>& bdy_data_xlo,
31 Vector<Vector<FArrayBox>>& bdy_data_xhi,
32 Vector<Vector<FArrayBox>>& bdy_data_ylo,
33 Vector<Vector<FArrayBox>>& bdy_data_yhi,
34 int& width, amrex::Real& start_bdy_time,
35 std::string bdry_time_varname,
36 amrex::GpuArray<amrex::GpuArray<bool, AMREX_SPACEDIM*2>,BdyVars::NumTypes+1>&);
37
38/** \brief helper function to initialize state from netcdf */
39void
41 FArrayBox& temp_fab, FArrayBox& salt_fab,
42 FArrayBox& x_vel_fab, FArrayBox& y_vel_fab,
43 FArrayBox& ubar_fab, FArrayBox& vbar_fab,
44 const Vector<FArrayBox>& NC_temp_fab,
45 const Vector<FArrayBox>& NC_salt_fab,
46 const Vector<FArrayBox>& NC_xvel_fab,
47 const Vector<FArrayBox>& NC_yvel_fab,
48 const Vector<FArrayBox>& NC_ubar_fab,
49 const Vector<FArrayBox>& NC_vbar_fab);
50
51/** \brief helper function to read bathymetry from netcdf */
52void
53read_bathymetry_from_netcdf (int lev, const Box& domain, const std::string& fname,
54 FArrayBox& NC_h_fab,
55 FArrayBox& NC_pm_fab, FArrayBox& NC_pn_fab,
56 FArrayBox& NC_xr_fab, FArrayBox& NC_yr_fab,
57 FArrayBox& NC_xu_fab, FArrayBox& NC_yu_fab,
58 FArrayBox& NC_xv_fab, FArrayBox& NC_yv_fab,
59 FArrayBox& NC_xp_fab, FArrayBox& NC_yp_fab);
60
61/** \brief helper function to read coriolis factor from netcdf */
62void
63read_coriolis_from_netcdf (int lev, const Box& domain, const std::string& fname, FArrayBox& NC_fcor_fab);
64
65/** \brief helper function to read sea surface height from netcdf */
66void
67read_zeta_from_netcdf (int lev, const Box& domain, const std::string& fname,
68 FArrayBox& NC_zeta_fab);
69
70/** \brief helper function to read climatology nudging from netcdf */
71void
72read_clim_nudg_coeff_from_netcdf (int lev, const Box& domain, const std::string& fname,
73 bool do_m2_clim_nudg,
74 bool do_m3_clim_nudg,
75 bool do_temp_clim_nudg,
76 bool do_salt_clim_nudg,
77 FArrayBox& NC_M2NC_fab,
78 FArrayBox& NC_M3NC_fab,
79 FArrayBox& NC_TempNC_fab,
80 FArrayBox& NC_SaltNC_fab);
81
82/** \brief helper function to read in vector of data from netcdf */
83//template <typename DType>
84void read_vec_from_netcdf (int lev, const std::string& fname, const std::string& field_name, amrex::Vector<int>& vec_dat);
85
86/**
87 * @param lev Integer specifying the current level
88 */
89void
91{
92 // *** FArrayBox's at this level for holding the INITIAL data
93 Vector<FArrayBox> NC_temp_fab ; NC_temp_fab.resize(num_boxes_at_level[lev]);
94 Vector<FArrayBox> NC_salt_fab ; NC_salt_fab.resize(num_boxes_at_level[lev]);
95 Vector<FArrayBox> NC_xvel_fab ; NC_xvel_fab.resize(num_boxes_at_level[lev]);
96 Vector<FArrayBox> NC_yvel_fab ; NC_yvel_fab.resize(num_boxes_at_level[lev]);
97 Vector<FArrayBox> NC_ubar_fab ; NC_ubar_fab.resize(num_boxes_at_level[lev]);
98 Vector<FArrayBox> NC_vbar_fab ; NC_vbar_fab.resize(num_boxes_at_level[lev]);
99
100 for (int idx = 0; idx < num_boxes_at_level[lev]; idx++)
101 {
102 read_data_from_netcdf(lev, boxes_at_level[lev][idx], nc_init_file[lev][idx],
103 NC_temp_fab[idx], NC_salt_fab[idx],
104 NC_xvel_fab[idx], NC_yvel_fab[idx],
105 NC_ubar_fab[idx], NC_vbar_fab[idx]);
106 }
107
108
109 MultiFab mf_temp(*cons_new[lev], make_alias, Temp_comp, 1);
110 MultiFab mf_salt(*cons_new[lev], make_alias, Salt_comp, 1);
111
112#ifdef _OPENMP
113#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
114#endif
115 {
116 // Don't tile this since we are operating on full FABs in this routine
117 for ( MFIter mfi(*cons_new[lev], false); mfi.isValid(); ++mfi )
118 {
119 // Define fabs for holding the initial data
120 FArrayBox &temp_fab = mf_temp[mfi];
121 FArrayBox &salt_fab = mf_salt[mfi];
122 FArrayBox &xvel_fab = (*xvel_new[lev])[mfi];
123 FArrayBox &yvel_fab = (*yvel_new[lev])[mfi];
124 FArrayBox &ubar_fab = (*vec_ubar[lev])[mfi];
125 FArrayBox &vbar_fab = (*vec_vbar[lev])[mfi];
126
127 init_state_from_netcdf(lev, temp_fab, salt_fab,
128 xvel_fab, yvel_fab,
129 ubar_fab, vbar_fab,
130 NC_temp_fab, NC_salt_fab,
131 NC_xvel_fab, NC_yvel_fab,
132 NC_ubar_fab, NC_vbar_fab);
133 } // mf
134 } // omp
135}
136
137/**
138 * @param lev Integer specifying the current level
139 */
140void
142{
143 // *** FArrayBox's at this level for holding the INITIAL data
144 Vector<FArrayBox> NC_zeta_fab ; NC_zeta_fab.resize(num_boxes_at_level[lev]);
145
146 for (int idx = 0; idx < num_boxes_at_level[lev]; idx++)
147 {
148 read_zeta_from_netcdf(lev,boxes_at_level[lev][idx], nc_init_file[lev][idx],
149 NC_zeta_fab[idx]);
150
151#ifdef _OPENMP
152#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
153#endif
154 {
155 // Don't tile this since we are operating on full FABs in this routine
156 for ( MFIter mfi(*cons_new[lev], false); mfi.isValid(); ++mfi )
157 {
158 FArrayBox &zeta_fab = (*vec_zeta[lev])[mfi];
159
160 //
161 // FArrayBox to FArrayBox copy does "copy on intersection"
162 // This only works here because we have broadcast the FArrayBox of data from the netcdf file to all ranks
163 //
164
165 zeta_fab.template copy<RunOn::Device>(NC_zeta_fab[idx],0,0,1);
166 } // mf
167 } // omp
168 } // idx
169 vec_zeta[lev]->FillBoundary(geom[lev].periodicity());
170 (*physbcs[lev])(*vec_zeta[lev],*vec_mskr[lev].get(),0,1,vec_zeta[lev]->nGrowVect(),t_new[lev],BCVars::zeta_bc);
171// (*physbcs[lev])(*vec_zeta[lev],*vec_mskr[lev].get(),1,1,vec_zeta[lev]->nGrowVect(),t_new[lev],BCVars::zeta_bc);
172// (*physbcs[lev])(*vec_zeta[lev],*vec_mskr[lev].get(),2,1,vec_zeta[lev]->nGrowVect(),t_new[lev],BCVars::zeta_bc);
173
175 Real told = t_new[lev];
177 }
178 if (lev>0) {
179 FillPatch(lev, t_old[lev], *vec_zeta[lev], GetVecOfPtrs(vec_zeta), BCVars::zeta_bc, BdyVars::zeta,
180 0, false,false);
181 }
182// fill_from_bdyfiles(*vec_zeta[lev], *vec_mskr[lev], told, BCVars::zeta_bc,BdyVars::zeta,1,1);
183// fill_from_bdyfiles(*vec_zeta[lev], *vec_mskr[lev], told, BCVars::zeta_bc,BdyVars::zeta,2,2);
184}
185/**
186 * @param lev Integer specifying the current level
187 */
188void
190{
191 // *** FArrayBox's at this level for holding the INITIAL data
192 Vector<FArrayBox> NC_h_fab ; NC_h_fab.resize(num_boxes_at_level[lev]);
193 Vector<FArrayBox> NC_pm_fab ; NC_pm_fab.resize(num_boxes_at_level[lev]);
194 Vector<FArrayBox> NC_pn_fab ; NC_pn_fab.resize(num_boxes_at_level[lev]);
195
196 Vector<FArrayBox> NC_xr_fab ; NC_xr_fab.resize(num_boxes_at_level[lev]);
197 Vector<FArrayBox> NC_yr_fab ; NC_yr_fab.resize(num_boxes_at_level[lev]);
198 Vector<FArrayBox> NC_xu_fab ; NC_xu_fab.resize(num_boxes_at_level[lev]);
199 Vector<FArrayBox> NC_yu_fab ; NC_yu_fab.resize(num_boxes_at_level[lev]);
200 Vector<FArrayBox> NC_xv_fab ; NC_xv_fab.resize(num_boxes_at_level[lev]);
201 Vector<FArrayBox> NC_yv_fab ; NC_yv_fab.resize(num_boxes_at_level[lev]);
202 Vector<FArrayBox> NC_xp_fab ; NC_xp_fab.resize(num_boxes_at_level[lev]);
203 Vector<FArrayBox> NC_yp_fab ; NC_yp_fab.resize(num_boxes_at_level[lev]);
204
205 for (int idx = 0; idx < num_boxes_at_level[lev]; idx++)
206 {
208 NC_h_fab[idx],
209 NC_pm_fab[idx], NC_pn_fab[idx],
210 NC_xr_fab[idx], NC_yr_fab[idx],
211 NC_xu_fab[idx], NC_yu_fab[idx],
212 NC_xv_fab[idx], NC_yv_fab[idx],
213 NC_xp_fab[idx], NC_yp_fab[idx]);
214
215#ifdef _OPENMP
216#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
217#endif
218 {
219 // Don't tile this since we are operating on full FABs in this routine
220 for ( MFIter mfi(*cons_new[lev], false); mfi.isValid(); ++mfi )
221 {
222 FArrayBox &h_fab = (*vec_hOfTheConfusingName[lev])[mfi];
223 FArrayBox &pm_fab = (*vec_pm[lev])[mfi];
224 FArrayBox &pn_fab = (*vec_pn[lev])[mfi];
225 FArrayBox &xr_fab = (*vec_xr[lev])[mfi];
226 FArrayBox &yr_fab = (*vec_yr[lev])[mfi];
227 FArrayBox &xu_fab = (*vec_xu[lev])[mfi];
228 FArrayBox &yu_fab = (*vec_yu[lev])[mfi];
229 FArrayBox &xv_fab = (*vec_xv[lev])[mfi];
230 FArrayBox &yv_fab = (*vec_yv[lev])[mfi];
231 FArrayBox &xp_fab = (*vec_xp[lev])[mfi];
232 FArrayBox &yp_fab = (*vec_yp[lev])[mfi];
233
234 //
235 // FArrayBox to FArrayBox copy does "copy on intersection"
236 // This only works here because we have broadcast the FArrayBox of data from the netcdf file to all ranks
237 //
238
239 // Copy into both components of h
240 h_fab.template copy<RunOn::Device>(NC_h_fab[idx],0,0,1);
241 h_fab.template copy<RunOn::Device>(NC_h_fab[idx],0,1,1);
242
243 pm_fab.template copy<RunOn::Device>(NC_pm_fab[idx]);
244 pn_fab.template copy<RunOn::Device>(NC_pn_fab[idx]);
245
246 xr_fab.template copy<RunOn::Device>(NC_xr_fab[idx]);
247 yr_fab.template copy<RunOn::Device>(NC_yr_fab[idx]);
248 xu_fab.template copy<RunOn::Device>(NC_xu_fab[idx]);
249 yu_fab.template copy<RunOn::Device>(NC_yu_fab[idx]);
250 xv_fab.template copy<RunOn::Device>(NC_xv_fab[idx]);
251 yv_fab.template copy<RunOn::Device>(NC_yv_fab[idx]);
252 xp_fab.template copy<RunOn::Device>(NC_xp_fab[idx]);
253 yp_fab.template copy<RunOn::Device>(NC_yp_fab[idx]);
254 } // mf
255 } // omp
256 } // idx
257
258 const double dummy_time = 0.0_rt;
259 // Unconditional foextrap will overwrite periodicity, but EnforcePeriodicity will
260 // be called on h afterwards
261 FillPatch(lev,dummy_time,*vec_hOfTheConfusingName[lev],GetVecOfPtrs(vec_hOfTheConfusingName),
263 BdyVars::null,0,false,true,1);
264 FillPatch(lev,dummy_time,*vec_hOfTheConfusingName[lev],GetVecOfPtrs(vec_hOfTheConfusingName),
266 BdyVars::null,1,false,true,1);
267
268 if (lev > 0) {
269 FillPatch(lev,dummy_time,*vec_pm[lev],GetVecOfPtrs(vec_pm),
271 BdyVars::null,0,false,true);
272 FillPatch(lev,dummy_time,*vec_pn[lev],GetVecOfPtrs(vec_pn),
274 BdyVars::null,0,false,true);
275 }
276
277
278 int ng = vec_pm[lev]->nGrow();
279
280 const auto& dom_lo = amrex::lbound(geom[lev].Domain());
281 const auto& dom_hi = amrex::ubound(geom[lev].Domain());
282
283 //
284 // We need values of pm and pn outside the domain so we fill
285 // them here with foextrap
286 //
287 // We first fill interior ghost cells because we will need to extrapolate
288 // from ghost cells inside the domain to ghost cells outside the domain
289 //
290 vec_pm[lev]->FillBoundary(geom[lev].periodicity());
291 vec_pn[lev]->FillBoundary(geom[lev].periodicity());
292
293 vec_xr[lev]->FillBoundary(geom[lev].periodicity());
294 vec_yr[lev]->FillBoundary(geom[lev].periodicity());
295 vec_xu[lev]->FillBoundary(geom[lev].periodicity());
296 vec_yu[lev]->FillBoundary(geom[lev].periodicity());
297 vec_xv[lev]->FillBoundary(geom[lev].periodicity());
298 vec_yv[lev]->FillBoundary(geom[lev].periodicity());
299 vec_xp[lev]->FillBoundary(geom[lev].periodicity());
300 vec_yp[lev]->FillBoundary(geom[lev].periodicity());
301
302 for ( MFIter mfi(*vec_pm[lev]); mfi.isValid(); ++mfi )
303 {
304 Box bx = mfi.tilebox();
305
306 auto pm_fab = vec_pm[lev]->array(mfi);
307 auto pn_fab = vec_pn[lev]->array(mfi);
308
309 Box gbx_lox = adjCellLo(bx,0,ng); gbx_lox.grow(1,ng); gbx_lox.setBig (0,dom_lo.x-2);
310 Box gbx_hix = adjCellHi(bx,0,ng); gbx_hix.grow(1,ng); gbx_hix.setSmall(0,dom_hi.x+2);
311 Box gbx_loy = adjCellLo(bx,1,ng); gbx_loy.grow(0,ng); gbx_loy.setBig (1,dom_lo.y-2);
312 Box gbx_hiy = adjCellHi(bx,1,ng); gbx_hiy.grow(0,ng); gbx_hiy.setSmall(1,dom_hi.y+2);
313
314 // if (gbx_lox.ok()) amrex::AllPrint() << "GBX_XLO " << gbx_lox << std::endl;
315 // if (gbx_hix.ok()) amrex::AllPrint() << "GBX_XHI " << gbx_hix << std::endl;
316 // if (gbx_loy.ok()) amrex::AllPrint() << "GBX_YLO " << gbx_loy << std::endl;
317 // if (gbx_hiy.ok()) amrex::AllPrint() << "GBX_YHI " << gbx_hiy << std::endl;
318
319 if (gbx_lox.ok()) {
320 ParallelFor(gbx_lox, [=] AMREX_GPU_DEVICE (int i, int j, int k)
321 {
322 pm_fab(i,j,k,0) = pm_fab(dom_lo.x-1,j,k,0);
323 pn_fab(i,j,k,0) = pn_fab(dom_lo.x-1,j,k,0);
324 });
325 }
326 if (gbx_hix.ok()) {
327 ParallelFor(gbx_hix, [=] AMREX_GPU_DEVICE (int i, int j, int k)
328 {
329 pm_fab(i,j,k,0) = pm_fab(dom_hi.x+1,j,k,0);
330 pn_fab(i,j,k,0) = pn_fab(dom_hi.x+1,j,k,0);
331 });
332 }
333 if (gbx_loy.ok()) {
334 ParallelFor(gbx_loy, [=] AMREX_GPU_DEVICE (int i, int j, int k)
335 {
336 pm_fab(i,j,k,0) = pm_fab(i,dom_lo.y-1,k,0);
337 pn_fab(i,j,k,0) = pn_fab(i,dom_lo.y-1,k,0);
338 });
339 }
340 if (gbx_hiy.ok()) {
341 ParallelFor(gbx_hiy, [=] AMREX_GPU_DEVICE (int i, int j, int k)
342 {
343 pm_fab(i,j,k,0) = pm_fab(i,dom_hi.y+1,k,0);
344 pn_fab(i,j,k,0) = pn_fab(i,dom_hi.y+1,k,0);
345 });
346 }
347 } // mfi
348
349 vec_hOfTheConfusingName[lev]->FillBoundary(geom[lev].periodicity());
350}
351
352/**
353 * @param lev Integer specifying the current level
354 */
355void
357{
358 // *** FArrayBox's at this level for holding the INITIAL data
359 Vector<FArrayBox> NC_fcor_fab ; NC_fcor_fab.resize(num_boxes_at_level[lev]);
360
361 for (int idx = 0; idx < num_boxes_at_level[lev]; idx++)
362 {
364 NC_fcor_fab[idx]);
365
366#ifdef _OPENMP
367#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
368#endif
369 {
370 // Don't tile this since we are operating on full FABs in this routine
371 for ( MFIter mfi(*cons_new[lev], false); mfi.isValid(); ++mfi )
372 {
373 FArrayBox &fcor_fab = (*vec_fcor[lev])[mfi];
374
375 //
376 // FArrayBox to FArrayBox copy does "copy on intersection"
377 // This only works here because we have broadcast the FArrayBox of data from the netcdf file to all ranks
378 //
379
380 fcor_fab.template copy<RunOn::Device>(NC_fcor_fab[idx]);
381 } // mf
382 } // omp
383 } // idx
384 vec_fcor[lev]->FillBoundary(geom[lev].periodicity());
385}
386
387/**
388 * @param lev Integer specifying the current level
389 */
390void
392{
393 // *** FArrayBox's at this level for holding the INITIAL data
394 Vector<FArrayBox> NC_mskr_fab ; NC_mskr_fab.resize(num_boxes_at_level[lev]);
395 Vector<FArrayBox> NC_msku_fab ; NC_msku_fab.resize(num_boxes_at_level[lev]);
396 Vector<FArrayBox> NC_mskv_fab ; NC_mskv_fab.resize(num_boxes_at_level[lev]);
397
398 for (int idx = 0; idx < num_boxes_at_level[lev]; idx++)
399 {
400 read_masks_from_netcdf(lev,boxes_at_level[lev][idx], nc_grid_file[lev][idx],
401 NC_mskr_fab[idx],NC_msku_fab[idx],
402 NC_mskv_fab[idx]);
403
404#ifdef _OPENMP
405#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
406#endif
407 {
408 // Don't tile this since we are operating on full FABs in this routine
409 for ( MFIter mfi(*cons_new[lev], false); mfi.isValid(); ++mfi )
410 {
411 FArrayBox &mskr_fab = (*vec_mskr[lev])[mfi];
412 FArrayBox &msku_fab = (*vec_msku[lev])[mfi];
413 FArrayBox &mskv_fab = (*vec_mskv[lev])[mfi];
414
415 //
416 // FArrayBox to FArrayBox copy does "copy on intersection"
417 // This only works here because we have broadcast the FArrayBox of data from the netcdf file to all ranks
418 //
419
420 mskr_fab.template copy<RunOn::Device>(NC_mskr_fab[idx]);
421 msku_fab.template copy<RunOn::Device>(NC_msku_fab[idx]);
422 mskv_fab.template copy<RunOn::Device>(NC_mskv_fab[idx]);
423 } // mf
424 } // omp
425 } // idx
426
427 update_mskp(lev);
428 vec_mskr[lev]->FillBoundary(geom[lev].periodicity());
429 vec_msku[lev]->FillBoundary(geom[lev].periodicity());
430 vec_mskv[lev]->FillBoundary(geom[lev].periodicity());
431 vec_mskp[lev]->FillBoundary(geom[lev].periodicity());
432}
433
434/**
435 * @param lev Integer specifying the current level
436 */
437void
439{
440 if (nc_bdry_file.empty()) {
441 amrex::Error("NetCDF boundary file name must be provided via input");
442 }
443
444 int width = 1;
445 for (int ifile = 0; ifile < nc_bdry_file.size(); ifile++) {
446 amrex::Real bdy_time_interval_this = read_bdry_from_netcdf(geom[0].Domain(), nc_bdry_file[ifile],
448 width, start_bdy_time,
451 if (ifile == 0) {
452 bdy_time_interval = bdy_time_interval_this;
453 } else {
454 AMREX_ALWAYS_ASSERT(bdy_time_interval == bdy_time_interval_this);
455 }
456 }
457
458 amrex::Print() << "Read in boundary data with width " << width;
459}
460
461/**
462 * \brief Helper function to initialize state and velocity data in a Fab from a REMORAdataset.
463 *
464 * @param lev Integer specifying current level
465 * @param state_fab FArrayBox object holding the state data we initialize
466 * @param temp_fab FArrayBox object holding the temperature data we initialize
467 * @param salt_fab FArrayBox object holding the salt data we initialize
468 * @param x_vel_fab FArrayBox object holding the x-velocity data we initialize
469 * @param y_vel_fab FArrayBox object holding the y-velocity data we initialize
470 * @param ubar_fab FArrayBox object holding the ubar data we initialize
471 * @param vbar_fab FArrayBox object holding the vbar data we initialize
472 * @param zeta_fab FArrayBox object holding the zeta data we initialize
473 * @param NC_temp_fab Vector of FArrayBox objects with the REMORA dataset specifying temperature
474 * @param NC_salt_fab Vector of FArrayBox objects with the REMORA dataset specifying salinity
475 * @param NC_xvel_fab Vector of FArrayBox objects with the REMORA dataset specifying x-velocity
476 * @param NC_yvel_fab Vector of FArrayBox objects with the REMORA dataset specifying y-velocity
477 * @param NC_ubar_fab Vector of FArrayBox objects with the REMORA dataset specifying ubar
478 * @param NC_vbar_fab Vector of FArrayBox objects with the REMORA dataset specifying vbar
479 * @param NC_zeta_fab Vector of FArrayBox objects with the REMORA dataset specifying zeta
480 */
481void
483 FArrayBox& temp_fab, FArrayBox& salt_fab,
484 FArrayBox& x_vel_fab, FArrayBox& y_vel_fab,
485 FArrayBox& ubar_fab, FArrayBox& vbar_fab,
486 const Vector<FArrayBox>& NC_temp_fab,
487 const Vector<FArrayBox>& NC_salt_fab,
488 const Vector<FArrayBox>& NC_xvel_fab,
489 const Vector<FArrayBox>& NC_yvel_fab,
490 const Vector<FArrayBox>& NC_ubar_fab,
491 const Vector<FArrayBox>& NC_vbar_fab)
492{
493 int nboxes = NC_xvel_fab.size();
494 for (int idx = 0; idx < nboxes; idx++)
495 {
496 //
497 // FArrayBox to FArrayBox copy does "copy on intersection"
498 // This only works here because we have broadcast the FArrayBox of data from the netcdf file to all ranks
499 //
500 temp_fab.template copy<RunOn::Device>(NC_temp_fab[idx]);
501 salt_fab.template copy<RunOn::Device>(NC_salt_fab[idx]);
502 x_vel_fab.template copy<RunOn::Device>(NC_xvel_fab[idx]);
503 y_vel_fab.template copy<RunOn::Device>(NC_yvel_fab[idx]);
504 ubar_fab.template copy<RunOn::Device>(NC_ubar_fab[idx],0,0,1);
505 vbar_fab.template copy<RunOn::Device>(NC_vbar_fab[idx],0,0,1);
506 } // idx
507}
508
509/**
510 * @param lev Integer specifying the current level
511 */
512void
514{
515 // *** FArrayBox's at this level for holding the INITIAL data
516 Vector<FArrayBox> NC_M2NC_fab ; NC_M2NC_fab.resize(num_boxes_at_level[lev]);
517 Vector<FArrayBox> NC_M3NC_fab ; NC_M3NC_fab.resize(num_boxes_at_level[lev]);
518 Vector<FArrayBox> NC_TempNC_fab ; NC_TempNC_fab.resize(num_boxes_at_level[lev]);
519 Vector<FArrayBox> NC_SaltNC_fab ; NC_SaltNC_fab.resize(num_boxes_at_level[lev]);
520
521 for (int idx = 0; idx < num_boxes_at_level[lev]; idx++)
522 {
528 NC_M2NC_fab[idx],NC_M3NC_fab[idx],
529 NC_TempNC_fab[idx],NC_SaltNC_fab[idx]);
530
531#ifdef _OPENMP
532#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
533#endif
534 {
535 // Don't tile this since we are operating on full FABs in this routine
536 for ( MFIter mfi(*cons_new[lev], false); mfi.isValid(); ++mfi )
537 {
539 FArrayBox &ubarNC_fab = (*vec_nudg_coeff[BdyVars::ubar][lev])[mfi];
540 ubarNC_fab.template copy<RunOn::Device>(NC_M2NC_fab[idx]);
541 FArrayBox &vbarNC_fab = (*vec_nudg_coeff[BdyVars::vbar][lev])[mfi];
542 vbarNC_fab.template copy<RunOn::Device>(NC_M2NC_fab[idx]);
543 }
545 FArrayBox &uNC_fab = (*vec_nudg_coeff[BdyVars::u][lev])[mfi];
546 uNC_fab.template copy<RunOn::Device>(NC_M3NC_fab[idx]);
547 FArrayBox &vNC_fab = (*vec_nudg_coeff[BdyVars::v][lev])[mfi];
548 vNC_fab.template copy<RunOn::Device>(NC_M3NC_fab[idx]);
549 }
551 FArrayBox &TempNC_fab = (*vec_nudg_coeff[BdyVars::t][lev])[mfi];
552 TempNC_fab.template copy<RunOn::Device>(NC_TempNC_fab[idx]);
553 }
555 FArrayBox &SaltNC_fab = (*vec_nudg_coeff[BdyVars::s][lev])[mfi];
556 SaltNC_fab.template copy<RunOn::Device>(NC_SaltNC_fab[idx]);
557 }
558
559 } // mf
560 } // omp
561 } // idx
562
564 vec_nudg_coeff[BdyVars::ubar][lev]->FillBoundary(geom[lev].periodicity());
565 vec_nudg_coeff[BdyVars::vbar][lev]->FillBoundary(geom[lev].periodicity());
568 }
570 vec_nudg_coeff[BdyVars::u][lev]->FillBoundary(geom[lev].periodicity());
571 vec_nudg_coeff[BdyVars::v][lev]->FillBoundary(geom[lev].periodicity());
574 }
576 vec_nudg_coeff[BdyVars::t][lev]->FillBoundary(geom[lev].periodicity());
578 }
580 vec_nudg_coeff[BdyVars::s][lev]->FillBoundary(geom[lev].periodicity());
582 }
583}
584
585/**
586 * @param[in ] lev level to read in river data
587 */
588void
590{
591 amrex::Vector<int> river_pos_x;
592 amrex::Vector<int> river_pos_y;
593 amrex::Vector<int> river_direction_tmp;
594
595 std::string river_x_name = "river_Xposition";
596 std::string river_y_name = "river_Eposition";
597 std::string river_dir_name = "river_direction";
598
599 read_vec_from_netcdf(lev, nc_riv_file, river_x_name, river_pos_x);
600 read_vec_from_netcdf(lev, nc_riv_file, river_y_name, river_pos_y);
601 read_vec_from_netcdf(lev, nc_riv_file, river_dir_name, river_direction_tmp);
602
603 int nriv = river_pos_x.size();
604 amrex::Gpu::DeviceVector<int> xpos_d(nriv);
605 amrex::Gpu::DeviceVector<int> ypos_d(nriv);
606 river_direction.resize(nriv);
607#ifdef AMREX_USE_GPU
608 Gpu::htod_memcpy(xpos_d.data(), river_pos_x.data(), sizeof(int)*nriv);
609 Gpu::htod_memcpy(ypos_d.data(), river_pos_y.data(), sizeof(int)*nriv);
610 Gpu::htod_memcpy(river_direction.data(), river_direction_tmp.data(), sizeof(int)*nriv);
611#else
612 std::memcpy(xpos_d.data(), river_pos_x.data(), sizeof(int)*nriv);
613 std::memcpy(ypos_d.data(), river_pos_y.data(), sizeof(int)*nriv);
614 std::memcpy(river_direction.data(), river_direction_tmp.data(), sizeof(int)*nriv);
615#endif
616 const int* xpos_ptr = xpos_d.data();
617 const int* ypos_ptr = ypos_d.data();
618
619 for (amrex::MFIter mfi(*(vec_river_position[lev]).get(),true); mfi.isValid(); ++mfi) {
620 amrex::Box bx = mfi.growntilebox(amrex::IntVect(NGROW,NGROW,0));
621 auto river_pos = vec_river_position[lev]->array(mfi);
622 ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int ) {
623 for (int iriv=0; iriv < nriv; iriv++) {
624 int xriv = xpos_ptr[iriv]-1;
625 int yriv = ypos_ptr[iriv]-1;
626 if (i==xriv && j==yriv) {
627 river_pos(i,j,0) = iriv;
628 }
629 }
630 });
631 }
632}
633
634/**
635 * @param[inout] mf multifab of data to convert
636 */
637void
639 Real inv_days_to_inv_s = 1.0_rt / (3600._rt * 24._rt);
640
641 for ( MFIter mfi(*mf, TilingIfNotGPU()); mfi.isValid(); ++mfi )
642 {
643 Array4<Real> const& arr = mf->array(mfi);
644 Box bx = mfi.growntilebox(IntVect(NGROW,NGROW,0));
645 ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
646 arr(i,j,k) *= inv_days_to_inv_s;
647 });
648 }
649
650}
651
652#endif // REMORA_USE_NETCDF
#define NGROW
#define Temp_comp
#define Salt_comp
void init_state_from_netcdf(int lev, FArrayBox &temp_fab, FArrayBox &salt_fab, FArrayBox &x_vel_fab, FArrayBox &y_vel_fab, FArrayBox &ubar_fab, FArrayBox &vbar_fab, const Vector< FArrayBox > &NC_temp_fab, const Vector< FArrayBox > &NC_salt_fab, const Vector< FArrayBox > &NC_xvel_fab, const Vector< FArrayBox > &NC_yvel_fab, const Vector< FArrayBox > &NC_ubar_fab, const Vector< FArrayBox > &NC_vbar_fab)
helper function to initialize state from netcdf
void read_coriolis_from_netcdf(int lev, const Box &domain, const std::string &fname, FArrayBox &NC_fcor_fab)
helper function to read coriolis factor from netcdf
void read_vec_from_netcdf(int lev, const std::string &fname, const std::string &field_name, amrex::Vector< int > &vec_dat)
helper function to read in vector of data from netcdf
void read_data_from_netcdf(int, const Box &domain, const std::string &fname, FArrayBox &NC_temp_fab, FArrayBox &NC_salt_fab, FArrayBox &NC_xvel_fab, FArrayBox &NC_yvel_fab, FArrayBox &NC_ubar_fab, FArrayBox &NC_vbar_fab)
helper function for reading in initial state data from netcdf
void read_zeta_from_netcdf(int lev, const Box &domain, const std::string &fname, FArrayBox &NC_zeta_fab)
helper function to read sea surface height from netcdf
Real read_bdry_from_netcdf(const Box &domain, const std::string &fname, Vector< Vector< FArrayBox > > &bdy_data_xlo, Vector< Vector< FArrayBox > > &bdy_data_xhi, Vector< Vector< FArrayBox > > &bdy_data_ylo, Vector< Vector< FArrayBox > > &bdy_data_yhi, int &width, amrex::Real &start_bdy_time, std::string bdry_time_varname, amrex::GpuArray< amrex::GpuArray< bool, AMREX_SPACEDIM *2 >, BdyVars::NumTypes+1 > &)
helper function for reading boundary data from netcdf
void read_masks_from_netcdf(int, const Box &domain, const std::string &fname, FArrayBox &NC_mskr_fab, FArrayBox &NC_msku_fab, FArrayBox &NC_mskv_fab)
helper function for reading in land-sea masks from netcdf
void read_bathymetry_from_netcdf(int lev, const Box &domain, const std::string &fname, FArrayBox &NC_h_fab, FArrayBox &NC_pm_fab, FArrayBox &NC_pn_fab, FArrayBox &NC_xr_fab, FArrayBox &NC_yr_fab, FArrayBox &NC_xu_fab, FArrayBox &NC_yu_fab, FArrayBox &NC_xv_fab, FArrayBox &NC_yv_fab, FArrayBox &NC_xp_fab, FArrayBox &NC_yp_fab)
helper function to read bathymetry from netcdf
void read_clim_nudg_coeff_from_netcdf(int lev, const Box &domain, const std::string &fname, bool do_m2_clim_nudg, bool do_m3_clim_nudg, bool do_temp_clim_nudg, bool do_salt_clim_nudg, FArrayBox &NC_M2NC_fab, FArrayBox &NC_M3NC_fab, FArrayBox &NC_TempNC_fab, FArrayBox &NC_SaltNC_fab)
helper function to read climatology nudging from netcdf
void init_bathymetry_from_netcdf(int lev)
Bathymetry data initialization from NetCDF file.
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_fcor
coriolis factor (2D)
Definition REMORA.H:363
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_hOfTheConfusingName
Bathymetry data (2D, positive valued, h in ROMS)
Definition REMORA.H:232
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_pm
horizontal scaling factor: 1 / dx (2D)
Definition REMORA.H:358
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_yv
y_grid on v-points (2D)
Definition REMORA.H:378
amrex::Vector< amrex::MultiFab * > cons_new
multilevel data container for current step's scalar data: temperature, salinity, passive scalar
Definition REMORA.H:220
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_mskr
land/sea mask at cell centers (2D)
Definition REMORA.H:349
std::string bdry_time_varname
Name of time field for boundary data.
Definition REMORA.H:1294
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_yp
y_grid on psi-points (2D)
Definition REMORA.H:383
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_xr
x_grid on rho points (2D)
Definition REMORA.H:366
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_xv
x_grid on v-points (2D)
Definition REMORA.H:376
amrex::Vector< amrex::Vector< amrex::Box > > boxes_at_level
the boxes specified at each level by tagging criteria
Definition REMORA.H:1107
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_msku
land/sea mask at x-faces (2D)
Definition REMORA.H:351
amrex::Vector< amrex::Vector< amrex::FArrayBox > > bdy_data_yhi
Vectors (over time) of Vector (over variables) of FArrayBoxs for holding Northern boundary data from ...
Definition REMORA.H:1012
void init_data_from_netcdf(int lev)
Problem initialization from NetCDF file.
void init_masks_from_netcdf(int lev)
Mask data initialization from NetCDF file.
amrex::Vector< amrex::MultiFab * > yvel_new
multilevel data container for current step's y velocities (v in ROMS)
Definition REMORA.H:224
amrex::Real bdy_time_interval
Interval between boundary data times.
Definition REMORA.H:1017
void init_bdry_from_netcdf()
Boundary data initialization from NetCDF file.
amrex::Real start_bdy_time
Start time in the time series of boundary data.
Definition REMORA.H:1015
amrex::Vector< int > num_boxes_at_level
how many boxes specified at each level by tagging criteria
Definition REMORA.H:1103
amrex::Vector< amrex::MultiFab * > xvel_new
multilevel data container for current step's x velocities (u in ROMS)
Definition REMORA.H:222
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_mskp
land/sea mask at cell corners (2D)
Definition REMORA.H:355
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_mskv
land/sea mask at y-faces (2D)
Definition REMORA.H:353
amrex::Vector< std::unique_ptr< REMORAPhysBCFunct > > physbcs
Vector (over level) of functors to apply physical boundary conditions.
Definition REMORA.H:1124
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.
amrex::Vector< amrex::Vector< amrex::FArrayBox > > bdy_data_ylo
Vectors (over time) of Vector (over variables) of FArrayBoxs for holding Southern boundary data from ...
Definition REMORA.H:1010
std::string nc_clim_coeff_file
NetCDF climatology coefficient file.
Definition REMORA.H:1291
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.
void init_riv_pos_from_netcdf(int lev)
static amrex::Vector< std::string > nc_bdry_file
NetCDF boundary data.
Definition REMORA.H:48
std::string nc_riv_file
Definition REMORA.H:1286
void update_mskp(int lev)
Set psi-point mask to be consistent with rho-point mask.
void init_zeta_from_netcdf(int lev)
Sea-surface height data initialization from NetCDF file.
void init_coriolis_from_netcdf(int lev)
Coriolis parameter data initialization from NetCDF file.
amrex::Vector< amrex::Real > t_new
new time at each level
Definition REMORA.H:1114
static SolverChoice solverChoice
Container for algorithmic choices.
Definition REMORA.H:1237
amrex::Vector< std::unique_ptr< amrex::iMultiFab > > vec_river_position
iMultiFab for river positions; contents are indices of rivers
Definition REMORA.H:1052
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_xp
x_grid on psi-points (2D)
Definition REMORA.H:381
amrex::Vector< amrex::Vector< amrex::FArrayBox > > bdy_data_xhi
Vectors (over time) of Vector (over variables) of FArrayBoxs for holding Eastern boundary data from f...
Definition REMORA.H:1008
static amrex::Vector< amrex::Vector< std::string > > nc_grid_file
NetCDF grid file.
Definition REMORA.H:50
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_zeta
free surface height (2D)
Definition REMORA.H:346
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
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_vbar
barotropic y velocity (2D)
Definition REMORA.H:344
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_ubar
barotropic x velocity (2D)
Definition REMORA.H:342
void init_clim_nudg_coeff_from_netcdf(int lev)
Climatology nudging coefficient initialization from NetCDF file.
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_yu
y_grid on u-points (2D)
Definition REMORA.H:373
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_xu
x_grid on u-points (2D)
Definition REMORA.H:371
amrex::GpuArray< amrex::GpuArray< bool, AMREX_SPACEDIM *2 >, BdyVars::NumTypes+1 > phys_bc_need_data
These are flags that indicate whether we need to read in boundary data from file.
Definition REMORA.H:1150
amrex::Vector< amrex::Vector< std::unique_ptr< amrex::MultiFab > > > vec_nudg_coeff
Climatology nudging coefficients.
Definition REMORA.H:417
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_pn
horizontal scaling factor: 1 / dy (2D)
Definition REMORA.H:360
amrex::Vector< amrex::Vector< amrex::FArrayBox > > bdy_data_xlo
Vectors (over time) of Vector (over variables) of FArrayBoxs for holding Western boundary data from f...
Definition REMORA.H:1006
static amrex::Vector< amrex::Vector< std::string > > nc_init_file
NetCDF initialization file.
Definition REMORA.H:49
amrex::Vector< amrex::Real > t_old
old time at each level
Definition REMORA.H:1116
void convert_inv_days_to_inv_s(amrex::MultiFab *)
Convert data in a multifab from inverse days to inverse seconds.
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_yr
y_grid on rho points (2D)
Definition REMORA.H:368