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 full-domain high resolution bathymetry from netcdf */
62void
63read_bathymetry_full_domain_from_netcdf (const Box& domain, const std::string& fname,
64 FArrayBox& NC_h_fab, IntVect ngrow);
65
66/** \brief helper function to read coriolis factor from netcdf */
67void
68read_coriolis_from_netcdf (int lev, const Box& domain, const std::string& fname, FArrayBox& NC_fcor_fab);
69
70/** \brief helper function to read sea surface height from netcdf */
71void
72read_zeta_from_netcdf (int lev, const Box& domain, const std::string& fname,
73 FArrayBox& NC_zeta_fab);
74
75/** \brief helper function to read climatology nudging from netcdf */
76void
77read_clim_nudg_coeff_from_netcdf (int lev, const Box& domain, const std::string& fname,
78 bool do_m2_clim_nudg,
79 bool do_m3_clim_nudg,
80 bool do_temp_clim_nudg,
81 bool do_salt_clim_nudg,
82 FArrayBox& NC_M2NC_fab,
83 FArrayBox& NC_M3NC_fab,
84 FArrayBox& NC_TempNC_fab,
85 FArrayBox& NC_SaltNC_fab);
86
87/** \brief helper function to read in vector of data from netcdf */
88void read_vec_from_netcdf (int lev, const amrex::Vector<std::string>& fnames, const std::string& field_name, amrex::Vector<int>& vec_dat);
89
90/**
91 * @param lev Integer specifying the current level
92 */
93void
95{
96 // *** FArrayBox's at this level for holding the INITIAL data
97 Vector<FArrayBox> NC_temp_fab ; NC_temp_fab.resize(num_boxes_at_level[lev]);
98 Vector<FArrayBox> NC_salt_fab ; NC_salt_fab.resize(num_boxes_at_level[lev]);
99 Vector<FArrayBox> NC_xvel_fab ; NC_xvel_fab.resize(num_boxes_at_level[lev]);
100 Vector<FArrayBox> NC_yvel_fab ; NC_yvel_fab.resize(num_boxes_at_level[lev]);
101 Vector<FArrayBox> NC_ubar_fab ; NC_ubar_fab.resize(num_boxes_at_level[lev]);
102 Vector<FArrayBox> NC_vbar_fab ; NC_vbar_fab.resize(num_boxes_at_level[lev]);
103
104 for (int idx = 0; idx < num_boxes_at_level[lev]; idx++)
105 {
106 read_data_from_netcdf(lev, boxes_at_level[lev][idx], nc_init_file[lev][idx],
107 NC_temp_fab[idx], NC_salt_fab[idx],
108 NC_xvel_fab[idx], NC_yvel_fab[idx],
109 NC_ubar_fab[idx], NC_vbar_fab[idx]);
110 }
111
112
113 MultiFab mf_temp(*cons_new[lev], make_alias, Temp_comp, 1);
114 MultiFab mf_salt(*cons_new[lev], make_alias, Salt_comp, 1);
115
116#ifdef _OPENMP
117#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
118#endif
119 {
120 // Don't tile this since we are operating on full FABs in this routine
121 for ( MFIter mfi(*cons_new[lev], false); mfi.isValid(); ++mfi )
122 {
123 // Define fabs for holding the initial data
124 FArrayBox &temp_fab = mf_temp[mfi];
125 FArrayBox &salt_fab = mf_salt[mfi];
126 FArrayBox &xvel_fab = (*xvel_new[lev])[mfi];
127 FArrayBox &yvel_fab = (*yvel_new[lev])[mfi];
128 FArrayBox &ubar_fab = (*vec_ubar[lev])[mfi];
129 FArrayBox &vbar_fab = (*vec_vbar[lev])[mfi];
130
131 init_state_from_netcdf(lev, temp_fab, salt_fab,
132 xvel_fab, yvel_fab,
133 ubar_fab, vbar_fab,
134 NC_temp_fab, NC_salt_fab,
135 NC_xvel_fab, NC_yvel_fab,
136 NC_ubar_fab, NC_vbar_fab);
137 } // mf
138 } // omp
139
140 if (nscalar > 1) {
141 cons_new[lev]->setVal(0.0_rt, Tracer_comp + 1, nscalar - 1, cons_new[lev]->nGrowVect());
142 }
143}
144
145/**
146 * @param lev Integer specifying the current level
147 */
148void
150{
151 // *** FArrayBox's at this level for holding the INITIAL data
152 Vector<FArrayBox> NC_zeta_fab ; NC_zeta_fab.resize(num_boxes_at_level[lev]);
153
154 for (int idx = 0; idx < num_boxes_at_level[lev]; idx++)
155 {
156 read_zeta_from_netcdf(lev,boxes_at_level[lev][idx], nc_init_file[lev][idx],
157 NC_zeta_fab[idx]);
158
159#ifdef _OPENMP
160#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
161#endif
162 {
163 // Don't tile this since we are operating on full FABs in this routine
164 for ( MFIter mfi(*cons_new[lev], false); mfi.isValid(); ++mfi )
165 {
166 FArrayBox &zeta_fab = (*vec_zeta[lev])[mfi];
167
168 //
169 // FArrayBox to FArrayBox copy does "copy on intersection"
170 // This only works here because we have broadcast the FArrayBox of data from the netcdf file to all ranks
171 //
172
173 zeta_fab.template copy<RunOn::Device>(NC_zeta_fab[idx],0,0,1);
174 } // mf
175 } // omp
176 } // idx
177
178 vec_zeta[lev]->FillBoundary(geom[lev].periodicity());
179 (*physbcs[lev])(*vec_zeta[lev],*vec_mskr[lev].get(),0,1,vec_zeta[lev]->nGrowVect(),t_new[lev],zeta_bc(),0,*vec_zeta[lev],*vec_msku[lev],*vec_mskv[lev]);
180// (*physbcs[lev])(*vec_zeta[lev],*vec_mskr[lev].get(),1,1,vec_zeta[lev]->nGrowVect(),t_new[lev],BCVars::zeta_bc);
181// (*physbcs[lev])(*vec_zeta[lev],*vec_mskr[lev].get(),2,1,vec_zeta[lev]->nGrowVect(),t_new[lev],BCVars::zeta_bc);
182
184 Real told = t_new[lev];
185 fill_from_bdyfiles(lev, *vec_zeta[lev], *vec_mskr[lev], told, zeta_bc(),BdyVars::zeta,0,0,
186 *vec_zeta[lev]);
187 }
188 if (lev>0) {
189 FillPatch(lev, t_old[lev], *vec_zeta[lev], GetVecOfPtrs(vec_zeta), zeta_bc(), BdyVars::zeta,
190 0, false,false,0,0,0.0,*vec_zeta[lev]);
191 }
192// fill_from_bdyfiles(lev, *vec_zeta[lev], *vec_mskr[lev], told, BCVars::zeta_bc,BdyVars::zeta,1,1);
193// fill_from_bdyfiles(lev, *vec_zeta[lev], *vec_mskr[lev], told, BCVars::zeta_bc,BdyVars::zeta,2,2);
194}
195/**
196 * @param lev Integer specifying the current level
197 */
198void
200{
201 // *** FArrayBox's at this level for holding the INITIAL data
202 Vector<FArrayBox> NC_h_fab ; NC_h_fab.resize(num_boxes_at_level[lev]);
203 Vector<FArrayBox> NC_pm_fab ; NC_pm_fab.resize(num_boxes_at_level[lev]);
204 Vector<FArrayBox> NC_pn_fab ; NC_pn_fab.resize(num_boxes_at_level[lev]);
205
206 Vector<FArrayBox> NC_xr_fab ; NC_xr_fab.resize(num_boxes_at_level[lev]);
207 Vector<FArrayBox> NC_yr_fab ; NC_yr_fab.resize(num_boxes_at_level[lev]);
208 Vector<FArrayBox> NC_xu_fab ; NC_xu_fab.resize(num_boxes_at_level[lev]);
209 Vector<FArrayBox> NC_yu_fab ; NC_yu_fab.resize(num_boxes_at_level[lev]);
210 Vector<FArrayBox> NC_xv_fab ; NC_xv_fab.resize(num_boxes_at_level[lev]);
211 Vector<FArrayBox> NC_yv_fab ; NC_yv_fab.resize(num_boxes_at_level[lev]);
212 Vector<FArrayBox> NC_xp_fab ; NC_xp_fab.resize(num_boxes_at_level[lev]);
213 Vector<FArrayBox> NC_yp_fab ; NC_yp_fab.resize(num_boxes_at_level[lev]);
214
215 for (int idx = 0; idx < num_boxes_at_level[lev]; idx++)
216 {
218 NC_h_fab[idx],
219 NC_pm_fab[idx], NC_pn_fab[idx],
220 NC_xr_fab[idx], NC_yr_fab[idx],
221 NC_xu_fab[idx], NC_yu_fab[idx],
222 NC_xv_fab[idx], NC_yv_fab[idx],
223 NC_xp_fab[idx], NC_yp_fab[idx]);
224
225#ifdef _OPENMP
226#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
227#endif
228 {
229 // Don't tile this since we are operating on full FABs in this routine
230 for ( MFIter mfi(*cons_new[lev], false); mfi.isValid(); ++mfi )
231 {
232 FArrayBox &h_fab = (*vec_h[lev])[mfi];
233 FArrayBox &pm_fab = (*vec_pm[lev])[mfi];
234 FArrayBox &pn_fab = (*vec_pn[lev])[mfi];
235 FArrayBox &xr_fab = (*vec_xr[lev])[mfi];
236 FArrayBox &yr_fab = (*vec_yr[lev])[mfi];
237 FArrayBox &xu_fab = (*vec_xu[lev])[mfi];
238 FArrayBox &yu_fab = (*vec_yu[lev])[mfi];
239 FArrayBox &xv_fab = (*vec_xv[lev])[mfi];
240 FArrayBox &yv_fab = (*vec_yv[lev])[mfi];
241 FArrayBox &xp_fab = (*vec_xp[lev])[mfi];
242 FArrayBox &yp_fab = (*vec_yp[lev])[mfi];
243
244 //
245 // FArrayBox to FArrayBox copy does "copy on intersection"
246 // This only works here because we have broadcast the FArrayBox of data from the netcdf file to all ranks
247 //
248
249 // Copy into both components of h
250 h_fab.template copy<RunOn::Device>(NC_h_fab[idx],0,0,1);
251 h_fab.template copy<RunOn::Device>(NC_h_fab[idx],0,1,1);
252
253 pm_fab.template copy<RunOn::Device>(NC_pm_fab[idx]);
254 pn_fab.template copy<RunOn::Device>(NC_pn_fab[idx]);
255
256 xr_fab.template copy<RunOn::Device>(NC_xr_fab[idx]);
257 yr_fab.template copy<RunOn::Device>(NC_yr_fab[idx]);
258 xu_fab.template copy<RunOn::Device>(NC_xu_fab[idx]);
259 yu_fab.template copy<RunOn::Device>(NC_yu_fab[idx]);
260 xv_fab.template copy<RunOn::Device>(NC_xv_fab[idx]);
261 yv_fab.template copy<RunOn::Device>(NC_yv_fab[idx]);
262 xp_fab.template copy<RunOn::Device>(NC_xp_fab[idx]);
263 yp_fab.template copy<RunOn::Device>(NC_yp_fab[idx]);
264 } // mf
265 } // omp
266 } // idx
267
268 const double dummy_time = 0.0_rt;
269 // Unconditional foextrap will overwrite periodicity, but EnforcePeriodicity will
270 // be called on h afterwards
271 FillPatch(lev,dummy_time,*vec_h[lev],GetVecOfPtrs(vec_h),
273 BdyVars::null,0,false,true,1);
274 FillPatch(lev,dummy_time,*vec_h[lev],GetVecOfPtrs(vec_h),
276 BdyVars::null,1,false,true,1);
277
278 if (lev > 0) {
279 FillPatch(lev,dummy_time,*vec_pm[lev],GetVecOfPtrs(vec_pm),
281 BdyVars::null,0,false,true);
282 FillPatch(lev,dummy_time,*vec_pn[lev],GetVecOfPtrs(vec_pn),
284 BdyVars::null,0,false,true);
285 }
286
287
288 int ng = vec_pm[lev]->nGrow();
289
290 const auto& dom_lo = amrex::lbound(geom[lev].Domain());
291 const auto& dom_hi = amrex::ubound(geom[lev].Domain());
292
293 //
294 // We need values of pm and pn outside the domain so we fill
295 // them here with foextrap
296 //
297 // We first fill interior ghost cells because we will need to extrapolate
298 // from ghost cells inside the domain to ghost cells outside the domain
299 //
300 vec_pm[lev]->FillBoundary(geom[lev].periodicity());
301 vec_pn[lev]->FillBoundary(geom[lev].periodicity());
302
303 vec_xr[lev]->FillBoundary(geom[lev].periodicity());
304 vec_yr[lev]->FillBoundary(geom[lev].periodicity());
305 vec_xu[lev]->FillBoundary(geom[lev].periodicity());
306 vec_yu[lev]->FillBoundary(geom[lev].periodicity());
307 vec_xv[lev]->FillBoundary(geom[lev].periodicity());
308 vec_yv[lev]->FillBoundary(geom[lev].periodicity());
309 vec_xp[lev]->FillBoundary(geom[lev].periodicity());
310 vec_yp[lev]->FillBoundary(geom[lev].periodicity());
311
312 for ( MFIter mfi(*vec_pm[lev]); mfi.isValid(); ++mfi )
313 {
314 Box bx = mfi.tilebox();
315
316 auto pm_fab = vec_pm[lev]->array(mfi);
317 auto pn_fab = vec_pn[lev]->array(mfi);
318
319 Box gbx_lox = adjCellLo(bx,0,ng); gbx_lox.grow(1,ng); gbx_lox.setBig (0,dom_lo.x-2);
320 Box gbx_hix = adjCellHi(bx,0,ng); gbx_hix.grow(1,ng); gbx_hix.setSmall(0,dom_hi.x+2);
321 Box gbx_loy = adjCellLo(bx,1,ng); gbx_loy.grow(0,ng); gbx_loy.setBig (1,dom_lo.y-2);
322 Box gbx_hiy = adjCellHi(bx,1,ng); gbx_hiy.grow(0,ng); gbx_hiy.setSmall(1,dom_hi.y+2);
323
324 // if (gbx_lox.ok()) amrex::AllPrint() << "GBX_XLO " << gbx_lox << std::endl;
325 // if (gbx_hix.ok()) amrex::AllPrint() << "GBX_XHI " << gbx_hix << std::endl;
326 // if (gbx_loy.ok()) amrex::AllPrint() << "GBX_YLO " << gbx_loy << std::endl;
327 // if (gbx_hiy.ok()) amrex::AllPrint() << "GBX_YHI " << gbx_hiy << std::endl;
328
329 if (gbx_lox.ok()) {
330 ParallelFor(gbx_lox, [=] AMREX_GPU_DEVICE (int i, int j, int k)
331 {
332 pm_fab(i,j,k,0) = pm_fab(dom_lo.x-1,j,k,0);
333 pn_fab(i,j,k,0) = pn_fab(dom_lo.x-1,j,k,0);
334 });
335 }
336 if (gbx_hix.ok()) {
337 ParallelFor(gbx_hix, [=] AMREX_GPU_DEVICE (int i, int j, int k)
338 {
339 pm_fab(i,j,k,0) = pm_fab(dom_hi.x+1,j,k,0);
340 pn_fab(i,j,k,0) = pn_fab(dom_hi.x+1,j,k,0);
341 });
342 }
343 if (gbx_loy.ok()) {
344 ParallelFor(gbx_loy, [=] AMREX_GPU_DEVICE (int i, int j, int k)
345 {
346 pm_fab(i,j,k,0) = pm_fab(i,dom_lo.y-1,k,0);
347 pn_fab(i,j,k,0) = pn_fab(i,dom_lo.y-1,k,0);
348 });
349 }
350 if (gbx_hiy.ok()) {
351 ParallelFor(gbx_hiy, [=] AMREX_GPU_DEVICE (int i, int j, int k)
352 {
353 pm_fab(i,j,k,0) = pm_fab(i,dom_hi.y+1,k,0);
354 pn_fab(i,j,k,0) = pn_fab(i,dom_hi.y+1,k,0);
355 });
356 }
357 } // mfi
358
359 vec_h[lev]->FillBoundary(geom[lev].periodicity());
360}
361
362/**
363 * @param lev Integer specifying the current level
364 */
365void
367{
368 // *** FArrayBox's at this level for holding the INITIAL data
369 Vector<FArrayBox> NC_fcor_fab ; NC_fcor_fab.resize(num_boxes_at_level[lev]);
370
371 for (int idx = 0; idx < num_boxes_at_level[lev]; idx++)
372 {
374 NC_fcor_fab[idx]);
375
376#ifdef _OPENMP
377#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
378#endif
379 {
380 // Don't tile this since we are operating on full FABs in this routine
381 for ( MFIter mfi(*cons_new[lev], false); mfi.isValid(); ++mfi )
382 {
383 FArrayBox &fcor_fab = (*vec_fcor[lev])[mfi];
384
385 //
386 // FArrayBox to FArrayBox copy does "copy on intersection"
387 // This only works here because we have broadcast the FArrayBox of data from the netcdf file to all ranks
388 //
389
390 fcor_fab.template copy<RunOn::Device>(NC_fcor_fab[idx]);
391 } // mf
392 } // omp
393 } // idx
394 vec_fcor[lev]->FillBoundary(geom[lev].periodicity());
395}
396
397/**
398 * @param lev Integer specifying the current level
399 */
400void
402{
403 // *** FArrayBox's at this level for holding the INITIAL data
404 Vector<FArrayBox> NC_mskr_fab ; NC_mskr_fab.resize(num_boxes_at_level[lev]);
405 Vector<FArrayBox> NC_msku_fab ; NC_msku_fab.resize(num_boxes_at_level[lev]);
406 Vector<FArrayBox> NC_mskv_fab ; NC_mskv_fab.resize(num_boxes_at_level[lev]);
407
408 for (int idx = 0; idx < num_boxes_at_level[lev]; idx++)
409 {
410 read_masks_from_netcdf(lev,boxes_at_level[lev][idx], nc_grid_file[lev][idx],
411 NC_mskr_fab[idx],NC_msku_fab[idx],
412 NC_mskv_fab[idx]);
413
414#ifdef _OPENMP
415#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
416#endif
417 {
418 // Don't tile this since we are operating on full FABs in this routine
419 for ( MFIter mfi(*cons_new[lev], false); mfi.isValid(); ++mfi )
420 {
421 FArrayBox &mskr_fab = (*vec_mskr[lev])[mfi];
422 FArrayBox &msku_fab = (*vec_msku[lev])[mfi];
423 FArrayBox &mskv_fab = (*vec_mskv[lev])[mfi];
424
425 //
426 // FArrayBox to FArrayBox copy does "copy on intersection"
427 // This only works here because we have broadcast the FArrayBox of data from the netcdf file to all ranks
428 //
429
430 mskr_fab.template copy<RunOn::Device>(NC_mskr_fab[idx]);
431 msku_fab.template copy<RunOn::Device>(NC_msku_fab[idx]);
432 mskv_fab.template copy<RunOn::Device>(NC_mskv_fab[idx]);
433 } // mf
434 } // omp
435 } // idx
436
437 update_mskp(lev);
438 vec_mskr[lev]->FillBoundary(geom[lev].periodicity());
439 vec_msku[lev]->FillBoundary(geom[lev].periodicity());
440 vec_mskv[lev]->FillBoundary(geom[lev].periodicity());
441 vec_mskp[lev]->FillBoundary(geom[lev].periodicity());
442}
443
444/**
445 * @param lev Integer specifying the current level
446 */
447void
449{
450 if (nc_bdry_file.empty() || nc_bdry_file[0].empty()) {
451 amrex::Error("NetCDF boundary file name must be provided via input");
452 }
453
454 amrex::Vector<std::string> field_name = {"u", "v", "temp", "salt", "ubar", "vbar", "zeta"};
455 amrex::Vector<IntVect > index_types = {IntVect(1,0,0), IntVect(0,1,0),
456 IntVect(0,0,0), IntVect(0,0,0),
457 IntVect(1,0,0), IntVect(0,1,0),
458 IntVect(0,0,0)};
459 std::vector<bool > is_2d = {false, false, false, false, true, true, true};
460
461 amrex::Print() << "DOING INIT AT LEVEL " << lev << std::endl;
462 int rx = 1; int ry = 1;
463 if (lev > 0) {
464 for (int k = lev-1; k >= 0; k--) {
465 rx *= ref_ratio[k][0];
466 ry *= ref_ratio[k][1];
467 }
468 }
469 for (int ivar = 0; ivar < BdyVars::NumTypes; ivar++) {
470 boundary_series[lev].push_back(std::unique_ptr<NCTimeSeriesBoundary>(new NCTimeSeriesBoundary(lev, geom, nc_bdry_file, field_name[ivar],
472 index_types[ivar],
473 &phys_bc_need_data[ivar], is_2d[ivar], rx, ry)));
474 boundary_series[lev][ivar]->Initialize();
475 }
476}
477
478/**
479 * \brief Helper function to initialize state and velocity data in a Fab from a REMORAdataset.
480 *
481 * @param lev Integer specifying current level
482 * @param state_fab FArrayBox object holding the state data we initialize
483 * @param temp_fab FArrayBox object holding the temperature data we initialize
484 * @param salt_fab FArrayBox object holding the salt data we initialize
485 * @param x_vel_fab FArrayBox object holding the x-velocity data we initialize
486 * @param y_vel_fab FArrayBox object holding the y-velocity data we initialize
487 * @param ubar_fab FArrayBox object holding the ubar data we initialize
488 * @param vbar_fab FArrayBox object holding the vbar data we initialize
489 * @param zeta_fab FArrayBox object holding the zeta data we initialize
490 * @param NC_temp_fab Vector of FArrayBox objects with the REMORA dataset specifying temperature
491 * @param NC_salt_fab Vector of FArrayBox objects with the REMORA dataset specifying salinity
492 * @param NC_xvel_fab Vector of FArrayBox objects with the REMORA dataset specifying x-velocity
493 * @param NC_yvel_fab Vector of FArrayBox objects with the REMORA dataset specifying y-velocity
494 * @param NC_ubar_fab Vector of FArrayBox objects with the REMORA dataset specifying ubar
495 * @param NC_vbar_fab Vector of FArrayBox objects with the REMORA dataset specifying vbar
496 * @param NC_zeta_fab Vector of FArrayBox objects with the REMORA dataset specifying zeta
497 */
498void
500 FArrayBox& temp_fab, FArrayBox& salt_fab,
501 FArrayBox& x_vel_fab, FArrayBox& y_vel_fab,
502 FArrayBox& ubar_fab, FArrayBox& vbar_fab,
503 const Vector<FArrayBox>& NC_temp_fab,
504 const Vector<FArrayBox>& NC_salt_fab,
505 const Vector<FArrayBox>& NC_xvel_fab,
506 const Vector<FArrayBox>& NC_yvel_fab,
507 const Vector<FArrayBox>& NC_ubar_fab,
508 const Vector<FArrayBox>& NC_vbar_fab)
509{
510 int nboxes = NC_xvel_fab.size();
511 for (int idx = 0; idx < nboxes; idx++)
512 {
513 //
514 // FArrayBox to FArrayBox copy does "copy on intersection"
515 // This only works here because we have broadcast the FArrayBox of data from the netcdf file to all ranks
516 //
517 temp_fab.template copy<RunOn::Device>(NC_temp_fab[idx]);
518 salt_fab.template copy<RunOn::Device>(NC_salt_fab[idx]);
519 x_vel_fab.template copy<RunOn::Device>(NC_xvel_fab[idx]);
520 y_vel_fab.template copy<RunOn::Device>(NC_yvel_fab[idx]);
521 ubar_fab.template copy<RunOn::Device>(NC_ubar_fab[idx],0,0,1);
522 vbar_fab.template copy<RunOn::Device>(NC_vbar_fab[idx],0,0,1);
523 } // idx
524}
525
526/**
527 * @param lev Integer specifying the current level
528 */
529void
531{
532 // *** FArrayBox's at this level for holding the INITIAL data
533 Vector<FArrayBox> NC_M2NC_fab ; NC_M2NC_fab.resize(num_boxes_at_level[lev]);
534 Vector<FArrayBox> NC_M3NC_fab ; NC_M3NC_fab.resize(num_boxes_at_level[lev]);
535 Vector<FArrayBox> NC_TempNC_fab ; NC_TempNC_fab.resize(num_boxes_at_level[lev]);
536 Vector<FArrayBox> NC_SaltNC_fab ; NC_SaltNC_fab.resize(num_boxes_at_level[lev]);
537
538 for (int idx = 0; idx < num_boxes_at_level[lev]; idx++)
539 {
545 NC_M2NC_fab[idx],NC_M3NC_fab[idx],
546 NC_TempNC_fab[idx],NC_SaltNC_fab[idx]);
547
548#ifdef _OPENMP
549#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
550#endif
551 {
552 // Don't tile this since we are operating on full FABs in this routine
553 for ( MFIter mfi(*cons_new[lev], false); mfi.isValid(); ++mfi )
554 {
556 FArrayBox &ubarNC_fab = (*vec_nudg_coeff[BdyVars::ubar][lev])[mfi];
557 ubarNC_fab.template copy<RunOn::Device>(NC_M2NC_fab[idx]);
558 FArrayBox &vbarNC_fab = (*vec_nudg_coeff[BdyVars::vbar][lev])[mfi];
559 vbarNC_fab.template copy<RunOn::Device>(NC_M2NC_fab[idx]);
560 }
562 FArrayBox &uNC_fab = (*vec_nudg_coeff[BdyVars::u][lev])[mfi];
563 uNC_fab.template copy<RunOn::Device>(NC_M3NC_fab[idx]);
564 FArrayBox &vNC_fab = (*vec_nudg_coeff[BdyVars::v][lev])[mfi];
565 vNC_fab.template copy<RunOn::Device>(NC_M3NC_fab[idx]);
566 }
568 FArrayBox &TempNC_fab = (*vec_nudg_coeff[BdyVars::t][lev])[mfi];
569 TempNC_fab.template copy<RunOn::Device>(NC_TempNC_fab[idx]);
570 }
572 FArrayBox &SaltNC_fab = (*vec_nudg_coeff[BdyVars::s][lev])[mfi];
573 SaltNC_fab.template copy<RunOn::Device>(NC_SaltNC_fab[idx]);
574 }
575
576 } // mf
577 } // omp
578 } // idx
579
581 vec_nudg_coeff[BdyVars::ubar][lev]->FillBoundary(geom[lev].periodicity());
582 vec_nudg_coeff[BdyVars::vbar][lev]->FillBoundary(geom[lev].periodicity());
585 }
587 vec_nudg_coeff[BdyVars::u][lev]->FillBoundary(geom[lev].periodicity());
588 vec_nudg_coeff[BdyVars::v][lev]->FillBoundary(geom[lev].periodicity());
591 }
593 vec_nudg_coeff[BdyVars::t][lev]->FillBoundary(geom[lev].periodicity());
595 }
597 vec_nudg_coeff[BdyVars::s][lev]->FillBoundary(geom[lev].periodicity());
599 }
600}
601
602/**
603 * @param[in ] lev level to read in river data
604 */
605void
607{
608 amrex::Vector<int> river_pos_x;
609 amrex::Vector<int> river_pos_y;
610 amrex::Vector<int> river_direction_tmp;
611
612 std::string river_x_name = "river_Xposition";
613 std::string river_y_name = "river_Eposition";
614 std::string river_dir_name = "river_direction";
615
616 read_vec_from_netcdf(lev, nc_riv_file, river_x_name, river_pos_x);
617 read_vec_from_netcdf(lev, nc_riv_file, river_y_name, river_pos_y);
618 read_vec_from_netcdf(lev, nc_riv_file, river_dir_name, river_direction_tmp);
619
620 int nriv = river_pos_x.size();
621 amrex::Gpu::DeviceVector<int> xpos_d(nriv);
622 amrex::Gpu::DeviceVector<int> ypos_d(nriv);
623 river_direction.resize(nriv);
624#ifdef AMREX_USE_GPU
625 Gpu::htod_memcpy(xpos_d.data(), river_pos_x.data(), sizeof(int)*nriv);
626 Gpu::htod_memcpy(ypos_d.data(), river_pos_y.data(), sizeof(int)*nriv);
627 Gpu::htod_memcpy(river_direction.data(), river_direction_tmp.data(), sizeof(int)*nriv);
628#else
629 std::memcpy(xpos_d.data(), river_pos_x.data(), sizeof(int)*nriv);
630 std::memcpy(ypos_d.data(), river_pos_y.data(), sizeof(int)*nriv);
631 std::memcpy(river_direction.data(), river_direction_tmp.data(), sizeof(int)*nriv);
632#endif
633 const int* xpos_ptr = xpos_d.data();
634 const int* ypos_ptr = ypos_d.data();
635
636 for (amrex::MFIter mfi(*(vec_river_position[lev]).get(),true); mfi.isValid(); ++mfi) {
637 amrex::Box bx = mfi.growntilebox(amrex::IntVect(NGROW,NGROW,0));
638 auto river_pos = vec_river_position[lev]->array(mfi);
639 ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int ) {
640 for (int iriv=0; iriv < nriv; iriv++) {
641 int xriv = xpos_ptr[iriv]-1;
642 int yriv = ypos_ptr[iriv]-1;
643 if (i==xriv && j==yriv) {
644 river_pos(i,j,0) = iriv;
645 }
646 }
647 });
648 }
649}
650
651/**
652 * @param[inout] mf multifab of data to convert
653 */
654void
656 Real inv_days_to_inv_s = 1.0_rt / (3600._rt * 24._rt);
657
658 for ( MFIter mfi(*mf, TilingIfNotGPU()); mfi.isValid(); ++mfi )
659 {
660 Array4<Real> const& arr = mf->array(mfi);
661 Box bx = mfi.growntilebox(IntVect(NGROW,NGROW,0));
662 ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
663 arr(i,j,k) *= inv_days_to_inv_s;
664 });
665 }
666
667}
668
669void
671{
672 if (nc_grid_file_hires.empty()) {
673 Abort("Must specify high-resolution grid file when initializing from NetCDF and hires_grid_level > 0");
674 }
675 Vector<FArrayBox> NC_h_fab ; NC_h_fab.resize(1);
677
678 // Don't tile this since we are operating on full FABs in this routine
679 for ( MFIter mfi(*vec_h_full_domain[hires_grid_level], false); mfi.isValid(); ++mfi )
680 {
681 FArrayBox &h_fab = (*vec_h_full_domain[hires_grid_level])[mfi];
682 h_fab.template copy<RunOn::Device>(NC_h_fab[0]);
683 }
684
685 // Average down to fill levels below hires_grid_level. Use a special average_down so
686 // grow cells get populated by averaged down fine data
687 for (int lev=hires_grid_level-1; lev >= 0; lev--) {
689 }
690}
691
692#endif // REMORA_USE_NETCDF
#define NGROW
#define Temp_comp
#define Tracer_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_vec_from_netcdf(int lev, const amrex::Vector< std::string > &fnames, const std::string &field_name, amrex::Vector< int > &vec_dat)
helper function to read in vector of data 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_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_bathymetry_full_domain_from_netcdf(const Box &domain, const std::string &fname, FArrayBox &NC_h_fab, IntVect ngrow)
helper function to read full-domain high resolution bathymetry 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
A class to hold and interpolate time series data read from a NetCDF file.
std::string nc_grid_file_hires
Grid file for high resolution bathymetry.
Definition REMORA.H:1465
amrex::Vector< std::string > nc_riv_file
NetCDF river file(s)
Definition REMORA.H:1474
int foextrap_periodic_bc() const noexcept
Definition REMORA.H:1059
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:409
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_h
Bathymetry data (2D, positive valued, h in ROMS)
Definition REMORA.H:254
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_pm
horizontal scaling factor: 1 / dx (2D)
Definition REMORA.H:404
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_yv
y_grid on v-points (2D)
Definition REMORA.H:424
amrex::Vector< amrex::MultiFab * > cons_new
multilevel data container for current step's scalar data: temperature, salinity, passive tracer
Definition REMORA.H:242
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_mskr
land/sea mask at cell centers (2D)
Definition REMORA.H:393
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_yp
y_grid on psi-points (2D)
Definition REMORA.H:429
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_xr
x_grid on rho points (2D)
Definition REMORA.H:412
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_xv
x_grid on v-points (2D)
Definition REMORA.H:422
amrex::Vector< amrex::Vector< amrex::Box > > boxes_at_level
the boxes specified at each level by tagging criteria
Definition REMORA.H:1266
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_msku
land/sea mask at x-faces (2D)
Definition REMORA.H:395
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:246
int zeta_bc() const noexcept
Definition REMORA.H:1057
amrex::Vector< amrex::IntVect > cum_ref_ratios
Cumulative refinement ratio between level 0 and level i.
Definition REMORA.H:1469
amrex::Vector< int > num_boxes_at_level
how many boxes specified at each level by tagging criteria
Definition REMORA.H:1262
amrex::Vector< amrex::MultiFab * > xvel_new
multilevel data container for current step's x velocities (u in ROMS)
Definition REMORA.H:244
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_mskp
land/sea mask at cell corners (2D)
Definition REMORA.H:399
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_mskv
land/sea mask at y-faces (2D)
Definition REMORA.H:397
amrex::Vector< std::unique_ptr< REMORAPhysBCFunct > > physbcs
Vector (over level) of functors to apply physical boundary conditions.
Definition REMORA.H:1283
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.
std::string nc_clim_coeff_file
NetCDF climatology coefficient file.
Definition REMORA.H:1479
amrex::Vector< std::string > bdry_time_name_byvar
Name of time fields for boundary data.
Definition REMORA.H:1484
void init_riv_pos_from_netcdf(int lev)
static amrex::Vector< std::string > nc_bdry_file
NetCDF boundary data.
Definition REMORA.H:51
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.
void fill_from_bdyfiles(int lev, amrex::MultiFab &mf_to_fill, const amrex::MultiFab &mf_mask, const amrex::Real time, const int bccomp, const int bdy_var_type, const int icomp_to_fill, const int icomp_calc=0, const amrex::MultiFab &mf_calc=amrex::MultiFab(), const amrex::Real=amrex::Real(0.0))
Fill boundary data from netcdf file.
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_h_full_domain
Bathymetry data on the whole domain at each potential level.
Definition REMORA.H:257
int nscalar
Number of passive scalars carried in the state.
Definition REMORA.H:1344
void average_down_with_grow_cells(int lev, amrex::Vector< std::unique_ptr< amrex::MultiFab > > &mf)
Average down from level lev+1 to lev in mf, including grow cells.
Definition REMORA.cpp:1636
amrex::Vector< amrex::Real > t_new
new time at each level
Definition REMORA.H:1273
void init_bdry_from_netcdf(int lev)
Boundary data initialization from NetCDF file.
static SolverChoice solverChoice
Container for algorithmic choices.
Definition REMORA.H:1403
amrex::Vector< std::unique_ptr< amrex::iMultiFab > > vec_river_position
iMultiFab for river positions; contents are indices of rivers
Definition REMORA.H:1211
amrex::Vector< amrex::Vector< std::unique_ptr< NCTimeSeriesBoundary > > > boundary_series
Vector over BdyVars of boundary series data containers.
Definition REMORA.H:1207
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_xp
x_grid on psi-points (2D)
Definition REMORA.H:427
static amrex::Vector< amrex::Vector< std::string > > nc_grid_file
NetCDF grid file.
Definition REMORA.H:53
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_zeta
free surface height (2D)
Definition REMORA.H:390
amrex::Gpu::DeviceVector< int > river_direction
Vector over rivers of river direction: 0: u-face; 1: v-face; 2: w-face.
Definition REMORA.H:1213
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_vbar
barotropic y velocity (2D)
Definition REMORA.H:388
amrex::Box nc_hires_grid_box
Box for the full domain at nc_hires_grid_level.
Definition REMORA.H:1467
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_ubar
barotropic x velocity (2D)
Definition REMORA.H:386
void init_clim_nudg_coeff_from_netcdf(int lev)
Climatology nudging coefficient initialization from NetCDF file.
void init_bathymetry_full_domain_from_netcdf()
Full domain bathymetry data initialization from NetCDF file.
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_yu
y_grid on u-points (2D)
Definition REMORA.H:419
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_xu
x_grid on u-points (2D)
Definition REMORA.H:417
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:1309
amrex::Vector< amrex::Vector< std::unique_ptr< amrex::MultiFab > > > vec_nudg_coeff
Climatology nudging coefficients.
Definition REMORA.H:463
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_pn
horizontal scaling factor: 1 / dy (2D)
Definition REMORA.H:406
static amrex::Vector< amrex::Vector< std::string > > nc_init_file
NetCDF initialization file.
Definition REMORA.H:52
amrex::Vector< amrex::Real > t_old
old time at each level
Definition REMORA.H:1275
int hires_grid_level
Which level the high resolution bathymetry is at.
Definition REMORA.H:1463
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:414