Canoe
Comprehensive Atmosphere N' Ocean Engine
pnetcdf.cpp
Go to the documentation of this file.
1 // C/C++
2 #include <cstdio>
3 #include <cstdlib>
4 #include <iomanip>
5 #include <iostream>
6 #include <sstream>
7 #include <stdexcept>
8 #include <string>
9 
10 // canoe
11 #include <configure.hpp>
12 #include <impl.hpp>
13 
14 // athena
15 #include <athena/athena.hpp>
16 #include <athena/coordinates/coordinates.hpp>
17 #include <athena/globals.hpp>
18 #include <athena/mesh/mesh.hpp>
19 #include <athena/outputs/user_outputs.hpp>
20 
21 // math
22 #include <climath/core.h>
23 
24 // harp
25 #include <harp/radiation.hpp>
26 #include <harp/radiation_band.hpp>
27 
28 // outputs
29 #include "output_utils.hpp"
30 
31 // Only proceed if PNETCDF output enabled
32 #ifdef PNETCDFOUTPUT
33 
34 // External library headers
35 #include <mpi.h>
36 #include <pnetcdf.h>
37 
38 #define ERR \
39  { \
40  if (err != NC_NOERR) { \
41  printf("Error at %s:%d : %s\n", __FILE__, __LINE__, \
42  ncmpi_strerror(err)); \
43  } \
44  }
45 
46 //----------------------------------------------------------------------------------------
47 // PnetcdfOutput constructor
48 // destructor - not needed for this derived class
49 
50 PnetcdfOutput::PnetcdfOutput(OutputParameters oparams) : OutputType(oparams) {}
51 
52 //----------------------------------------------------------------------------------------
54 // \brief Cycles over all MeshBlocks and writes OutputData in PNETCDF format
55 // One timestep per file
56 
57 void PnetcdfOutput::WriteOutputFile(Mesh *pm, ParameterInput *pin, bool flag) {
58  // create filename: "file_basename"+"."+"file_id"+"."+XXXXX+".nc",
59  // where XXXXX = 5-digit file_number
60  auto pmeta = MetadataTable::GetInstance();
61 
62  std::string fname;
63  char number[6];
64  int err;
65  snprintf(number, sizeof(number), "%05d", output_params.file_number);
66 
67  fname.assign(output_params.file_basename);
68  fname.append(".");
69  fname.append(output_params.file_id);
70  fname.append(".");
71  fname.append(number);
72  fname.append(".nc");
73 
74  // 0. reference radius for spherical polar geometry
75  float radius;
76  if (strcmp(COORDINATE_SYSTEM, "spherical_polar") == 0) {
77  try {
78  radius = pin->GetReal("problem", "radius");
79  } catch (std::runtime_error &e) {
80  std::stringstream msg;
81  msg << "### FATAL ERROR in PnetcdfOutput::WriteOutputFile" << std::endl
82  << "Spherical polar coordinate system must define a reference radius "
83  "in section <problem>"
84  << std::endl
85  << "Use radius = XXX " << std::endl;
86  ATHENA_ERROR(msg);
87  }
88  }
89 
90  // 1. open file for output
91  int ifile;
92  err = ncmpi_create(MPI_COMM_WORLD, fname.c_str(), NC_CLOBBER, MPI_INFO_NULL,
93  &ifile);
94  ERR
95 
96  // 2. coordinate structure
97  size_t nx1 = pm->mesh_size.nx1;
98  size_t nx2 = pm->mesh_size.nx2;
99  size_t nx3 = pm->mesh_size.nx3;
100  size_t nx1f = nx1;
101  if (nx1 > 1) nx1f++;
102  size_t nx2f = nx2;
103  if (nx2 > 1) nx2f++;
104  size_t nx3f = nx3;
105  if (nx3 > 1) nx3f++;
107 
108 #if ENABLE_HARP
109  size_t nrays = pm->my_blocks(0)->pimpl->prad->radiance.GetDim3();
110 #endif
111 
112  // 2. define coordinate
113  int idt, idx1, idx2, idx3, idx1f, idx2f, idx3f, iray;
114  int ndims = 0;
115 
116  ncmpi_def_dim(ifile, "time", NC_UNLIMITED, &idt);
117 
118  ncmpi_def_dim(ifile, "x1", nx1, &idx1);
119  ndims++;
120  if (nx1 > 1) {
121  ncmpi_def_dim(ifile, "x1f", nx1f, &idx1f);
122  ndims++;
123  }
124 
125  ncmpi_def_dim(ifile, "x2", nx2, &idx2);
126  ndims++;
127  if (nx2 > 1) {
128  ncmpi_def_dim(ifile, "x2f", nx2f, &idx2f);
129  ndims++;
130  }
131 
132  ncmpi_def_dim(ifile, "x3", nx3, &idx3);
133  ndims++;
134  if (nx3 > 1) {
135  ncmpi_def_dim(ifile, "x3f", nx3f, &idx3f);
136  ndims++;
137  }
138 
139 #if ENABLE_HARP
140  if (nrays > 0) {
141  ncmpi_def_dim(ifile, "ray", nrays, &iray);
142  ndims += 2;
143  }
144 #endif
145 
146  // 3. define variables
147  int ivt, ivx1, ivx2, ivx3, ivx1f, ivx2f, ivx3f, imu, iphi;
148 
149  ncmpi_def_var(ifile, "time", NC_FLOAT, 1, &idt, &ivt);
150  ncmpi_put_att_text(ifile, ivt, "axis", 1, "T");
151  ncmpi_put_att_text(ifile, ivt, "long_name", 4, "time");
152  if (strcmp(COORDINATE_SYSTEM, "spherical_polar") == 0) {
153  ncmpi_put_att_text(ifile, ivt, "units", 25, "seconds since 1-1-1 0:0:0");
154  } else {
155  ncmpi_put_att_text(ifile, ivt, "units", 7, "seconds");
156  }
157 
158  ncmpi_def_var(ifile, "x1", NC_FLOAT, 1, &idx1, &ivx1);
159  ncmpi_put_att_text(ifile, ivx1, "axis", 1, "Z");
160  if (strcmp(COORDINATE_SYSTEM, "spherical_polar") == 0) {
161  ncmpi_put_att_text(ifile, ivx1, "units", 6, "meters");
162  ncmpi_put_att_text(ifile, ivx1, "long_name", 8, "altitude");
163  } else {
164  ncmpi_put_att_text(ifile, ivx1, "units", 6, "meters");
165  ncmpi_put_att_text(ifile, ivx1, "long_name", 27,
166  "Z-coordinate at cell center");
167  }
168  if (nx1 > 1) {
169  ncmpi_def_var(ifile, "x1f", NC_FLOAT, 1, &idx1f, &ivx1f);
170  ncmpi_put_att_text(ifile, ivx1f, "axis", 1, "Z");
171  if (strcmp(COORDINATE_SYSTEM, "spherical_polar") == 0) {
172  ncmpi_put_att_text(ifile, ivx1f, "units", 6, "meters");
173  ncmpi_put_att_text(ifile, ivx1f, "long_name", 8, "altitude");
174  } else {
175  ncmpi_put_att_text(ifile, ivx1f, "units", 6, "meters");
176  ncmpi_put_att_text(ifile, ivx1f, "long_name", 25,
177  "Z-coordinate at cell face");
178  }
179  }
180 
181  ncmpi_def_var(ifile, "x2", NC_FLOAT, 1, &idx2, &ivx2);
182  ncmpi_put_att_text(ifile, ivx2, "axis", 1, "Y");
183  if (strcmp(COORDINATE_SYSTEM, "spherical_polar") == 0) {
184  ncmpi_put_att_text(ifile, ivx2, "units", 13, "degrees_north");
185  ncmpi_put_att_text(ifile, ivx2, "long_name", 8, "latitude");
186  } else {
187  ncmpi_put_att_text(ifile, ivx2, "units", 6, "meters");
188  ncmpi_put_att_text(ifile, ivx2, "long_name", 27,
189  "Y-coordinate at cell center");
190  }
191  if (nx2 > 1) {
192  ncmpi_def_var(ifile, "x2f", NC_FLOAT, 1, &idx2f, &ivx2f);
193  ncmpi_put_att_text(ifile, ivx2f, "axis", 1, "Y");
194  if (strcmp(COORDINATE_SYSTEM, "spherical_polar") == 0) {
195  ncmpi_put_att_text(ifile, ivx2f, "units", 13, "degrees_north");
196  ncmpi_put_att_text(ifile, ivx2f, "long_name", 8, "latitude");
197  } else {
198  ncmpi_put_att_text(ifile, ivx2f, "units", 6, "meters");
199  ncmpi_put_att_text(ifile, ivx2f, "long_name", 25,
200  "Y-coordinate at cell face");
201  }
202  }
203 
204  ncmpi_def_var(ifile, "x3", NC_FLOAT, 1, &idx3, &ivx3);
205  ncmpi_put_att_text(ifile, ivx3, "axis", 1, "X");
206  if (strcmp(COORDINATE_SYSTEM, "spherical_polar") == 0) {
207  ncmpi_put_att_text(ifile, ivx3, "units", 12, "degrees_east");
208  ncmpi_put_att_text(ifile, ivx3, "long_name", 9, "longitude");
209  } else {
210  ncmpi_put_att_text(ifile, ivx3, "units", 6, "meters");
211  ncmpi_put_att_text(ifile, ivx3, "long_name", 27,
212  "X-coordinate at cell center");
213  }
214  if (nx3 > 1) {
215  ncmpi_def_var(ifile, "x3f", NC_FLOAT, 1, &idx3f, &ivx3f);
216  ncmpi_put_att_text(ifile, ivx3f, "axis", 1, "X");
217  if (strcmp(COORDINATE_SYSTEM, "spherical_polar") == 0) {
218  ncmpi_put_att_text(ifile, ivx3f, "units", 12, "degrees_east");
219  ncmpi_put_att_text(ifile, ivx3f, "long_name", 9, "longitude");
220  } else {
221  ncmpi_put_att_text(ifile, ivx3f, "units", 6, "meters");
222  ncmpi_put_att_text(ifile, ivx3f, "long_name", 25,
223  "X-coordinate at cell face");
224  }
225  }
226 
227 #if ENABLE_HARP
228  if (nrays > 0) {
229  ncmpi_def_var(ifile, "mu_out", NC_FLOAT, 1, &iray, &imu);
230  ncmpi_put_att_text(ifile, imu, "units", 1, "1");
231  ncmpi_put_att_text(ifile, imu, "long_name", 18, "cosine polar angle");
232  ncmpi_def_var(ifile, "phi_out", NC_FLOAT, 1, &iray, &iphi);
233  ncmpi_put_att_text(ifile, iphi, "units", 3, "rad");
234  ncmpi_put_att_text(ifile, iphi, "long_name", 15, "azimuthal angle");
235  }
236 #endif
237 
238  MeshBlock *pmb = pm->my_blocks(0);
239  LoadOutputData(pmb);
240  OutputData *pdata = pfirst_data_;
241 
242  // count total variables (vector variables are expanded into flat scalars)
243  int total_vars = 0;
244  while (pdata != nullptr) {
245  std::string grid = pmeta->GetGridType(pdata->name);
246  int nvar = get_num_variables(grid, pdata->data);
247  total_vars += nvar;
248  pdata = pdata->pnext;
249  }
250 
251  int iaxis[4] = {idt, idx1, idx2, idx3};
252  int iaxis1[4] = {idt, idx1f, idx2, idx3};
253  int iaxis2[4] = {idt, idx1, idx2f, idx3};
254  int iaxis3[4] = {idt, idx1, idx2, idx3f};
255  int iaxisr[4] = {idt, iray, idx2, idx3};
256  int *var_ids = new int[total_vars];
257  int *ivar = var_ids;
258 
259  pdata = pfirst_data_;
260  while (pdata != nullptr) {
261  std::string grid = pmeta->GetGridType(pdata->name);
262  std::string name, attr;
263  std::vector<std::string> varnames;
264  int nvar = get_num_variables(grid, pdata->data);
265 
266  for (int n = 0; n < nvar; ++n) {
267  size_t pos = pdata->name.find('?');
268  if (nvar == 1) { // SCALARS
269  if (pos < pdata->name.length()) { // find '?'
270  varnames.push_back(pdata->name.substr(0, pos) +
271  pdata->name.substr(pos + 1));
272  } else {
273  varnames.push_back(pdata->name);
274  }
275  } else { // VECTORS
276  char c[16];
277  snprintf(c, sizeof(c), "%d", n + 1);
278  if (pos < pdata->name.length()) { // find '?'
279  varnames.push_back(pdata->name.substr(0, pos) + c +
280  pdata->name.substr(pos + 1));
281  } else {
282  varnames.push_back(pdata->name + c);
283  }
284  }
285  }
286 
287  for (int n = 0; n < nvar; ++n) {
288  name = varnames[n];
289  if (grid == "RCC") // radiation rays
290  ncmpi_def_var(ifile, name.c_str(), NC_FLOAT, 4, iaxisr, ivar);
291  else if (grid == "CCF")
292  ncmpi_def_var(ifile, name.c_str(), NC_FLOAT, 4, iaxis1, ivar);
293  else if ((grid == "CFC") && (nx2 > 1))
294  ncmpi_def_var(ifile, name.c_str(), NC_FLOAT, 4, iaxis2, ivar);
295  else if ((grid == "FCC") && (nx3 > 1))
296  ncmpi_def_var(ifile, name.c_str(), NC_FLOAT, 4, iaxis3, ivar);
297  else if (grid == "--C")
298  ncmpi_def_var(ifile, name.c_str(), NC_FLOAT, 2, iaxis, ivar);
299  else if (grid == "--F")
300  ncmpi_def_var(ifile, name.c_str(), NC_FLOAT, 2, iaxis1, ivar);
301  else if (grid == "---")
302  ncmpi_def_var(ifile, name.c_str(), NC_FLOAT, 1, iaxis, ivar);
303  else
304  ncmpi_def_var(ifile, name.c_str(), NC_FLOAT, 4, iaxis, ivar);
305 
306  // set units
307  attr = pmeta->GetUnits(pdata->name);
308  if (attr != "") {
309  ncmpi_put_att_text(ifile, *ivar, "units", attr.length(), attr.c_str());
310  }
311 
312  // set long_name
313  attr = pmeta->GetLongName(pdata->name);
314  if (attr != "") {
315  ncmpi_put_att_text(ifile, *ivar, "long_name", attr.length(),
316  attr.c_str());
317  }
318 
319  ivar++;
320  }
321  pdata = pdata->pnext;
322  }
323 
324  err = ncmpi_put_att_text(ifile, NC_GLOBAL, "Conventions", 6, "COARDS");
325  ERR if (strcmp(COORDINATE_SYSTEM, "spherical_polar") == 0) {
326  err = ncmpi_put_att_float(ifile, NC_GLOBAL, "PlanetRadius", NC_FLOAT, 1,
327  &radius);
328  ERR
329  }
330 
331  err = ncmpi_enddef(ifile);
332  ERR
333 
334  // 4. allocate data buffer, MPI requests and status
335  int max_ncells = 0,
336  nmb = 0;
337  for (int b = 0; b < pm->nblocal; ++b) {
338  MeshBlock *pmb = pm->my_blocks(b);
339  int nf1 = pmb->ie - pmb->is + 2 * NGHOST;
340  int nf2 = pmb->je - pmb->js + 2 * NGHOST;
341  int nf3 = pmb->ke - pmb->ks + 2 * NGHOST;
342  max_ncells = std::max(max_ncells, nf1 * nf2 * nf3);
343  nmb++; // number of meshblocks this rank
344  }
345  int nbufs = nmb * (ndims + total_vars);
346  int *reqs = new int[nbufs];
347  int *stts = new int[nbufs];
348  float **buf = new float *[nbufs];
349  buf[0] = new float[nbufs * max_ncells];
350  for (int i = 0; i < nbufs; ++i) buf[i] = buf[0] + i * max_ncells;
351  int *ir = reqs;
352  float **ib = buf;
353 
354  // 5. first meshblock writes time
355  float time = (float)pm->time;
356  MPI_Offset stime = 0, etime = 1;
357  err = ncmpi_put_vara_float_all(ifile, ivt, &stime, &etime, &time);
358  ERR
359 
360  ClearOutputData(); // required when LoadOutputData() is used.
361 
362  // Loop over MeshBlocks
363  // do {
364  for (int b = 0; b < pm->nblocal; ++b) {
365  MeshBlock *pmb = pm->my_blocks(b);
366  // set start/end array indices depending on whether ghost zones are included
367  out_is = pmb->is;
368  out_ie = pmb->ie;
369  out_js = pmb->js;
370  out_je = pmb->je;
371  out_ks = pmb->ks;
372  out_ke = pmb->ke;
373 
374  // \todo include_ghost zones probably doesn't work with grids other than CCC
375  if (output_params.include_ghost_zones) {
376  out_is -= NGHOST;
377  out_ie += NGHOST;
378  if (out_js != out_je) {
379  out_js -= NGHOST;
380  out_je += NGHOST;
381  }
382  if (out_ks != out_ke) {
383  out_ks -= NGHOST;
384  out_ke += NGHOST;
385  }
386  }
387 
388  LoadOutputData(pmb);
389 
390  // 6. each meshblock writes coordinates and variables
391  MPI_Offset ncells1 = out_ie - out_is + 1;
392  MPI_Offset ncells2 = out_je - out_js + 1;
393  MPI_Offset ncells3 = out_ke - out_ks + 1;
394  MPI_Offset nfaces1 = ncells1;
395  if (ncells1 > 1) nfaces1++;
396  MPI_Offset nfaces2 = ncells2;
397  if (ncells2 > 1) nfaces2++;
398  MPI_Offset nfaces3 = ncells3;
399  if (ncells3 > 1) nfaces3++;
400 
401  MPI_Offset start[4] = {0, ncells1 * pmb->loc.lx1, ncells2 * pmb->loc.lx2,
402  ncells3 * pmb->loc.lx3};
403  MPI_Offset count[4] = {1, ncells1, ncells2, ncells3};
404  MPI_Offset count1[4] = {1, nfaces1, ncells2, ncells3};
405  MPI_Offset count2[4] = {1, ncells1, nfaces2, ncells3};
406  MPI_Offset count3[4] = {1, ncells1, ncells2, nfaces3};
407  MPI_Offset startr[4] = {0, 0, ncells2 * pmb->loc.lx2,
408  ncells3 * pmb->loc.lx3};
409 
410 #if ENABLE_HARP
411  MPI_Offset countr[4] = {1, (MPI_Offset)nrays, ncells2, ncells3};
412 #endif
413 
414  // MPI_Offset start_ax1[2] = {0, ncells1*pmb->loc.lx1};
415  // MPI_Offset count_ax1[2] = {1, ncells1};
416  // MPI_Offset count_ax1f[2] = {1, nfaces1};
417 
418  if (strcmp(COORDINATE_SYSTEM, "spherical_polar") == 0) {
419  for (int i = out_is; i <= out_ie; ++i)
420  (*ib)[i - out_is] = (float)(pmb->pcoord->x1v(i)) - radius;
421  } else {
422  for (int i = out_is; i <= out_ie; ++i)
423  (*ib)[i - out_is] = (float)(pmb->pcoord->x1v(i));
424  }
425  err = ncmpi_iput_vara_float(ifile, ivx1, start + 1, count + 1, *ib++, ir++);
426  ERR
427 
428  if (nx1 > 1) {
429  if (strcmp(COORDINATE_SYSTEM, "spherical_polar") == 0) {
430  for (int i = out_is; i <= out_ie + 1; ++i)
431  (*ib)[i - out_is] = (float)(pmb->pcoord->x1f(i)) - radius;
432  } else {
433  for (int i = out_is; i <= out_ie + 1; ++i)
434  (*ib)[i - out_is] = (float)(pmb->pcoord->x1f(i));
435  }
436  err = ncmpi_iput_vara_float(ifile, ivx1f, start + 1, count1 + 1, *ib++,
437  ir++);
438  ERR
439  }
440 
441  if (strcmp(COORDINATE_SYSTEM, "spherical_polar") == 0) {
442  for (int j = out_js; j <= out_je; ++j)
443  (*ib)[j - out_js] = 90. - (float)rad2deg(pmb->pcoord->x2v(j));
444  } else {
445  for (int j = out_js; j <= out_je; ++j)
446  (*ib)[j - out_js] = (float)(pmb->pcoord->x2v(j));
447  }
448  err = ncmpi_iput_vara_float(ifile, ivx2, start + 2, count + 2, *ib++, ir++);
449  ERR
450 
451  if (nx2 > 1) {
452  if (strcmp(COORDINATE_SYSTEM, "spherical_polar") == 0) {
453  for (int j = out_js; j <= out_je + 1; ++j)
454  (*ib)[j - out_js] = 90. - (float)rad2deg(pmb->pcoord->x2f(j));
455  } else {
456  for (int j = out_js; j <= out_je + 1; ++j)
457  (*ib)[j - out_js] = (float)(pmb->pcoord->x2f(j));
458  }
459  err = ncmpi_iput_vara_float(ifile, ivx2f, start + 2, count2 + 2, *ib++,
460  ir++);
461  ERR
462  }
463 
464  if (strcmp(COORDINATE_SYSTEM, "spherical_polar") == 0) {
465  for (int k = out_ks; k <= out_ke; ++k)
466  (*ib)[k - out_ks] = (float)rad2deg(pmb->pcoord->x3v(k));
467  } else {
468  for (int k = out_ks; k <= out_ke; ++k)
469  (*ib)[k - out_ks] = (float)(pmb->pcoord->x3v(k));
470  }
471  err = ncmpi_iput_vara_float(ifile, ivx3, start + 3, count + 3, *ib++, ir++);
472  ERR
473 
474  if (nx3 > 1) {
475  if (strcmp(COORDINATE_SYSTEM, "spherical_polar") == 0) {
476  for (int k = out_ks; k <= out_ke + 1; ++k)
477  (*ib)[k - out_ks] = (float)rad2deg(pmb->pcoord->x3f(k));
478  } else {
479  for (int k = out_ks; k <= out_ke + 1; ++k)
480  (*ib)[k - out_ks] = (float)(pmb->pcoord->x3f(k));
481  }
482  err = ncmpi_iput_vara_float(ifile, ivx3f, start + 3, count3 + 3, *ib++,
483  ir++);
484  ERR
485  }
486 
487 #if ENABLE_HARP
488  if (nrays > 0) {
489  int m = 0;
490  for (auto p : pmb->pimpl->prad->bands) {
491  for (int n = 0; n < p->getNumOutgoingRays(); ++n)
492  (*ib)[m++] = (float)(p->getCosinePolarAngle(n));
493  }
494  err = ncmpi_iput_var_float(ifile, imu, *ib++, ir++);
495  ERR
496 
497  m = 0;
498  for (auto p : pmb->pimpl->prad->bands) {
499  for (int n = 0; n < p->getNumOutgoingRays(); ++n)
500  (*ib)[m++] = (float)(p->getAzimuthalAngle(n));
501  }
502  err = ncmpi_iput_var_float(ifile, iphi, *ib++, ir++);
503  ERR
504  }
505 #endif
506 
507  ivar = var_ids;
508  pdata = pfirst_data_;
509  while (pdata != nullptr) {
510  std::string grid = pmeta->GetGridType(pdata->name);
511  int nvar = get_num_variables(grid, pdata->data);
512 
513  if (grid == "RCC") { // radiation rays
514 #if ENABLE_HARP
515  for (int n = 0; n < nvar; n++) {
516  float *it = *ib;
517  for (int m = 0; m < nrays; ++m)
518  for (int j = out_js; j <= out_je; ++j)
519  for (int k = out_ks; k <= out_ke; ++k)
520  *it++ = (float)pdata->data(n, m, k, j);
521  err = ncmpi_iput_vara_float(ifile, *ivar++, startr, countr, *ib++,
522  ir++);
523  ERR
524  }
525 #endif
526  } else if (grid == "CCF") {
527  for (int n = 0; n < nvar; n++) {
528  float *it = *ib;
529  for (int i = out_is; i <= out_ie + 1; ++i)
530  for (int j = out_js; j <= out_je; ++j)
531  for (int k = out_ks; k <= out_ke; ++k)
532  *it++ = (float)pdata->data(n, k, j, i);
533  err =
534  ncmpi_iput_vara_float(ifile, *ivar++, start, count1, *ib++, ir++);
535  ERR
536  }
537  } else if ((grid == "CFC") && (nx2 > 1)) {
538  for (int n = 0; n < nvar; n++) {
539  float *it = *ib;
540  for (int i = out_is; i <= out_ie; ++i)
541  for (int j = out_js; j <= out_je + 1; ++j)
542  for (int k = out_ks; k <= out_ke; ++k)
543  *it++ = (float)pdata->data(n, k, j, i);
544  err =
545  ncmpi_iput_vara_float(ifile, *ivar++, start, count2, *ib++, ir++);
546  ERR
547  }
548  } else if ((grid == "FCC") && (nx3 > 1)) {
549  for (int n = 0; n < nvar; n++) {
550  float *it = *ib;
551  for (int i = out_is; i <= out_ie; ++i)
552  for (int j = out_js; j <= out_je; ++j)
553  for (int k = out_ks; k <= out_ke + 1; ++k)
554  *it++ = (float)pdata->data(n, k, j, i);
555  err =
556  ncmpi_iput_vara_float(ifile, *ivar++, start, count3, *ib++, ir++);
557  ERR
558  }
559  } else if (grid == "--C") {
560  for (int n = 0; n < nvar; n++) {
561  float *it = *ib;
562  for (int i = out_is; i <= out_ie; ++i)
563  *it++ = (float)pdata->data(n, i);
564  err =
565  ncmpi_iput_vara_float(ifile, *ivar++, start, count, *ib++, ir++);
566  ERR
567  }
568  } else if (grid == "--F") {
569  for (int n = 0; n < nvar; n++) {
570  float *it = *ib;
571  for (int i = out_is; i <= out_ie + 1; ++i)
572  *it++ = (float)pdata->data(n, i);
573  err =
574  ncmpi_iput_vara_float(ifile, *ivar++, start, count1, *ib++, ir++);
575  ERR
576  }
577  } else if (grid == "---") {
578  for (int n = 0; n < nvar; n++) {
579  **ib = (float)pdata->data(n, 0);
580  err =
581  ncmpi_iput_vara_float(ifile, *ivar++, start, count, *ib++, ir++);
582  ERR
583  }
584  } else {
585  for (int n = 0; n < nvar; n++) {
586  float *it = *ib;
587  for (int i = out_is; i <= out_ie; ++i)
588  for (int j = out_js; j <= out_je; ++j)
589  for (int k = out_ks; k <= out_ke; ++k)
590  *it++ = (float)pdata->data(n, k, j, i);
591  err =
592  ncmpi_iput_vara_float(ifile, *ivar++, start, count, *ib++, ir++);
593  ERR
594  }
595  }
596 
597  pdata = pdata->pnext;
598  }
599 
600  ClearOutputData(); // required when LoadOutputData() is used.
601  }
602  //} while (LoadOutputData(pmb)); // end loop over MeshBLocks
603 
604  // 7. wait for all writings to complete
605  ncmpi_wait_all(ifile, nbufs, reqs, stts);
606  for (int i = 0; i < nbufs; ++i)
607  if (stts[i] != NC_NOERR) {
608  std::stringstream msg;
609  msg << "### FATAL ERROR in PnetcdfOutput::WriteOutputFile" << std::endl
610  << "Nonblocking write fails on request " << i << std::endl
611  << ncmpi_strerror(stts[i]);
612  ATHENA_ERROR(msg);
613  }
614 
615  // 8. close nc file
616  ncmpi_close(ifile);
617 
618  // 9. clear data buffers
619  delete[] buf[0];
620  delete[] buf;
621  delete[] var_ids;
622  delete[] reqs;
623  delete[] stts;
624 
625  // 10. increment counters
626  output_params.file_number++;
627  output_params.next_time += output_params.dt;
628  pin->SetInteger(output_params.block_name, "file_number",
629  output_params.file_number);
630  pin->SetReal(output_params.block_name, "next_time", output_params.next_time);
631 
632  return;
633 }
634 
635 #endif // PNETCDFOUTPUT
MetadataTable::GetInstance
static MetadataTable const * GetInstance()
Definition: output_utils.cpp:57
max
double max(double x1, double x2, double x3)
Definition: core.h:19
rad2deg
double rad2deg(double phi)
Definition: core.h:40
fit_ammonia.e
e
Definition: fit_ammonia.py:8
fit_ammonia.b
b
Definition: fit_ammonia.py:8
fit_ammonia.p
p
Definition: fit_ammonia.py:8
make_plots.j
j
Definition: make_plots.py:26
make_plots.time
time
Definition: make_plots.py:9
get_num_variables
int get_num_variables(std::string grid, AthenaArray< Real > const &data)
Definition: output_utils.cpp:12
output_utils.hpp
radiation_band.hpp