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,0,*vec_zeta[lev],*vec_msku[lev],*vec_mskv[lev]);
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 *vec_zeta[lev]);
178 }
179 if (lev>0) {
180 FillPatch(lev, t_old[lev], *vec_zeta[lev], GetVecOfPtrs(vec_zeta), BCVars::zeta_bc, BdyVars::zeta,
181 0, false,false,0,0,0.0,*vec_zeta[lev]);
182 }
183// fill_from_bdyfiles(*vec_zeta[lev], *vec_mskr[lev], told, BCVars::zeta_bc,BdyVars::zeta,1,1);
184// fill_from_bdyfiles(*vec_zeta[lev], *vec_mskr[lev], told, BCVars::zeta_bc,BdyVars::zeta,2,2);
185}
186/**
187 * @param lev Integer specifying the current level
188 */
189void
191{
192 // *** FArrayBox's at this level for holding the INITIAL data
193 Vector<FArrayBox> NC_h_fab ; NC_h_fab.resize(num_boxes_at_level[lev]);
194 Vector<FArrayBox> NC_pm_fab ; NC_pm_fab.resize(num_boxes_at_level[lev]);
195 Vector<FArrayBox> NC_pn_fab ; NC_pn_fab.resize(num_boxes_at_level[lev]);
196
197 Vector<FArrayBox> NC_xr_fab ; NC_xr_fab.resize(num_boxes_at_level[lev]);
198 Vector<FArrayBox> NC_yr_fab ; NC_yr_fab.resize(num_boxes_at_level[lev]);
199 Vector<FArrayBox> NC_xu_fab ; NC_xu_fab.resize(num_boxes_at_level[lev]);
200 Vector<FArrayBox> NC_yu_fab ; NC_yu_fab.resize(num_boxes_at_level[lev]);
201 Vector<FArrayBox> NC_xv_fab ; NC_xv_fab.resize(num_boxes_at_level[lev]);
202 Vector<FArrayBox> NC_yv_fab ; NC_yv_fab.resize(num_boxes_at_level[lev]);
203 Vector<FArrayBox> NC_xp_fab ; NC_xp_fab.resize(num_boxes_at_level[lev]);
204 Vector<FArrayBox> NC_yp_fab ; NC_yp_fab.resize(num_boxes_at_level[lev]);
205
206 for (int idx = 0; idx < num_boxes_at_level[lev]; idx++)
207 {
209 NC_h_fab[idx],
210 NC_pm_fab[idx], NC_pn_fab[idx],
211 NC_xr_fab[idx], NC_yr_fab[idx],
212 NC_xu_fab[idx], NC_yu_fab[idx],
213 NC_xv_fab[idx], NC_yv_fab[idx],
214 NC_xp_fab[idx], NC_yp_fab[idx]);
215
216#ifdef _OPENMP
217#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
218#endif
219 {
220 // Don't tile this since we are operating on full FABs in this routine
221 for ( MFIter mfi(*cons_new[lev], false); mfi.isValid(); ++mfi )
222 {
223 FArrayBox &h_fab = (*vec_h[lev])[mfi];
224 FArrayBox &pm_fab = (*vec_pm[lev])[mfi];
225 FArrayBox &pn_fab = (*vec_pn[lev])[mfi];
226 FArrayBox &xr_fab = (*vec_xr[lev])[mfi];
227 FArrayBox &yr_fab = (*vec_yr[lev])[mfi];
228 FArrayBox &xu_fab = (*vec_xu[lev])[mfi];
229 FArrayBox &yu_fab = (*vec_yu[lev])[mfi];
230 FArrayBox &xv_fab = (*vec_xv[lev])[mfi];
231 FArrayBox &yv_fab = (*vec_yv[lev])[mfi];
232 FArrayBox &xp_fab = (*vec_xp[lev])[mfi];
233 FArrayBox &yp_fab = (*vec_yp[lev])[mfi];
234
235 //
236 // FArrayBox to FArrayBox copy does "copy on intersection"
237 // This only works here because we have broadcast the FArrayBox of data from the netcdf file to all ranks
238 //
239
240 // Copy into both components of h
241 h_fab.template copy<RunOn::Device>(NC_h_fab[idx],0,0,1);
242 h_fab.template copy<RunOn::Device>(NC_h_fab[idx],0,1,1);
243
244 pm_fab.template copy<RunOn::Device>(NC_pm_fab[idx]);
245 pn_fab.template copy<RunOn::Device>(NC_pn_fab[idx]);
246
247 xr_fab.template copy<RunOn::Device>(NC_xr_fab[idx]);
248 yr_fab.template copy<RunOn::Device>(NC_yr_fab[idx]);
249 xu_fab.template copy<RunOn::Device>(NC_xu_fab[idx]);
250 yu_fab.template copy<RunOn::Device>(NC_yu_fab[idx]);
251 xv_fab.template copy<RunOn::Device>(NC_xv_fab[idx]);
252 yv_fab.template copy<RunOn::Device>(NC_yv_fab[idx]);
253 xp_fab.template copy<RunOn::Device>(NC_xp_fab[idx]);
254 yp_fab.template copy<RunOn::Device>(NC_yp_fab[idx]);
255 } // mf
256 } // omp
257 } // idx
258
259 const double dummy_time = 0.0_rt;
260 // Unconditional foextrap will overwrite periodicity, but EnforcePeriodicity will
261 // be called on h afterwards
262 FillPatch(lev,dummy_time,*vec_h[lev],GetVecOfPtrs(vec_h),
264 BdyVars::null,0,false,true,1);
265 FillPatch(lev,dummy_time,*vec_h[lev],GetVecOfPtrs(vec_h),
267 BdyVars::null,1,false,true,1);
268
269 if (lev > 0) {
270 FillPatch(lev,dummy_time,*vec_pm[lev],GetVecOfPtrs(vec_pm),
272 BdyVars::null,0,false,true);
273 FillPatch(lev,dummy_time,*vec_pn[lev],GetVecOfPtrs(vec_pn),
275 BdyVars::null,0,false,true);
276 }
277
278
279 int ng = vec_pm[lev]->nGrow();
280
281 const auto& dom_lo = amrex::lbound(geom[lev].Domain());
282 const auto& dom_hi = amrex::ubound(geom[lev].Domain());
283
284 //
285 // We need values of pm and pn outside the domain so we fill
286 // them here with foextrap
287 //
288 // We first fill interior ghost cells because we will need to extrapolate
289 // from ghost cells inside the domain to ghost cells outside the domain
290 //
291 vec_pm[lev]->FillBoundary(geom[lev].periodicity());
292 vec_pn[lev]->FillBoundary(geom[lev].periodicity());
293
294 vec_xr[lev]->FillBoundary(geom[lev].periodicity());
295 vec_yr[lev]->FillBoundary(geom[lev].periodicity());
296 vec_xu[lev]->FillBoundary(geom[lev].periodicity());
297 vec_yu[lev]->FillBoundary(geom[lev].periodicity());
298 vec_xv[lev]->FillBoundary(geom[lev].periodicity());
299 vec_yv[lev]->FillBoundary(geom[lev].periodicity());
300 vec_xp[lev]->FillBoundary(geom[lev].periodicity());
301 vec_yp[lev]->FillBoundary(geom[lev].periodicity());
302
303 for ( MFIter mfi(*vec_pm[lev]); mfi.isValid(); ++mfi )
304 {
305 Box bx = mfi.tilebox();
306
307 auto pm_fab = vec_pm[lev]->array(mfi);
308 auto pn_fab = vec_pn[lev]->array(mfi);
309
310 Box gbx_lox = adjCellLo(bx,0,ng); gbx_lox.grow(1,ng); gbx_lox.setBig (0,dom_lo.x-2);
311 Box gbx_hix = adjCellHi(bx,0,ng); gbx_hix.grow(1,ng); gbx_hix.setSmall(0,dom_hi.x+2);
312 Box gbx_loy = adjCellLo(bx,1,ng); gbx_loy.grow(0,ng); gbx_loy.setBig (1,dom_lo.y-2);
313 Box gbx_hiy = adjCellHi(bx,1,ng); gbx_hiy.grow(0,ng); gbx_hiy.setSmall(1,dom_hi.y+2);
314
315 // if (gbx_lox.ok()) amrex::AllPrint() << "GBX_XLO " << gbx_lox << std::endl;
316 // if (gbx_hix.ok()) amrex::AllPrint() << "GBX_XHI " << gbx_hix << std::endl;
317 // if (gbx_loy.ok()) amrex::AllPrint() << "GBX_YLO " << gbx_loy << std::endl;
318 // if (gbx_hiy.ok()) amrex::AllPrint() << "GBX_YHI " << gbx_hiy << std::endl;
319
320 if (gbx_lox.ok()) {
321 ParallelFor(gbx_lox, [=] AMREX_GPU_DEVICE (int i, int j, int k)
322 {
323 pm_fab(i,j,k,0) = pm_fab(dom_lo.x-1,j,k,0);
324 pn_fab(i,j,k,0) = pn_fab(dom_lo.x-1,j,k,0);
325 });
326 }
327 if (gbx_hix.ok()) {
328 ParallelFor(gbx_hix, [=] AMREX_GPU_DEVICE (int i, int j, int k)
329 {
330 pm_fab(i,j,k,0) = pm_fab(dom_hi.x+1,j,k,0);
331 pn_fab(i,j,k,0) = pn_fab(dom_hi.x+1,j,k,0);
332 });
333 }
334 if (gbx_loy.ok()) {
335 ParallelFor(gbx_loy, [=] AMREX_GPU_DEVICE (int i, int j, int k)
336 {
337 pm_fab(i,j,k,0) = pm_fab(i,dom_lo.y-1,k,0);
338 pn_fab(i,j,k,0) = pn_fab(i,dom_lo.y-1,k,0);
339 });
340 }
341 if (gbx_hiy.ok()) {
342 ParallelFor(gbx_hiy, [=] AMREX_GPU_DEVICE (int i, int j, int k)
343 {
344 pm_fab(i,j,k,0) = pm_fab(i,dom_hi.y+1,k,0);
345 pn_fab(i,j,k,0) = pn_fab(i,dom_hi.y+1,k,0);
346 });
347 }
348 } // mfi
349
350 vec_h[lev]->FillBoundary(geom[lev].periodicity());
351}
352
353/**
354 * @param lev Integer specifying the current level
355 */
356void
358{
359 // *** FArrayBox's at this level for holding the INITIAL data
360 Vector<FArrayBox> NC_fcor_fab ; NC_fcor_fab.resize(num_boxes_at_level[lev]);
361
362 for (int idx = 0; idx < num_boxes_at_level[lev]; idx++)
363 {
365 NC_fcor_fab[idx]);
366
367#ifdef _OPENMP
368#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
369#endif
370 {
371 // Don't tile this since we are operating on full FABs in this routine
372 for ( MFIter mfi(*cons_new[lev], false); mfi.isValid(); ++mfi )
373 {
374 FArrayBox &fcor_fab = (*vec_fcor[lev])[mfi];
375
376 //
377 // FArrayBox to FArrayBox copy does "copy on intersection"
378 // This only works here because we have broadcast the FArrayBox of data from the netcdf file to all ranks
379 //
380
381 fcor_fab.template copy<RunOn::Device>(NC_fcor_fab[idx]);
382 } // mf
383 } // omp
384 } // idx
385 vec_fcor[lev]->FillBoundary(geom[lev].periodicity());
386}
387
388/**
389 * @param lev Integer specifying the current level
390 */
391void
393{
394 // *** FArrayBox's at this level for holding the INITIAL data
395 Vector<FArrayBox> NC_mskr_fab ; NC_mskr_fab.resize(num_boxes_at_level[lev]);
396 Vector<FArrayBox> NC_msku_fab ; NC_msku_fab.resize(num_boxes_at_level[lev]);
397 Vector<FArrayBox> NC_mskv_fab ; NC_mskv_fab.resize(num_boxes_at_level[lev]);
398
399 for (int idx = 0; idx < num_boxes_at_level[lev]; idx++)
400 {
401 read_masks_from_netcdf(lev,boxes_at_level[lev][idx], nc_grid_file[lev][idx],
402 NC_mskr_fab[idx],NC_msku_fab[idx],
403 NC_mskv_fab[idx]);
404
405#ifdef _OPENMP
406#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
407#endif
408 {
409 // Don't tile this since we are operating on full FABs in this routine
410 for ( MFIter mfi(*cons_new[lev], false); mfi.isValid(); ++mfi )
411 {
412 FArrayBox &mskr_fab = (*vec_mskr[lev])[mfi];
413 FArrayBox &msku_fab = (*vec_msku[lev])[mfi];
414 FArrayBox &mskv_fab = (*vec_mskv[lev])[mfi];
415
416 //
417 // FArrayBox to FArrayBox copy does "copy on intersection"
418 // This only works here because we have broadcast the FArrayBox of data from the netcdf file to all ranks
419 //
420
421 mskr_fab.template copy<RunOn::Device>(NC_mskr_fab[idx]);
422 msku_fab.template copy<RunOn::Device>(NC_msku_fab[idx]);
423 mskv_fab.template copy<RunOn::Device>(NC_mskv_fab[idx]);
424 } // mf
425 } // omp
426 } // idx
427
428 update_mskp(lev);
429 vec_mskr[lev]->FillBoundary(geom[lev].periodicity());
430 vec_msku[lev]->FillBoundary(geom[lev].periodicity());
431 vec_mskv[lev]->FillBoundary(geom[lev].periodicity());
432 vec_mskp[lev]->FillBoundary(geom[lev].periodicity());
433}
434
435/**
436 * @param lev Integer specifying the current level
437 */
438void
440{
441 if (nc_bdry_file.empty()) {
442 amrex::Error("NetCDF boundary file name must be provided via input");
443 }
444
445 int width = 1;
446
447 amrex::Vector<std::string> field_name = {"u", "v", "temp", "salt", "ubar", "vbar", "zeta"};
448 amrex::Vector<IntVect > index_types = {IntVect(1,0,0), IntVect(0,1,0),
449 IntVect(0,0,0), IntVect(0,0,0),
450 IntVect(1,0,0), IntVect(0,1,0),
451 IntVect(0,0,0)};
452 std::vector<bool > is_2d = {false, false, false, false, true, true, true};
453 for (int ivar = 0; ivar < BdyVars::NumTypes; ivar++) {
454 boundary_series.push_back(new NCTimeSeriesBoundary(nc_bdry_file, field_name[ivar],
456 geom[0].Domain(), index_types[ivar],
457 &phys_bc_need_data[ivar], is_2d[ivar]));
458 boundary_series[ivar]->Initialize();
459 }
460 amrex::Print() << "Read in boundary data with width " << width << std::endl;
461}
462
463/**
464 * \brief Helper function to initialize state and velocity data in a Fab from a REMORAdataset.
465 *
466 * @param lev Integer specifying current level
467 * @param state_fab FArrayBox object holding the state data we initialize
468 * @param temp_fab FArrayBox object holding the temperature data we initialize
469 * @param salt_fab FArrayBox object holding the salt data we initialize
470 * @param x_vel_fab FArrayBox object holding the x-velocity data we initialize
471 * @param y_vel_fab FArrayBox object holding the y-velocity data we initialize
472 * @param ubar_fab FArrayBox object holding the ubar data we initialize
473 * @param vbar_fab FArrayBox object holding the vbar data we initialize
474 * @param zeta_fab FArrayBox object holding the zeta data we initialize
475 * @param NC_temp_fab Vector of FArrayBox objects with the REMORA dataset specifying temperature
476 * @param NC_salt_fab Vector of FArrayBox objects with the REMORA dataset specifying salinity
477 * @param NC_xvel_fab Vector of FArrayBox objects with the REMORA dataset specifying x-velocity
478 * @param NC_yvel_fab Vector of FArrayBox objects with the REMORA dataset specifying y-velocity
479 * @param NC_ubar_fab Vector of FArrayBox objects with the REMORA dataset specifying ubar
480 * @param NC_vbar_fab Vector of FArrayBox objects with the REMORA dataset specifying vbar
481 * @param NC_zeta_fab Vector of FArrayBox objects with the REMORA dataset specifying zeta
482 */
483void
485 FArrayBox& temp_fab, FArrayBox& salt_fab,
486 FArrayBox& x_vel_fab, FArrayBox& y_vel_fab,
487 FArrayBox& ubar_fab, FArrayBox& vbar_fab,
488 const Vector<FArrayBox>& NC_temp_fab,
489 const Vector<FArrayBox>& NC_salt_fab,
490 const Vector<FArrayBox>& NC_xvel_fab,
491 const Vector<FArrayBox>& NC_yvel_fab,
492 const Vector<FArrayBox>& NC_ubar_fab,
493 const Vector<FArrayBox>& NC_vbar_fab)
494{
495 int nboxes = NC_xvel_fab.size();
496 for (int idx = 0; idx < nboxes; idx++)
497 {
498 //
499 // FArrayBox to FArrayBox copy does "copy on intersection"
500 // This only works here because we have broadcast the FArrayBox of data from the netcdf file to all ranks
501 //
502 temp_fab.template copy<RunOn::Device>(NC_temp_fab[idx]);
503 salt_fab.template copy<RunOn::Device>(NC_salt_fab[idx]);
504 x_vel_fab.template copy<RunOn::Device>(NC_xvel_fab[idx]);
505 y_vel_fab.template copy<RunOn::Device>(NC_yvel_fab[idx]);
506 ubar_fab.template copy<RunOn::Device>(NC_ubar_fab[idx],0,0,1);
507 vbar_fab.template copy<RunOn::Device>(NC_vbar_fab[idx],0,0,1);
508 } // idx
509}
510
511/**
512 * @param lev Integer specifying the current level
513 */
514void
516{
517 // *** FArrayBox's at this level for holding the INITIAL data
518 Vector<FArrayBox> NC_M2NC_fab ; NC_M2NC_fab.resize(num_boxes_at_level[lev]);
519 Vector<FArrayBox> NC_M3NC_fab ; NC_M3NC_fab.resize(num_boxes_at_level[lev]);
520 Vector<FArrayBox> NC_TempNC_fab ; NC_TempNC_fab.resize(num_boxes_at_level[lev]);
521 Vector<FArrayBox> NC_SaltNC_fab ; NC_SaltNC_fab.resize(num_boxes_at_level[lev]);
522
523 for (int idx = 0; idx < num_boxes_at_level[lev]; idx++)
524 {
530 NC_M2NC_fab[idx],NC_M3NC_fab[idx],
531 NC_TempNC_fab[idx],NC_SaltNC_fab[idx]);
532
533#ifdef _OPENMP
534#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
535#endif
536 {
537 // Don't tile this since we are operating on full FABs in this routine
538 for ( MFIter mfi(*cons_new[lev], false); mfi.isValid(); ++mfi )
539 {
541 FArrayBox &ubarNC_fab = (*vec_nudg_coeff[BdyVars::ubar][lev])[mfi];
542 ubarNC_fab.template copy<RunOn::Device>(NC_M2NC_fab[idx]);
543 FArrayBox &vbarNC_fab = (*vec_nudg_coeff[BdyVars::vbar][lev])[mfi];
544 vbarNC_fab.template copy<RunOn::Device>(NC_M2NC_fab[idx]);
545 }
547 FArrayBox &uNC_fab = (*vec_nudg_coeff[BdyVars::u][lev])[mfi];
548 uNC_fab.template copy<RunOn::Device>(NC_M3NC_fab[idx]);
549 FArrayBox &vNC_fab = (*vec_nudg_coeff[BdyVars::v][lev])[mfi];
550 vNC_fab.template copy<RunOn::Device>(NC_M3NC_fab[idx]);
551 }
553 FArrayBox &TempNC_fab = (*vec_nudg_coeff[BdyVars::t][lev])[mfi];
554 TempNC_fab.template copy<RunOn::Device>(NC_TempNC_fab[idx]);
555 }
557 FArrayBox &SaltNC_fab = (*vec_nudg_coeff[BdyVars::s][lev])[mfi];
558 SaltNC_fab.template copy<RunOn::Device>(NC_SaltNC_fab[idx]);
559 }
560
561 } // mf
562 } // omp
563 } // idx
564
566 vec_nudg_coeff[BdyVars::ubar][lev]->FillBoundary(geom[lev].periodicity());
567 vec_nudg_coeff[BdyVars::vbar][lev]->FillBoundary(geom[lev].periodicity());
570 }
572 vec_nudg_coeff[BdyVars::u][lev]->FillBoundary(geom[lev].periodicity());
573 vec_nudg_coeff[BdyVars::v][lev]->FillBoundary(geom[lev].periodicity());
576 }
578 vec_nudg_coeff[BdyVars::t][lev]->FillBoundary(geom[lev].periodicity());
580 }
582 vec_nudg_coeff[BdyVars::s][lev]->FillBoundary(geom[lev].periodicity());
584 }
585}
586
587/**
588 * @param[in ] lev level to read in river data
589 */
590void
592{
593 amrex::Vector<int> river_pos_x;
594 amrex::Vector<int> river_pos_y;
595 amrex::Vector<int> river_direction_tmp;
596
597 std::string river_x_name = "river_Xposition";
598 std::string river_y_name = "river_Eposition";
599 std::string river_dir_name = "river_direction";
600
601 read_vec_from_netcdf(lev, nc_riv_file, river_x_name, river_pos_x);
602 read_vec_from_netcdf(lev, nc_riv_file, river_y_name, river_pos_y);
603 read_vec_from_netcdf(lev, nc_riv_file, river_dir_name, river_direction_tmp);
604
605 int nriv = river_pos_x.size();
606 amrex::Gpu::DeviceVector<int> xpos_d(nriv);
607 amrex::Gpu::DeviceVector<int> ypos_d(nriv);
608 river_direction.resize(nriv);
609#ifdef AMREX_USE_GPU
610 Gpu::htod_memcpy(xpos_d.data(), river_pos_x.data(), sizeof(int)*nriv);
611 Gpu::htod_memcpy(ypos_d.data(), river_pos_y.data(), sizeof(int)*nriv);
612 Gpu::htod_memcpy(river_direction.data(), river_direction_tmp.data(), sizeof(int)*nriv);
613#else
614 std::memcpy(xpos_d.data(), river_pos_x.data(), sizeof(int)*nriv);
615 std::memcpy(ypos_d.data(), river_pos_y.data(), sizeof(int)*nriv);
616 std::memcpy(river_direction.data(), river_direction_tmp.data(), sizeof(int)*nriv);
617#endif
618 const int* xpos_ptr = xpos_d.data();
619 const int* ypos_ptr = ypos_d.data();
620
621 for (amrex::MFIter mfi(*(vec_river_position[lev]).get(),true); mfi.isValid(); ++mfi) {
622 amrex::Box bx = mfi.growntilebox(amrex::IntVect(NGROW,NGROW,0));
623 auto river_pos = vec_river_position[lev]->array(mfi);
624 ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int ) {
625 for (int iriv=0; iriv < nriv; iriv++) {
626 int xriv = xpos_ptr[iriv]-1;
627 int yriv = ypos_ptr[iriv]-1;
628 if (i==xriv && j==yriv) {
629 river_pos(i,j,0) = iriv;
630 }
631 }
632 });
633 }
634}
635
636/**
637 * @param[inout] mf multifab of data to convert
638 */
639void
641 Real inv_days_to_inv_s = 1.0_rt / (3600._rt * 24._rt);
642
643 for ( MFIter mfi(*mf, TilingIfNotGPU()); mfi.isValid(); ++mfi )
644 {
645 Array4<Real> const& arr = mf->array(mfi);
646 Box bx = mfi.growntilebox(IntVect(NGROW,NGROW,0));
647 ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
648 arr(i,j,k) *= inv_days_to_inv_s;
649 });
650 }
651
652}
653
654#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
A class to hold and interpolate time series data read from a NetCDF file.
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:377
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_h
Bathymetry data (2D, positive valued, h in ROMS)
Definition REMORA.H:241
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_pm
horizontal scaling factor: 1 / dx (2D)
Definition REMORA.H:372
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_yv
y_grid on v-points (2D)
Definition REMORA.H:392
amrex::Vector< amrex::MultiFab * > cons_new
multilevel data container for current step's scalar data: temperature, salinity, passive scalar
Definition REMORA.H:229
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_mskr
land/sea mask at cell centers (2D)
Definition REMORA.H:363
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_yp
y_grid on psi-points (2D)
Definition REMORA.H:397
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_xr
x_grid on rho points (2D)
Definition REMORA.H:380
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_xv
x_grid on v-points (2D)
Definition REMORA.H:390
amrex::Vector< amrex::Vector< amrex::Box > > boxes_at_level
the boxes specified at each level by tagging criteria
Definition REMORA.H:1150
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_msku
land/sea mask at x-faces (2D)
Definition REMORA.H:365
amrex::Vector< NCTimeSeriesBoundary * > boundary_series
Vector over BdyVars of boundary series data containers.
Definition REMORA.H:1091
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:233
void init_bdry_from_netcdf()
Boundary data initialization from NetCDF file.
amrex::Vector< int > num_boxes_at_level
how many boxes specified at each level by tagging criteria
Definition REMORA.H:1146
amrex::Vector< amrex::MultiFab * > xvel_new
multilevel data container for current step's x velocities (u in ROMS)
Definition REMORA.H:231
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_mskp
land/sea mask at cell corners (2D)
Definition REMORA.H:369
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_mskv
land/sea mask at y-faces (2D)
Definition REMORA.H:367
amrex::Vector< std::unique_ptr< REMORAPhysBCFunct > > physbcs
Vector (over level) of functors to apply physical boundary conditions.
Definition REMORA.H:1167
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:1340
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.
amrex::Vector< std::string > bdry_time_name_byvar
Name of time fields for boundary data.
Definition REMORA.H:1345
void init_riv_pos_from_netcdf(int lev)
static amrex::Vector< std::string > nc_bdry_file
NetCDF boundary data.
Definition REMORA.H:51
std::string nc_riv_file
Definition REMORA.H:1335
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:1157
static SolverChoice solverChoice
Container for algorithmic choices.
Definition REMORA.H:1280
amrex::Vector< std::unique_ptr< amrex::iMultiFab > > vec_river_position
iMultiFab for river positions; contents are indices of rivers
Definition REMORA.H:1095
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_xp
x_grid on psi-points (2D)
Definition REMORA.H:395
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:360
amrex::Gpu::DeviceVector< int > river_direction
Vector over rivers of river direction: 0: u-face; 1: v-face; 2: w-face.
Definition REMORA.H:1097
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_vbar
barotropic y velocity (2D)
Definition REMORA.H:358
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_ubar
barotropic x velocity (2D)
Definition REMORA.H:356
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:387
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_xu
x_grid on u-points (2D)
Definition REMORA.H:385
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:1193
amrex::Vector< amrex::Vector< std::unique_ptr< amrex::MultiFab > > > vec_nudg_coeff
Climatology nudging coefficients.
Definition REMORA.H:431
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_pn
horizontal scaling factor: 1 / dy (2D)
Definition REMORA.H:374
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:1159
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:382