ticeModelVec.hh - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
(HTM) git clone git://src.adamsgaard.dk/pism
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
ticeModelVec.hh (22895B)
---
1 // Copyright (C) 2008--2019 Ed Bueler, Constantine Khroulev, and David Maxwell
2 //
3 // This file is part of PISM.
4 //
5 // PISM is free software; you can redistribute it and/or modify it under the
6 // terms of the GNU General Public License as published by the Free Software
7 // Foundation; either version 3 of the License, or (at your option) any later
8 // version.
9 //
10 // PISM is distributed in the hope that it will be useful, but WITHOUT ANY
11 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
12 // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
13 // details.
14 //
15 // You should have received a copy of the GNU General Public License
16 // along with PISM; if not, write to the Free Software
17 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
18
19 #ifndef __IceModelVec_hh
20 #define __IceModelVec_hh
21
22 #include <initializer_list>
23 #include <memory>
24 #include <cstdint> // uint64_t
25
26 #include <petscvec.h>
27 #include <gsl/gsl_interp.h> // gsl_interp_accel
28
29 #include "VariableMetadata.hh"
30 #include "pism/util/petscwrappers/Viewer.hh"
31 #include "Vector2.hh"
32 #include "StarStencil.hh"
33 #include "pism/util/petscwrappers/DM.hh"
34 #include "pism/util/petscwrappers/Vec.hh"
35 #include "pism/util/IceGrid.hh"
36 #include "pism/util/io/IO_Flags.hh"
37 #include "pism/pism_config.hh" // Pism_DEBUG
38 #include "pism/util/interpolation.hh" // InterpolationType
39
40 namespace pism {
41
42 class IceGrid;
43 class File;
44
45 //! What "kind" of a vector to create: with or without ghosts.
46 enum IceModelVecKind {WITHOUT_GHOSTS=0, WITH_GHOSTS=1};
47
48 struct Range {
49 double min, max;
50 };
51
52 // NB: Do not change the order of elements in this struct. IceModelVec2S::box() and
53 // IceModelVec2Int::int_box() depend on it.
54 template <typename T>
55 struct BoxStencil {
56 T ij, n, nw, w, sw, s, se, e, ne;
57 };
58
59 class PetscAccessible {
60 public:
61 virtual ~PetscAccessible() {}
62 virtual void begin_access() const = 0;
63 virtual void end_access() const = 0;
64 };
65
66 //! Makes sure that we call begin_access() and end_access() for all accessed IceModelVecs.
67 class AccessList {
68 public:
69 AccessList();
70 AccessList(std::initializer_list<const PetscAccessible *> vecs);
71 AccessList(const PetscAccessible &v);
72 ~AccessList();
73 void add(const PetscAccessible &v);
74 void add(const std::vector<const PetscAccessible*> vecs);
75 private:
76 std::vector<const PetscAccessible*> m_vecs;
77 };
78
79 /*!
80 * Interpolation helper. Does not check if points needed for interpolation are within the current
81 * processor's sub-domain.
82 */
83 template<class F, typename T>
84 T interpolate(const F &field, double x, double y) {
85 auto grid = field.grid();
86
87 int i_left = 0, i_right = 0, j_bottom = 0, j_top = 0;
88 grid->compute_point_neighbors(x, y, i_left, i_right, j_bottom, j_top);
89
90 auto w = grid->compute_interp_weights(x, y);
91
92 return (w[0] * field(i_left, j_bottom) +
93 w[1] * field(i_right, j_bottom) +
94 w[2] * field(i_right, j_top) +
95 w[3] * field(i_left, j_top));
96 }
97
98 //! \brief Abstract class for reading, writing, allocating, and accessing a
99 //! DA-based PETSc Vec (2D and 3D fields) from within IceModel.
100 /*!
101 @anchor icemodelvec_use
102
103 This class represents 2D and 3D fields in PISM. Its methods common to all
104 the derived classes can be split (roughly) into six kinds:
105
106 - memory allocation (create)
107 - point-wise access (begin_access(), end_access())
108 - arithmetic (range(), norm(), add(), shift(), scale(), set(), ...)
109 - setting or reading metadata (set_attrs(), metadata())
110 - file input/output (read, write, regrid)
111 - tracking whether a field was updated (get_state_counter(), inc_state_counter())
112
113 ## Memory allocation
114
115 Creating an IceModelVec... object does not allocate memory for storing it
116 (some IceModelVecs serve as "references" and don't have their own storage).
117 To complete IceModelVec... creation, use the "create()" method:
118
119 \code
120 IceModelVec2S var;
121 ierr = var.create(grid, "var_name", WITH_GHOSTS); CHKERRQ(ierr);
122 // var is ready to use
123 \endcode
124
125 ("WITH_GHOSTS" means "can be used in computations using map-plane neighbors
126 of grid points.)
127
128 It is usually a good idea to set variable metadata right after creating it.
129 The method set_attrs() is used throughout PISM to set commonly used
130 attributes.
131
132 ## Point-wise access
133
134 PETSc performs some pointer arithmetic magic to allow convenient indexing of
135 grid point values. Because of this one needs to surround the code using row,
136 column or level indexes with begin_access() and end_access() calls:
137
138 \code
139 double foo;
140 int i = 0, j = 0;
141 IceModelVec2S var;
142 // assume that var was allocated
143 ierr = var.begin_access(); CHKERRQ(ierr);
144 foo = var(i,j) * 2;
145 ierr = var.end_access(); CHKERRQ(ierr);
146 \endcode
147
148 Please see [this page](@ref computational_grid) for a discussion of the
149 organization of PISM's computational grid and examples of for-loops you will
150 probably put between begin_access() and end_access().
151
152 To ensure that ghost values are up to date add the following call
153 before the code using ghosts:
154
155 \code
156 ierr = var.update_ghosts(); CHKERRQ(ierr);
157 \endcode
158
159 ## Reading and writing variables
160
161 PISM can read variables either from files with data on a grid matching the
162 current grid (read()) or, using bilinear interpolation, from files
163 containing data on a different (but compatible) grid (regrid()).
164
165 To write a field to a "prepared" NetCDF file, use write(). (A file is prepared
166 if it contains all the necessary dimensions, coordinate variables and global
167 metadata.)
168
169 If you need to "prepare" a file, do:
170 \code
171 File file(grid.com, PISM_NETCDF3);
172 io::prepare_for_output(file, *grid.ctx());
173 \endcode
174
175 A note about NetCDF write performance: due to limitations of the NetCDF
176 (classic, version 3) format, it is significantly faster to
177 \code
178 for (all variables)
179 var.define(...);
180
181 for (all variables)
182 var.write(...);
183 \endcode
184
185 as opposed to
186
187 \code
188 for (all variables) {
189 var.define(...);
190 var.write(...);
191 }
192 \endcode
193
194 IceModelVec::define() is here so that we can use the first approach.
195
196 ## Tracking if a field changed
197
198 It is possible to track if a certain field changed with the help of
199 get_state_counter() and inc_state_counter() methods.
200
201 For example, PISM's SIA code re-computes the smoothed bed only if the bed
202 deformation code updated it:
203
204 \code
205 if (bed->get_state_counter() > bed_state_counter) {
206 ierr = bed_smoother->preprocess_bed(...); CHKERRQ(ierr);
207 bed_state_counter = bed->get_state_counter();
208 }
209 \endcode
210
211 The state counter is **not** updated automatically. For the code snippet above
212 to work, a bed deformation model has to call inc_state_counter() after an
213 update.
214 */
215 class IceModelVec : public PetscAccessible {
216 public:
217 IceModelVec();
218 virtual ~IceModelVec();
219
220 typedef std::shared_ptr<IceModelVec> Ptr;
221 typedef std::shared_ptr<const IceModelVec> ConstPtr;
222
223
224 virtual bool was_created() const;
225 IceGrid::ConstPtr grid() const;
226 unsigned int ndims() const;
227 std::vector<int> shape() const;
228 //! \brief Returns the number of degrees of freedom per grid point.
229 unsigned int ndof() const;
230 unsigned int stencil_width() const;
231 std::vector<double> levels() const;
232
233 virtual Range range() const;
234 double norm(int n) const;
235 std::vector<double> norm_all(int n) const;
236 virtual void add(double alpha, const IceModelVec &x);
237 virtual void squareroot();
238 virtual void shift(double alpha);
239 virtual void scale(double alpha);
240 // This is used in Python code (as a local-to-global replacement),
241 // but we should be able to get rid of it.
242 void copy_to_vec(petsc::DM::Ptr destination_da, Vec destination) const;
243 void copy_from_vec(Vec source);
244 virtual void copy_from(const IceModelVec &source);
245 Vec vec();
246 petsc::DM::Ptr dm() const;
247 virtual void set_name(const std::string &name);
248 const std::string& get_name() const;
249 void set_attrs(const std::string &pism_intent,
250 const std::string &long_name,
251 const std::string &units,
252 const std::string &glaciological_units,
253 const std::string &standard_name,
254 unsigned int component);
255 virtual void read_attributes(const std::string &filename, int component = 0);
256 virtual void define(const File &nc, IO_Type default_type = PISM_DOUBLE) const;
257
258 void read(const std::string &filename, unsigned int time);
259 void read(const File &nc, unsigned int time);
260
261 void write(const std::string &filename) const;
262 void write(const File &nc) const;
263
264 void regrid(const std::string &filename, RegriddingFlag flag,
265 double default_value = 0.0);
266 void regrid(const File &nc, RegriddingFlag flag,
267 double default_value = 0.0);
268
269 virtual void begin_access() const;
270 virtual void end_access() const;
271 virtual void update_ghosts();
272 virtual void update_ghosts(IceModelVec &destination) const;
273
274 petsc::Vec::Ptr allocate_proc0_copy() const;
275 void put_on_proc0(Vec onp0) const;
276 void get_from_proc0(Vec onp0);
277
278 void set(double c);
279
280 SpatialVariableMetadata& metadata(unsigned int N = 0);
281
282 const SpatialVariableMetadata& metadata(unsigned int N = 0) const;
283
284 int state_counter() const;
285 void inc_state_counter();
286 void set_time_independent(bool flag);
287
288 protected:
289
290 //! If true, report range when regridding.
291 bool m_report_range;
292
293 void global_to_local(petsc::DM::Ptr dm, Vec source, Vec destination) const;
294 virtual void read_impl(const File &nc, unsigned int time);
295 virtual void regrid_impl(const File &nc, RegriddingFlag flag,
296 double default_value = 0.0);
297 virtual void write_impl(const File &nc) const;
298
299 std::vector<double> m_zlevels;
300
301 //! Internal storage
302 petsc::Vec m_v;
303 std::string m_name;
304
305 //! stores metadata (NetCDF variable attributes)
306 std::vector<SpatialVariableMetadata> m_metadata;
307
308 IceGrid::ConstPtr m_grid;
309
310 unsigned int m_dof; //!< number of "degrees of freedom" per grid point
311 unsigned int m_da_stencil_width; //!< stencil width supported by the DA
312 bool m_has_ghosts; //!< m_has_ghosts == true means "has ghosts"
313 petsc::DM::Ptr m_da; //!< distributed mesh manager (DM)
314
315 bool m_begin_end_access_use_dof;
316
317 //! It is a map, because a temporary IceModelVec can be used to view
318 //! different quantities
319 mutable std::map<std::string,petsc::Viewer::Ptr> m_map_viewers;
320
321 mutable void *m_array; // will be cast to double** or double*** in derived classes
322
323 mutable int m_access_counter; // used in begin_access() and end_access()
324 int m_state_counter; //!< Internal IceModelVec "revision number"
325
326 InterpolationType m_interpolation_type;
327
328 virtual void checkCompatibility(const char *function, const IceModelVec &other) const;
329
330 //! \brief Check the array indices and warn if they are out of range.
331 void check_array_indices(int i, int j, unsigned int k) const;
332 void reset_attrs(unsigned int N);
333 NormType int_to_normtype(int input) const;
334
335 void get_dof(petsc::DM::Ptr da_result, Vec result, unsigned int n,
336 unsigned int count=1) const;
337 void set_dof(petsc::DM::Ptr da_source, Vec source, unsigned int n,
338 unsigned int count=1);
339 private:
340 size_t size() const;
341 // disable copy constructor and the assignment operator:
342 IceModelVec(const IceModelVec &other);
343 IceModelVec& operator=(const IceModelVec&);
344 public:
345 //! Dump an IceModelVec to a file. *This is for debugging only.*
346 //! Uses const char[] to make it easier to call it from gdb.
347 void dump(const char filename[]) const;
348
349 uint64_t fletcher64() const;
350 std::string checksum() const;
351 void print_checksum(const char *prefix = "") const;
352
353 typedef pism::AccessList AccessList;
354 protected:
355 void put_on_proc0(Vec parallel, Vec onp0) const;
356 void get_from_proc0(Vec onp0, Vec parallel);
357 };
358
359 bool set_contains(const std::set<std::string> &S, const IceModelVec &field);
360
361 class IceModelVec2S;
362
363 /** Class for a 2d DA-based Vec.
364
365 As for the difference between IceModelVec2 and IceModelVec2S, the
366 former can store fields with more than 1 "degree of freedom" per grid
367 point (such as 2D fields on the "staggered" grid, with the first
368 degree of freedom corresponding to the i-offset and second to
369 j-offset). */
370 class IceModelVec2 : public IceModelVec {
371 public:
372 IceModelVec2();
373 IceModelVec2(IceGrid::ConstPtr grid, const std::string &short_name,
374 IceModelVecKind ghostedp, unsigned int stencil_width, int dof);
375
376 typedef std::shared_ptr<IceModelVec2> Ptr;
377 typedef std::shared_ptr<const IceModelVec2> ConstPtr;
378
379 static Ptr To2D(IceModelVec::Ptr input);
380
381 virtual void view(int viewer_size) const;
382 virtual void view(petsc::Viewer::Ptr v1, petsc::Viewer::Ptr v2) const;
383 // component-wise access:
384 virtual void get_component(unsigned int n, IceModelVec2S &result) const;
385 virtual void set_component(unsigned int n, const IceModelVec2S &source);
386 inline double& operator() (int i, int j, int k);
387 inline const double& operator() (int i, int j, int k) const;
388 void create(IceGrid::ConstPtr grid, const std::string &short_name,
389 IceModelVecKind ghostedp, unsigned int stencil_width, int dof);
390 protected:
391 virtual void read_impl(const File &nc, const unsigned int time);
392 virtual void regrid_impl(const File &nc, RegriddingFlag flag,
393 double default_value = 0.0);
394 virtual void write_impl(const File &nc) const;
395 };
396
397 //! A "fat" storage vector for combining related fields (such as SSAFEM coefficients).
398 template<typename T>
399 class IceModelVec2Fat : public IceModelVec2 {
400 public:
401 IceModelVec2Fat() {
402 m_dof = sizeof(T) / sizeof(double);
403 m_begin_end_access_use_dof = false;
404 }
405
406 void create(IceGrid::ConstPtr grid, const std::string &short_name,
407 IceModelVecKind ghostedp, unsigned int stencil_width = 1) {
408
409 m_name = short_name;
410
411 IceModelVec2::create(grid, short_name, ghostedp, stencil_width, m_dof);
412 }
413
414 inline T& operator()(int i, int j) {
415 #if (Pism_DEBUG==1)
416 check_array_indices(i, j, 0);
417 #endif
418 return static_cast<T**>(m_array)[j][i];
419 }
420
421 inline const T& operator()(int i, int j) const {
422 #if (Pism_DEBUG==1)
423 check_array_indices(i, j, 0);
424 #endif
425 return static_cast<T**>(m_array)[j][i];
426 }
427
428 };
429
430
431 class IceModelVec2V;
432
433 /** A class for storing and accessing scalar 2D fields.
434 IceModelVec2S is just IceModelVec2 with "dof == 1" */
435 class IceModelVec2S : public IceModelVec2 {
436 friend class IceModelVec2V;
437 friend class IceModelVec2Stag;
438 public:
439 IceModelVec2S();
440 IceModelVec2S(IceGrid::ConstPtr grid, const std::string &name,
441 IceModelVecKind ghostedp, int width = 1);
442
443 typedef std::shared_ptr<IceModelVec2S> Ptr;
444 typedef std::shared_ptr<const IceModelVec2S> ConstPtr;
445
446 static Ptr To2DScalar(IceModelVec::Ptr input);
447
448 /*!
449 * Interpolation helper. See the pism::interpolate() for details.
450 */
451 double interpolate(double x, double y) const {
452 return pism::interpolate<IceModelVec2S, double>(*this, x, y);
453 }
454
455 // does not need a copy constructor, because it does not add any new data members
456 using IceModelVec2::create;
457 void create(IceGrid::ConstPtr grid, const std::string &name,
458 IceModelVecKind ghostedp, int width = 1);
459 virtual void copy_from(const IceModelVec &source);
460 double** get_array();
461 virtual void set_to_magnitude(const IceModelVec2S &v_x, const IceModelVec2S &v_y);
462 virtual void set_to_magnitude(const IceModelVec2V &input);
463 virtual void mask_by(const IceModelVec2S &M, double fill = 0.0);
464 virtual void add(double alpha, const IceModelVec &x);
465 virtual void add(double alpha, const IceModelVec &x, IceModelVec &result) const;
466 virtual double sum() const;
467 virtual double min() const;
468 virtual double max() const;
469 virtual double absmax() const;
470 virtual double diff_x(int i, int j) const;
471 virtual double diff_y(int i, int j) const;
472 virtual double diff_x_p(int i, int j) const;
473 virtual double diff_y_p(int i, int j) const;
474
475 //! Provides access (both read and write) to the internal double array.
476 /*!
477 Note that i corresponds to the x direction and j to the y.
478 */
479 inline double& operator() (int i, int j);
480 inline const double& operator()(int i, int j) const;
481 inline StarStencil<double> star(int i, int j) const;
482 inline BoxStencil<double> box(int i, int j) const;
483 };
484
485
486 //! \brief A simple class "hiding" the fact that the mask is stored as
487 //! floating-point scalars (instead of integers).
488 class IceModelVec2Int : public IceModelVec2S {
489 public:
490 IceModelVec2Int();
491 IceModelVec2Int(IceGrid::ConstPtr grid, const std::string &name,
492 IceModelVecKind ghostedp, int width = 1);
493
494 typedef std::shared_ptr<IceModelVec2Int> Ptr;
495 typedef std::shared_ptr<const IceModelVec2Int> ConstPtr;
496
497 inline int as_int(int i, int j) const;
498 inline StarStencil<int> int_star(int i, int j) const;
499 inline BoxStencil<int> int_box(int i, int j) const;
500 };
501
502 /** Class for storing and accessing 2D vector fields used in IceModel.
503 IceModelVec2V is IceModelVec2 with "dof == 2". (Plus some extra methods, of course.)
504 */
505 class IceModelVec2V : public IceModelVec2 {
506 public:
507 IceModelVec2V();
508 IceModelVec2V(IceGrid::ConstPtr grid, const std::string &short_name,
509 IceModelVecKind ghostedp, unsigned int stencil_width = 1);
510 ~IceModelVec2V();
511
512 typedef std::shared_ptr<IceModelVec2V> Ptr;
513 typedef std::shared_ptr<const IceModelVec2V> ConstPtr;
514
515 static Ptr ToVector(IceModelVec::Ptr input);
516
517 void create(IceGrid::ConstPtr grid, const std::string &short_name,
518 IceModelVecKind ghostedp, unsigned int stencil_width = 1);
519 virtual void copy_from(const IceModelVec &source);
520 virtual void add(double alpha, const IceModelVec &x);
521 virtual void add(double alpha, const IceModelVec &x, IceModelVec &result) const;
522
523 // I/O:
524 Vector2** get_array();
525 inline Vector2& operator()(int i, int j);
526 inline const Vector2& operator()(int i, int j) const;
527 inline StarStencil<Vector2> star(int i, int j) const;
528
529 /*!
530 * Interpolation helper. See the pism::interpolate() for details.
531 */
532 Vector2 interpolate(double x, double y) const {
533 return pism::interpolate<IceModelVec2V, Vector2>(*this, x, y);
534 }
535 };
536
537 //! \brief A class for storing and accessing internal staggered-grid 2D fields.
538 //! Uses dof=2 storage. This class is identical to IceModelVec2V, except that
539 //! components are not called `u` and `v` (to avoid confusion).
540 class IceModelVec2Stag : public IceModelVec2 {
541 public:
542 IceModelVec2Stag();
543 IceModelVec2Stag(IceGrid::ConstPtr grid, const std::string &short_name,
544 IceModelVecKind ghostedp, unsigned int stencil_width = 1);
545
546 typedef std::shared_ptr<IceModelVec2Stag> Ptr;
547 typedef std::shared_ptr<const IceModelVec2Stag> ConstPtr;
548
549 static Ptr ToStaggered(IceModelVec::Ptr input);
550
551 void create(IceGrid::ConstPtr grid, const std::string &short_name,
552 IceModelVecKind ghostedp, unsigned int stencil_width = 1);
553 virtual void staggered_to_regular(IceModelVec2S &result) const;
554 virtual void staggered_to_regular(IceModelVec2V &result) const;
555 virtual std::vector<double> absmaxcomponents() const;
556
557 //! Returns the values at interfaces of the cell i,j using the staggered grid.
558 /*! The ij member of the return value is set to 0, since it has no meaning in
559 this context.
560 */
561 inline StarStencil<double> star(int i, int j) const;
562 };
563
564 //! \brief A virtual class collecting methods common to ice and bedrock 3D
565 //! fields.
566 class IceModelVec3D : public IceModelVec {
567 public:
568 IceModelVec3D();
569 virtual ~IceModelVec3D();
570
571 void set_column(int i, int j, double c);
572 void set_column(int i, int j, const double *valsIN);
573 double* get_column(int i, int j);
574 const double* get_column(int i, int j) const;
575
576 // testing methods (for use from Python)
577 void set_column(int i, int j, const std::vector<double> &valsIN);
578 const std::vector<double> get_column_vector(int i, int j) const;
579
580 virtual double getValZ(int i, int j, double z) const;
581 virtual bool isLegalLevel(double z) const;
582
583 inline double& operator() (int i, int j, int k);
584 inline const double& operator() (int i, int j, int k) const;
585 protected:
586 void allocate(IceGrid::ConstPtr mygrid, const std::string &short_name,
587 IceModelVecKind ghostedp, const std::vector<double> &levels,
588 unsigned int stencil_width = 1);
589 private:
590 gsl_interp_accel *m_bsearch_accel;
591 };
592
593
594 //! Class for a 3d DA-based Vec for ice scalar quantities.
595 class IceModelVec3 : public IceModelVec3D {
596 public:
597 IceModelVec3();
598 IceModelVec3(IceGrid::ConstPtr mygrid, const std::string &short_name,
599 IceModelVecKind ghostedp,
600 unsigned int stencil_width = 1);
601
602 virtual ~IceModelVec3();
603
604 typedef std::shared_ptr<IceModelVec3> Ptr;
605 typedef std::shared_ptr<const IceModelVec3> ConstPtr;
606
607 static Ptr To3DScalar(IceModelVec::Ptr input);
608
609 void create(IceGrid::ConstPtr mygrid, const std::string &short_name,
610 IceModelVecKind ghostedp,
611 unsigned int stencil_width = 1);
612
613 void getHorSlice(Vec &gslice, double z) const; // used in iMmatlab.cc
614 void getHorSlice(IceModelVec2S &gslice, double z) const;
615 void getSurfaceValues(IceModelVec2S &gsurf, const IceModelVec2S &myH) const;
616
617 void sumColumns(IceModelVec2S &output, double A, double B) const;
618 };
619
620 /**
621 * Convert a PETSc Vec from the units in `from` into units in `to` (in place).
622 *
623 * @param v data to convert
624 * @param system unit system
625 * @param spec1 source unit specification string
626 * @param spec2 destination unit specification string
627 */
628 void convert_vec(Vec v, units::System::Ptr system,
629 const std::string &spec1, const std::string &spec2);
630
631 class IceModelVec2CellType;
632
633 /*!
634 * Average a scalar field from the staggered grid onto the regular grid by considering
635 * only ice-covered grid.
636 *
637 * If `include_floating_ice` is true, include floating ice, otherwise consider grounded
638 * icy cells only.
639 */
640 void staggered_to_regular(const IceModelVec2CellType &cell_type,
641 const IceModelVec2Stag &input,
642 bool include_floating_ice,
643 IceModelVec2S &result);
644
645 /*!
646 * Average a vector field from the staggered grid onto the regular grid by considering
647 * only ice-covered grid.
648 *
649 * If `include_floating_ice` is true, include floating ice, otherwise consider grounded
650 * icy cells only.
651 */
652 void staggered_to_regular(const IceModelVec2CellType &cell_type,
653 const IceModelVec2Stag &input,
654 bool include_floating_ice,
655 IceModelVec2V &result);
656
657 } // end of namespace pism
658
659 // include inline methods; contents are wrapped in namespace pism {...}
660 #include "IceModelVec_inline.hh"
661
662 #endif /* __IceModelVec_hh */
663