diff --git a/src/avalanche/functionObjects/gridfileWrite.H b/src/avalanche/functionObjects/gridfileWrite.H
index e30bbb73aba3b8af443409881edbcf8d99237951..635c58c2f8703eeb0dff00b9dde9a7e511ef1a0e 100644
--- a/src/avalanche/functionObjects/gridfileWrite.H
+++ b/src/avalanche/functionObjects/gridfileWrite.H
@@ -38,8 +38,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef gridfileWrite_H
-#define gridfileWrite_H
+#ifndef Foam_avalanche_gridfileWrite_H
+#define Foam_avalanche_gridfileWrite_H
 
 #include "regionFunctionObject.H"
 #include "wordRes.H"
diff --git a/src/avalanche/gistools/gridfile.C b/src/avalanche/gistools/gridfile.C
index c4d9239666255fb829dea5e6576f54193c92c689..cbbc76cf0609064812017920d57b26f72c16e549 100644
--- a/src/avalanche/gistools/gridfile.C
+++ b/src/avalanche/gistools/gridfile.C
@@ -36,19 +36,33 @@ Author
 #include <cmath>
 #include <regex>
 
+// Debugging output
+#undef  DebugInfo
+#define DebugInfo  if (debug_) this->log()
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 gridfile::gridfile()
+:
+    debug_(0),
+    v_(nullptr)
 {
-    debug = 0;
-    v_ = nullptr;
     clear();
 }
 
 
-gridfile::gridfile(const double &xllcenter, const double &yllcenter,
-                   const double &dx, const double &dy,
-                   const int &ncols, const int &nrows):
+gridfile::gridfile
+(
+    double xllcenter,
+    double yllcenter,
+    double dx,
+    double dy,
+    int ncols,
+    int nrows
+)
+:
+    debug_(0),
     xllcenter_(xllcenter),
     yllcenter_(yllcenter),
     dx_(dx),
@@ -58,7 +72,6 @@ gridfile::gridfile(const double &xllcenter, const double &yllcenter,
     NODATA_value_(-9999.),
     v_(nullptr)
 {
-    debug = 0;
     if (this->nrows_ > 0 && this->ncols_ > 0)
     {
         double * vv = new double [this->nrows_*this->ncols_];
@@ -90,11 +103,12 @@ gridfile::~gridfile()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-std::ostream &gridfile::log() const
+std::ostream& gridfile::log() const
 {
     return std::cout;
 }
 
+
 void gridfile::clear()
 {
     this->ncols_ = 0;
@@ -113,6 +127,7 @@ void gridfile::clear()
     v_ = nullptr;
 }
 
+
 std::string gridfile::info() const
 {
     std::ostringstream ss;
@@ -130,13 +145,13 @@ std::string gridfile::info() const
 }
 
 
-int gridfile::read(std::string filename)
+int gridfile::read(const std::string& fname)
 {
     clear();
 
-    this->filename_ = filename;
+    this->filename_ = fname;
 
-    std::ifstream pFile (filename);
+    std::ifstream pFile(filename_);
     if (!pFile.is_open())
     {
         return 0;
@@ -221,7 +236,7 @@ int gridfile::read(std::string filename)
     std::regex e("(\\S+)");
     while(!pFile.eof() && i < this->nrows_)
     {
-        if (debug) log() << "Reading line " << i << "/" << this->nrows_ << std::endl;
+        DebugInfo << "Reading line " << i << "/" << this->nrows_ << std::endl;
 
         std::getline(pFile, line);
         j = 0;
@@ -235,7 +250,7 @@ int gridfile::read(std::string filename)
             }
             catch(const std::invalid_argument &ia)
             {
-                if (debug) log() << "Invalid argument in record: " << ia.what() << std::endl;
+                DebugInfo << "Invalid argument in record: " << ia.what() << std::endl;
             }
             line = match.suffix().str();
         }
@@ -256,11 +271,12 @@ int gridfile::read(std::string filename)
     return 1;
 }
 
-int gridfile::write(std::string filename)
+
+int gridfile::write(const std::string& fname)
 {
-    this->filename_ = filename;
+    this->filename_ = fname;
 
-    std::ofstream pFile (filename);
+    std::ofstream pFile(filename_);
     if (!pFile.is_open())
     {
         return 0;
@@ -297,7 +313,7 @@ int gridfile::write(std::string filename)
     return 1;
 }
 
-double gridfile::interpolate(const double &x, const double &y) const
+double gridfile::interpolate(const double x, const double y) const
 {
     double dx = x-this->xllcenter_;
     double dy = y-this->yllcenter_;
@@ -324,7 +340,8 @@ double gridfile::interpolate(const double &x, const double &y) const
            v_[jpp][ipp]*(dy)*(dx);
 }
 
-double gridfile::interpolateNN(const double &x, const double &y) const
+
+double gridfile::interpolateNN(const double x, const double y) const
 {
     double dx = x-this->xllcenter_;
     double dy = y-this->yllcenter_;
diff --git a/src/avalanche/gistools/gridfile.H b/src/avalanche/gistools/gridfile.H
index b6658217faee2adc776eb65976f24b943625c6a4..a4b738f9e10b1160d7953e41d4db055796182a98 100644
--- a/src/avalanche/gistools/gridfile.H
+++ b/src/avalanche/gistools/gridfile.H
@@ -37,8 +37,8 @@ Author
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef gridfile_H
-#define gridfile_H
+#ifndef Foam_avalanche_gridfile_H
+#define Foam_avalanche_gridfile_H
 
 #include <iostream>
 
@@ -48,56 +48,60 @@ public:
 
     gridfile();
 
-    gridfile(const double &xllcenter, const double &yllcenter,
-             const double &dx, const double &dy,
-             const int &ncols, const int &nrows);
+    gridfile(double xllcenter, double yllcenter,
+             double dx, double dy,
+             int ncols,  int nrows);
 
     ~gridfile();
 
-    std::ostream &log() const;
+    // Logging (eg, stdout)
+    std::ostream& log() const;
 
     void clear();
 
     std::string info() const;
 
-    int read(std::string filename_);
+    int read(const std::string& fname);
 
-    int write(std::string filename_);
+    int write(const std::string& fname);
 
     //Bilinear interpolation
-    double interpolate(const double &x, const double &y) const;
+    double interpolate(const double x, const double y) const;
 
     //Nearest neighbor interpolation
-    double interpolateNN(const double &x, const double &y) const;
+    double interpolateNN(const double x, const double y) const;
 
-    inline const std::string &filename() const {return this->filename_;}
+    const std::string& filename() const noexcept { return filename_; }
 
-    inline const double &xllcenter() const {return xllcenter_;}
+    double xllcenter() const noexcept { return xllcenter_; }
 
-    inline const double &yllcenter() const {return yllcenter_;}
+    double yllcenter() const noexcept { return yllcenter_; }
 
-    inline const double &dx() const {return dx_;}
+    double dx() const noexcept { return dx_; }
 
-    inline const double &dy() const {return dy_;}
+    double dy() const noexcept { return dy_; }
 
-    inline const unsigned int &ncols() const {return ncols_;}
+    unsigned int ncols() const noexcept {return ncols_;}
 
-    inline const unsigned int &nrows() const {return nrows_;}
+    unsigned int nrows() const noexcept { return nrows_; }
 
-    inline const unsigned int &xCellCount() const {return ncols();}
+    unsigned int xCellCount() const noexcept { return ncols_; }
 
-    inline const unsigned int &yCellCount() const {return nrows();}
+    unsigned int yCellCount() const noexcept { return nrows_ ; }
 
-    inline const double &v(int xindex, int yindex) const {return v_[xindex][yindex];}
+    double v(int xindex, int yindex) const { return v_[xindex][yindex]; }
 
-    inline double &vRef(int xindex, int yindex) {return v_[xindex][yindex];}
+    double& vRef(int xindex, int yindex) { return v_[xindex][yindex]; }
 
-    inline double** vRef() {return v_;}
+    inline double** vRef() { return v_; }
+
+    double NODATA_value() const noexcept { return NODATA_value_; }
 
-    inline const double &NODATA_value() const {return NODATA_value_;}
 
 private:
-    int debug;
+
+    // Debug level
+    int debug_;
 
     //Last known filename
     std::string filename_;
diff --git a/src/avalanche/gistools/shapefile.C b/src/avalanche/gistools/shapefile.C
index d0fda0440c1f1809e98877b344a1bc7e8b180f8f..8de4200aeea30cd9f69749cc981d30705d6b316f 100644
--- a/src/avalanche/gistools/shapefile.C
+++ b/src/avalanche/gistools/shapefile.C
@@ -39,6 +39,11 @@ Author
 #include <cstring>
 #include <stdexcept>
 
+// Debugging output
+#undef  DebugInfo
+#define DebugInfo  if (debug_) this->log()
+
+
 // * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
 
 namespace
@@ -53,7 +58,7 @@ namespace
 #undef ntohl
 #endif
 
-static inline uint32_t ntohl(uint32_t u)
+inline uint32_t ntohl(uint32_t u)
 {
 #if (__BYTE_ORDER__ == __ORDER_BIG_ENDIAN__)
     return u;
@@ -76,7 +81,7 @@ static inline uint32_t ntohl(uint32_t u)
 
 
 template<class T>
-static std::string vec2string(const std::vector<T>& v)
+std::string vec2string(const std::vector<T>& v)
 {
     std::ostringstream ss;
 
@@ -92,6 +97,7 @@ static std::string vec2string(const std::vector<T>& v)
 
     return ss.str();
 }
+
 } // End anonymous namespace
 
 
@@ -122,12 +128,12 @@ int shapefile::addRecord()
     this->shapeparts_.push_back(0);
     this->shapepoints_.push_back(0);
 
-    this->px_.push_back(std::vector<double>());
-    this->py_.push_back(std::vector<double>());
-    this->pz_.push_back(std::vector<double>());
-    this->pm_.push_back(std::vector<double>());
-    this->pn_.push_back(std::vector<int>());
-    this->pt_.push_back(std::vector<int>());
+    this->px_.emplace_back();
+    this->py_.emplace_back();
+    this->pz_.emplace_back();
+    this->pm_.emplace_back();
+    this->pn_.emplace_back();
+    this->pt_.emplace_back();
 
     this->bxmin_.push_back(0);
     this->bxmax_.push_back(1);
@@ -138,8 +144,8 @@ int shapefile::addRecord()
     this->bmmin_.push_back(0);
     this->bmmax_.push_back(1);
 
-    this->v_.push_back(std::vector<double>());
-    this->c_.push_back(std::vector<std::string>());
+    this->v_.emplace_back();
+    this->c_.emplace_back();
 
     return this->recordcount_-1;
 }
@@ -192,7 +198,7 @@ int shapefile::addField(std::string fieldname, char fieldtype, int length, int d
         this->isnumeric_.push_back(false);
         for (std::vector<std::string> &cRef : c_)
         {
-            cRef.push_back(std::string());
+            cRef.emplace_back();
         }
         this->fieldIndex_.push_back(this->c_[0].size()-1);
     }
@@ -216,7 +222,7 @@ void shapefile::setField(int recordIndex, int fieldIndex, const std::string &val
 }
 
 
-const double &shapefile::numericfield(int recordIndex, int fieldIndex) const
+double shapefile::numericfield(int recordIndex, int fieldIndex) const
 {
     if(!this->isnumeric_[fieldIndex])
         throw std::invalid_argument("field is not numeric");
@@ -224,7 +230,7 @@ const double &shapefile::numericfield(int recordIndex, int fieldIndex) const
 }
 
 
-const std::string &shapefile::stringfield(int recordIndex, int fieldIndex) const
+const std::string& shapefile::stringfield(int recordIndex, int fieldIndex) const
 {
     if(this->isnumeric_[fieldIndex])
         throw std::invalid_argument("field is numeric");
@@ -232,7 +238,7 @@ const std::string &shapefile::stringfield(int recordIndex, int fieldIndex) const
 }
 
 
-std::ostream &shapefile::log() const
+std::ostream& shapefile::log() const
 {
     return std::cout;
 }
@@ -240,9 +246,9 @@ std::ostream &shapefile::log() const
 
 void shapefile::clear()
 {
-    if (debug_) log() << "shapefile cleared." << std::endl;
+    DebugInfo << "shapefile cleared." << std::endl;
 
-    filename_ = std::string("");
+    filename_.clear();
 
     shptype_ = NULLSHP;
     recordcount_ = 0;
@@ -336,7 +342,7 @@ void shapefile::calcBoundingBox()
         }
     }
 
-    if (debug_) log() << "BB ("<<  xmin_ << ", " << ymin_ << ") -- (" << xmax_ << ", " << ymax_ << ")" << std::endl;
+    DebugInfo <<  xmin_ << ", " << ymin_ << ") -- (" << xmax_ << ", " << ymax_ << ")" << std::endl;
 }
 
 
@@ -404,7 +410,7 @@ int shapefile::read(std::string filename)
     char* c = &buffer[0];
 
     int bufferlength = buffer.size();
-    if (debug_) log() << "\n\nShx-file is " << bufferlength << " byte long." << std::endl;
+    DebugInfo << "\n\nShx-file is " << bufferlength << " byte long." << std::endl;
 
     int pos = 100;
 
@@ -412,13 +418,12 @@ int shapefile::read(std::string filename)
     {
         uint32_t recordoffset = ntohl(*reinterpret_cast<uint32_t*>(c+pos));
         uint32_t recordlength = ntohl(*reinterpret_cast<uint32_t*>(c+pos+4));
-
-        if (debug_) log() << "Recordoffset: " << recordoffset << std::endl;
-        if (debug_) log() << "Recordlength: " << recordlength << std::endl;
-
         pos += 8;
 
-        if (debug_) log() << "read " << pos << "/" << bufferlength << std::endl;
+        DebugInfo
+            << "Recordoffset: " << recordoffset << std::endl
+            << "Recordlength: " << recordlength << std::endl
+            << "read " << pos << "/" << bufferlength << std::endl;
 
         recordPositions.push_back(recordoffset);
     }
@@ -436,7 +441,7 @@ int shapefile::read(std::string filename)
 
     bufferlength = buffer.size();
 
-    if (debug_) log() << "\n\nShp-file is " << bufferlength << " byte long." << std::endl;
+    DebugInfo << "\n\nShp-file is " << bufferlength << " byte long." << std::endl;
 
     uint32_t filecode = ntohl(*reinterpret_cast<uint32_t*>(c));
     uint32_t filelength = ntohl(*reinterpret_cast<uint32_t*>(c+24));
@@ -452,15 +457,17 @@ int shapefile::read(std::string filename)
     this->mmin_ = *reinterpret_cast<double*>(c+84);
     this->mmax_ = *reinterpret_cast<double*>(c+92);
 
-    if (debug_) log() << "Filecode: " << filecode << std::endl;
-    if (debug_) log() << "Filelength: " << filelength << std::endl;
-    if (debug_) log() << "Version: " << version << std::endl;
-    if (debug_) log() << "Shypetype: " << shptype_ << std::endl;
-    if (debug_) log() << "BB ("
-                     << this->xmin_ << ", " << this->ymin_ << ", "
-                     << this->zmin_ << ", " << this->mmin_ << ") -- ("
-                     << this->xmax_ << ", " << this->ymax_ << ", "
-                     << this->zmin_ << ", " << this->mmin_ << ")" << std::endl;
+    DebugInfo
+        << "Filecode: " << filecode << std::endl
+        << "Filelength: " << filelength << std::endl
+        << "Version: " << version << std::endl
+        << "Shypetype: " << shptype_ << std::endl
+        << "BB ("
+            << this->xmin_ << ", " << this->ymin_ << ", "
+            << this->zmin_ << ", " << this->mmin_ << ") -- ("
+            << this->xmax_ << ", " << this->ymax_ << ", "
+            << this->zmin_ << ", " << this->mmin_ << ")"
+        << std::endl;
 
     unsigned int i = 0;
     pos = 100;
@@ -479,25 +486,26 @@ int shapefile::read(std::string filename)
         uint32_t shptype2 = *reinterpret_cast<uint32_t*>(pos2);
         pos2 += 4;
 
-        if (debug_) log() << "Shapetype (global/local): " << this->shptype_ << "/" << shptype2 << std::endl;
+        DebugInfo << "Shapetype (global/local): " << this->shptype_ << "/" << shptype2 << std::endl;
         if (this->shptype_ != shptype2)
         {
             log() << "Uexpected shapetype. Ignoring local shapetype";
         }
 
-        if (debug_) log() << std::endl << "Recordnumber: " << recordnumber << std::endl;
-        if (debug_) log() << "Recordoffset: " << pos/2 << std::endl;
-        if (debug_) log() << "Recordlength: " << recordlength << std::endl;
-
-        if (debug_) log() << "Current pos is " << pos << std::endl;
+        DebugInfo
+            << std::endl
+            << "Recordnumber: " << recordnumber << std::endl
+            << "Recordoffset: " << pos/2 << std::endl
+            << "Recordlength: " << recordlength << std::endl
+            << "Current pos is " << pos << std::endl;
 
         if (recordlength < 3){
 
-            if (debug_) log() << "Recording is empty - not reading further." << std::endl;
+            DebugInfo << "Recording is empty - not reading further." << std::endl;
 
-            this->px_.push_back(std::vector<double>());
-            this->py_.push_back(std::vector<double>());
-            this->pn_.push_back(std::vector<int>());
+            this->px_.emplace_back();
+            this->py_.emplace_back();
+            this->pn_.emplace_back();
             this->shapeparts_.push_back(0);
             this->shapepoints_.push_back(0);
 
@@ -528,7 +536,7 @@ int shapefile::read(std::string filename)
             double yy = *reinterpret_cast<double*>(pos2);
             pos2 += 8;
 
-            if (debug_) log() << "p(" << xx << ", " << yy << ")" << std::endl;
+            DebugInfo << "p(" << xx << ", " << yy << ")" << std::endl;
 
             x.push_back(xx);
             y.push_back(yy);
@@ -538,7 +546,7 @@ int shapefile::read(std::string filename)
                 double zz = *reinterpret_cast<double*>(pos2);
                 pos2 += 8;
 
-                if (debug_) log() << ", z " << zz << std::endl;
+                DebugInfo << ", z " << zz << std::endl;
                 z.push_back(yy);
             }
 
@@ -548,7 +556,7 @@ int shapefile::read(std::string filename)
                 double mm = *reinterpret_cast<double*>(pos2);
                 pos2 += 8;
 
-                if (debug_) log() << "m " << mm << std::endl;
+                DebugInfo << "m " << mm << std::endl;
                 m.push_back(mm);
             }
         }
@@ -568,7 +576,7 @@ int shapefile::read(std::string filename)
             uint32_t numpoints = *reinterpret_cast<uint32_t*>(pos2);
             pos2 += 4;
 
-            if (debug_) log() << "# Point: " << numpoints << std::endl;
+            DebugInfo << "# Point: " << numpoints << std::endl;
 
             this->shapeparts_.push_back(1);
             this->shapepoints_.push_back(numpoints);
@@ -580,7 +588,7 @@ int shapefile::read(std::string filename)
                 double yy = *reinterpret_cast<double*>(pos2);
                 pos2 += 8;
 
-                if (debug_) log() << "p(" << xx << ", " << yy << ")" << std::endl;
+                DebugInfo << "p(" << xx << ", " << yy << ")" << std::endl;
 
                 x.push_back(xx);
                 y.push_back(yy);
@@ -598,7 +606,7 @@ int shapefile::read(std::string filename)
                     double zz = *reinterpret_cast<double*>(pos2);
                     pos2 += 8;
 
-                    if (debug_) log() << "z " << zz << std::endl;
+                    DebugInfo << "z " << zz << std::endl;
 
                     z.push_back(zz);
                 }
@@ -617,13 +625,11 @@ int shapefile::read(std::string filename)
                     double mm = *reinterpret_cast<double*>(pos2);
                     pos2 += 8;
 
-                    if (debug_) log() << "m " << mm << std::endl;
+                    DebugInfo << "m " << mm << std::endl;
 
                     m.push_back(mm);
                 }
             }
-
-
         }
         else if (this->shptype_ == POLYLINE ||
                  this->shptype_ == POLYLINEM ||
@@ -650,8 +656,9 @@ int shapefile::read(std::string filename)
             uint32_t numpoints = *reinterpret_cast<uint32_t*>(pos2);
             pos2 += 4;
 
-            if (debug_) log() << "# Parts: " << numparts << std::endl;
-            if (debug_) log() << "# Point: " << numpoints << std::endl;
+            DebugInfo
+                << "# Parts: " << numparts << std::endl
+                << "# Point: " << numpoints << std::endl;
 
             this->shapeparts_.push_back(numparts);
             this->shapepoints_.push_back(numpoints);
@@ -661,7 +668,7 @@ int shapefile::read(std::string filename)
                 int j = *reinterpret_cast<uint32_t*>(pos2);
                 pos2 += 4;
 
-                if (debug_) log() << "Part start: " << j << std::endl;
+                DebugInfo << "Part start: " << j << std::endl;
 
                 parts.push_back(j);
             }
@@ -673,7 +680,7 @@ int shapefile::read(std::string filename)
                     int j = *reinterpret_cast<uint32_t*>(pos2);
                     pos2 += 4;
 
-                    if (debug_) log() << "Part type: " << j << std::endl;
+                    DebugInfo << "Part type: " << j << std::endl;
 
                     parttypes.push_back(j);
                 }
@@ -686,7 +693,7 @@ int shapefile::read(std::string filename)
                 double yy = *reinterpret_cast<double*>(pos2);
                 pos2 += 8;
 
-                if (debug_) log() << "p(" << xx << ", " << yy << ")" << std::endl;
+                DebugInfo << "p(" << xx << ", " << yy << ")" << std::endl;
 
                 x.push_back(xx);
                 y.push_back(yy);
@@ -706,7 +713,7 @@ int shapefile::read(std::string filename)
                     double zz = *reinterpret_cast<double*>(pos2);
                     pos2 += 8;
 
-                    if (debug_) log() << "z " << zz << std::endl;
+                    DebugInfo << "z " << zz << std::endl;
 
                     z.push_back(zz);
                 }
@@ -727,7 +734,7 @@ int shapefile::read(std::string filename)
                     double mm = *reinterpret_cast<double*>(pos2);
                     pos2 += 8;
 
-                    if (debug_) log() << "m " << mm << std::endl;
+                    DebugInfo << "m " << mm << std::endl;
 
                     m.push_back(mm);
                 }
@@ -736,7 +743,7 @@ int shapefile::read(std::string filename)
         else
         {
             //Not implemented shape
-            if (debug_) log() << "Shape type \"" << shptype_ << "\" not implemented - not reading further." << std::endl;
+            DebugInfo << "Shape type \"" << shptype_ << "\" not implemented - not reading further." << std::endl;
             return 1;
         }
 
@@ -753,12 +760,12 @@ int shapefile::read(std::string filename)
             break;
         }
 
-        if (debug_) log() << "read " << pos << "/" << bufferlength << std::endl;
+        DebugInfo << "read " << pos << "/" << bufferlength << std::endl;
     }
 
     this->recordcount_ = recordPositions.size();
 
-    if (debug_) log() << "\nFinished reading shp-file. Found " << this->recordcount_ << ". Closing file.\n" << std::endl;
+    DebugInfo << "\nFinished reading shp-file. Found " << this->recordcount_ << ". Closing file.\n" << std::endl;
 
     input.close();
 
@@ -772,7 +779,7 @@ int shapefile::read(std::string filename)
     c = &buffer[0];
 
     bufferlength = buffer.size();
-    if (debug_) log() << "\n\nDbf-file is " << bufferlength << " byte long." << std::endl;
+    DebugInfo << "\n\nDbf-file is " << bufferlength << " byte long." << std::endl;
 
     this->dbase_versionbyte_ = c[0];
     this->dbase_datestring_[0] = c[1];
@@ -783,9 +790,10 @@ int shapefile::read(std::string filename)
     uint16_t headerlength = *reinterpret_cast<uint16_t*>(c+8);
     uint16_t recordlength = *reinterpret_cast<uint16_t*>(c+10);
 
-    if (debug_) log() << "Recordnumber: " << recordnumber << std::endl;
-    if (debug_) log() << "Headerlength: " << headerlength << std::endl;
-    if (debug_) log() << "Recordlenght: " << recordlength << std::endl;
+    DebugInfo
+        << "Recordnumber: " << recordnumber << std::endl
+        << "Headerlength: " << headerlength << std::endl
+        << "Recordlenght: " << recordlength << std::endl;
 
 
     i = 0;
@@ -800,21 +808,23 @@ int shapefile::read(std::string filename)
 
         strncpy(fieldname, c+pos, 11);
         strncpy(fieldtype, c+pos+11, 1);
-        if (debug_) log() << "\n\n\nFieldname: " << fieldname << std::endl;
-        if (debug_) log() << "Fieldtype: " << fieldtype << std::endl;
+        DebugInfo
+            << "\n\n\nFieldname: " << fieldname << std::endl
+            << "Fieldtype: " << fieldtype << std::endl;
 
 
         uint32_t fieldpos = *reinterpret_cast<uint32_t*>(c+pos+12);
         uint16_t fieldlength = *reinterpret_cast<uint8_t*>(c+pos+16);
         uint16_t fielddec = *reinterpret_cast<uint8_t*>(c+pos+17);
 
-        if (debug_) log() << "Fieldpos: " << fieldpos << std::endl;
-        if (debug_) log() << "Fieldlength: " << fieldlength << std::endl;
-        if (debug_) log() << "Fielddec: " << fielddec << std::endl;
+        DebugInfo
+            << "Fieldpos: " << fieldpos << std::endl
+            << "Fieldlength: " << fieldlength << std::endl
+            << "Fielddec: " << fielddec << std::endl;
 
 
         uint16_t termbit = *reinterpret_cast<uint8_t*>(c+pos+32);
-        if (debug_) log() << "termbit: " << termbit << std::endl;
+        DebugInfo << "termbit: " << termbit << std::endl;
 
         pos += 32;
 
@@ -841,7 +851,9 @@ int shapefile::read(std::string filename)
 
     this->fieldcount_ = i;
 
-    if (debug_) log() << std::endl << std::endl << "Fieldcount: " << fieldcount_ << std::endl;
+    DebugInfo
+        << std::endl << std::endl
+        << "Fieldcount: " << fieldcount_ << std::endl;
 
 
     pos = headerlength;
@@ -854,8 +866,9 @@ int shapefile::read(std::string filename)
         std::vector<double> values;
         std::vector<std::string> strings;
 
-        if (debug_) log() << std::endl;
-        if (debug_) log() << "Reading record #" << i << std::endl << std::endl;
+        DebugInfo
+            << std::endl
+            << "Reading record #" << i << std::endl << std::endl;
 
         for(unsigned int j=0; j<fieldcount_; j++)
         {
@@ -881,32 +894,32 @@ int shapefile::read(std::string filename)
                 }
                 catch(const std::invalid_argument &ia)
                 {
-                    if (debug_) log() << "Invalid argument in record: " << ia.what() << std::endl;
+                    DebugInfo << "Invalid argument in record: " << ia.what() << std::endl;
                 }
-                strings.push_back(std::string(""));
+                strings.emplace_back();
                 values.push_back(val);
 
-                if (debug_) log() << fn_[j] << " = " << "\"" << buf << "\"" << "/" << values[values.size()-1] << std::endl;
+                DebugInfo << fn_[j] << " = " << "\"" << buf << "\"" << "/" << values[values.size()-1] << std::endl;
             }
             else if (this->ft_[j] == 'C')
             {
-                strings.push_back(std::string(buf));
+                strings.emplace_back(buf);
                 values.push_back(NAN);
 
-                if (debug_) log() << fn_[j] << " = " << "\"" << buf << "\"" << "/" << strings[strings.size()-1] << std::endl;
+                DebugInfo << fn_[j] << " = " << "\"" << buf << "\"" << "/" << strings[strings.size()-1] << std::endl;
             }
             else if (this->ft_[j] == 'D')
             {
                 //Date in format YYYYMMDD
-                strings.push_back(std::string(buf));
+                strings.emplace_back(buf);
                 values.push_back(NAN);
 
-                if (debug_) log() << fn_[j] << " = " << "\"" << buf << "\"" << "/" << strings[strings.size()-1] << std::endl;
+                DebugInfo << fn_[j] << " = " << "\"" << buf << "\"" << "/" << strings[strings.size()-1] << std::endl;
             }
             else if (this->ft_[j] == 'L')
             {
                 //Boolean value YyTt -> 1 NnFf -> 0
-                strings.push_back(std::string(buf));
+                strings.emplace_back(buf);
                 if (buf[0] == 'Y' || buf[0] == 'y' || buf[0] == 'T' || buf[0] == 't')
                     values.push_back(1);
                 else
@@ -914,7 +927,7 @@ int shapefile::read(std::string filename)
             }
             else
             {
-                if (debug_) log() << "Value type \"" << this->ft_[j] << "\" not implemented - not reading further." << std::endl;
+                DebugInfo << "Value type \"" << this->ft_[j] << "\" not implemented - not reading further." << std::endl;
                 return 0;
             }
         }
@@ -923,14 +936,14 @@ int shapefile::read(std::string filename)
         this->c_.push_back(strings);
         i++;
 
-        if (debug_) log() << "read " << pos << "/" << bufferlength << std::endl;
+        DebugInfo << "read " << pos << "/" << bufferlength << std::endl;
 
     }
     while(i <= recordcount_ && pos < bufferlength);
 
     input.close();
 
-    if (debug_) log() << "Finished reading file." << std::endl;
+    DebugInfo << "Finished reading file." << std::endl;
     return 1;
 }
 
@@ -1032,7 +1045,7 @@ int shapefile::write(std::string filename)
         }
         else
         {
-            if (debug_) log() << "Shape type \"" << this->shptype_ << "\" not implemented - not writing further." << std::endl;
+            DebugInfo << "Shape type \"" << this->shptype_ << "\" not implemented - not writing further." << std::endl;
             return 0;
         }
     }
@@ -1209,7 +1222,7 @@ int shapefile::write(std::string filename)
         }
         else
         {
-            if (debug_) log() << "Shape type \"" << this->shptype_ << "\" not implemented - not writing further." << std::endl;
+            DebugInfo << "Shape type \"" << this->shptype_ << "\" not implemented - not writing further." << std::endl;
             return 0;
         }
 
@@ -1239,8 +1252,8 @@ int shapefile::write(std::string filename)
     output.write(reinterpret_cast<char*>(&this->mmax_), 8); //boundingbox
 
     for (unsigned int i=0; i<recordLengths.size(); i++)
-   {
-       if (debug_) log() << "Writing indices: " << recordPos[i] << "/" << recordLengths[i] << std::endl;
+    {
+        DebugInfo << "Writing indices: " << recordPos[i] << "/" << recordLengths[i] << std::endl;
 
        recordPos[i] = ntohl(recordPos[i]);
        recordLengths[i] = ntohl(recordLengths[i]);
@@ -1327,8 +1340,9 @@ int shapefile::write(std::string filename)
                     snprintf(buf, buflen+1, formatstr, int(v_[i][j]));
                 }
 
-                if (debug_) log() << formatstr << std::endl;
-                if (debug_) log() << v_[i][j] << "->" << "\"" << buf << "\"" << std::endl;
+                DebugInfo
+                    << formatstr << std::endl
+                    << v_[i][j] << "->" << "\"" << buf << "\"" << std::endl;
 
                 /*for(int k=0; k<buflen; k++)
                 {
@@ -1337,7 +1351,7 @@ int shapefile::write(std::string filename)
                     if (buf[k] == '\0')
                         buf[k] = '0';
                 }*/
-                if (debug_) log() << v_[i][j] << "->" << "\"" << buf << "\"" << std::endl;
+                DebugInfo << v_[i][j] << "->" << "\"" << buf << "\"" << std::endl;
 
                 output.write(buf, buflen); //write all values
             }
@@ -1350,14 +1364,16 @@ int shapefile::write(std::string filename)
                 char* buf = &(bufstring[0]);
                 //char buf[buflen+1] = {'\0'};
 
-                if (debug_) log() << c_[i].size() << std::endl;
-                if (debug_) log() << c_[i][j] << std::endl;
+                DebugInfo
+                    << c_[i].size() << std::endl
+                    << c_[i][j] << std::endl;
+
                 strncpy(buf, c_[i][j].c_str(), this->fl_[j]);
                 output.write(buf, buflen); //write all values
             }
             else
             {
-                if (debug_) log() << "Value type \"" << this->ft_[j] << "\" not implemented - not reading further." << std::endl;
+                DebugInfo << "Value type \"" << this->ft_[j] << "\" not implemented - not reading further." << std::endl;
                 return 0;
             }
         }
diff --git a/src/avalanche/gistools/shapefile.H b/src/avalanche/gistools/shapefile.H
index c27d3021a8ea0ebc50c153fc2fccf37e4ad0cc8b..36800d05dfd009640a0e75fe0a450a5316182e06 100644
--- a/src/avalanche/gistools/shapefile.H
+++ b/src/avalanche/gistools/shapefile.H
@@ -29,7 +29,7 @@ Class
 Description
     A minimalistic class for handling ESRI shapefiles. The class can read and
     write a subset of ESRI shapefiles. The class is limited to points,
-    polylines and ploygons. 
+    polylines and ploygons.
 
     en.wikipedia.org/wiki/Shapefile
 
@@ -43,9 +43,10 @@ Author
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef shapefile_H
-#define shapefile_H
+#ifndef Foam_avalanche_shapefile_H
+#define Foam_avalanche_shapefile_H
 
+#include <cstdint>  // For uint32_t
 #include <vector>
 #include <iostream>
 
@@ -83,47 +84,46 @@ public:
 
     shapefile(type shptype_ = NULLSHP);
 
-    inline const uint32_t &shptype() const {return this->shptype_;}
+    uint32_t shptype() const noexcept { return shptype_; }
 
-    inline const unsigned int &fieldcount() const {return this->fieldcount_;}
+    unsigned int fieldcount() const noexcept { return fieldcount_; }
 
-    inline const unsigned int &recordcount() const {return this->recordcount_;}
+    unsigned int recordcount() const noexcept { return recordcount_; }
 
-    inline const int &partcount(int recordIndex) const {return this->shapeparts_[recordIndex];}
+    int partcount(int recordIndex) const { return shapeparts_[recordIndex]; }
 
-    inline const int &pointcount(int recordIndex) const {return this->shapepoints_[recordIndex];}
+    int pointcount(int recordIndex) const { return shapepoints_[recordIndex]; }
 
     int pointcount(int recordIndex, int partIndex) const;
 
-    inline const int &parttype(int recordIndex, int partIndex) const {return this->pt_[recordIndex][partIndex];}
+    int parttype(int recordIndex, int partIndex) const { return pt_[recordIndex][partIndex]; }
 
-    inline const int &partStart(int recordIndex, int partIndex) const {return this->pn_[recordIndex][partIndex];}
+    int partStart(int recordIndex, int partIndex) const { return pn_[recordIndex][partIndex]; }
 
-    inline const double &x(int recordIndex, int pointIndex) const {return this->px_[recordIndex][pointIndex];}
+    double x(int recordIndex, int pointIndex) const { return px_[recordIndex][pointIndex]; }
 
-    inline const double &y(int recordIndex, int pointIndex) const {return this->py_[recordIndex][pointIndex];}
+    double y(int recordIndex, int pointIndex) const { return py_[recordIndex][pointIndex]; }
 
-    inline const double &z(int recordIndex, int pointIndex) const {return this->pz_[recordIndex][pointIndex];}
+    double z(int recordIndex, int pointIndex) const { return pz_[recordIndex][pointIndex]; }
 
-    inline const double &m(int recordIndex, int pointIndex) const {return this->pm_[recordIndex][pointIndex];}
+    double m(int recordIndex, int pointIndex) const { return pm_[recordIndex][pointIndex]; }
 
-    inline double &xRef(int recordIndex, int pointIndex) {return this->px_[recordIndex][pointIndex];}
+    double xRef(int recordIndex, int pointIndex) { return px_[recordIndex][pointIndex]; }
 
-    inline double &yRef(int recordIndex, int pointIndex) {return this->py_[recordIndex][pointIndex];}
+    double yRef(int recordIndex, int pointIndex) { return py_[recordIndex][pointIndex]; }
 
-    inline double &zRef(int recordIndex, int pointIndex) {return this->pz_[recordIndex][pointIndex];}
+    double zRef(int recordIndex, int pointIndex) { return pz_[recordIndex][pointIndex]; }
 
-    inline double &mRef(int recordIndex, int pointIndex) {return this->pm_[recordIndex][pointIndex];}
+    double mRef(int recordIndex, int pointIndex) { return pm_[recordIndex][pointIndex]; }
 
 
-    inline const std::string &fieldname(int fieldIndex) const {return this->fn_[fieldIndex];}
+    const std::string& fieldname(int fieldIndex) const { return fn_[fieldIndex]; }
 
-    inline bool isnumeric(int fieldIndex) const {return this->isnumeric_[fieldIndex];}
+    bool isnumeric(int fieldIndex) const { return isnumeric_[fieldIndex]; }
 
+    double numericfield(int recordIndex, int fieldIndex_) const;
 
-    const double &numericfield(int recordIndex, int fieldIndex_) const;
-
-    const std::string &stringfield(int recordIndex, int fieldIndex_) const;
+    const std::string& stringfield(int recordIndex, int fieldIndex_) const;
 
 
     int addRecord();
@@ -139,7 +139,7 @@ public:
     void setField(int recordIndex, int fieldIndex_, const std::string &value);
 
 
-    std::ostream &log() const;
+    std::ostream& log() const;
 
     void clear();
 
@@ -171,25 +171,25 @@ private:
     std::vector<int> shapeparts_;
 
     //Index of parts
-    std::vector<std::vector<int> > pn_;
+    std::vector<std::vector<int>> pn_;
 
     //Type of parts, only used in multipatch
-    std::vector<std::vector<int> > pt_;
+    std::vector<std::vector<int>> pt_;
 
     //Number of points in a record
     std::vector<int> shapepoints_;
 
     //x-values
-    std::vector<std::vector<double> > px_;
+    std::vector<std::vector<double>> px_;
 
     //y-values
-    std::vector<std::vector<double> > py_;
+    std::vector<std::vector<double>> py_;
 
     //z-values
-    std::vector<std::vector<double> > pz_;
+    std::vector<std::vector<double>> pz_;
 
     //m-values
-    std::vector<std::vector<double> > pm_;
+    std::vector<std::vector<double>> pm_;
 
     //Number of different fields
     unsigned int fieldcount_;
@@ -213,14 +213,14 @@ private:
     std::vector<int> fieldIndex_;
 
     //Field values
-    std::vector<std::vector<double> > v_;
-    std::vector<std::vector<std::string> > c_;
+    std::vector<std::vector<double>> v_;
+    std::vector<std::vector<std::string>> c_;
 
     //part boundary boxes;
-     std::vector<double> bxmin_, bxmax_;
-     std::vector<double> bymin_, bymax_;
-     std::vector<double> bzmin_, bzmax_;
-     std::vector<double> bmmin_, bmmax_;
+    std::vector<double> bxmin_, bxmax_;
+    std::vector<double> bymin_, bymax_;
+    std::vector<double> bzmin_, bzmax_;
+    std::vector<double> bmmin_, bmmax_;
 
     //Boundary box
     double xmin_, xmax_;