.. _program_listing_file_src_modules_PBSM3D.hpp:

Program Listing for File PBSM3D.hpp
===================================

|exhale_lsh| :ref:`Return to documentation for file <file_src_modules_PBSM3D.hpp>` (``src/modules/PBSM3D.hpp``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: cpp

   //
   // Canadian Hydrological Model - The Canadian Hydrological Model (CHM) is a novel
   // modular unstructured mesh based approach for hydrological modelling
   // Copyright (C) 2018 Christopher Marsh
   //
   // This file is part of Canadian Hydrological Model.
   //
   // Canadian Hydrological Model is free software: you can redistribute it and/or
   // modify
   // it under the terms of the GNU General Public License as published by
   // the Free Software Foundation, either version 3 of the License, or
   // (at your option) any later version.
   //
   // Canadian Hydrological Model is distributed in the hope that it will be useful,
   // but WITHOUT ANY WARRANTY; without even the implied warranty of
   // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   // GNU General Public License for more details.
   //
   // You should have received a copy of the GNU General Public License
   // along with Canadian Hydrological Model.  If not, see
   // <http://www.gnu.org/licenses/>.
   //
   
   #pragma once
   
   //#define BOOST_MATH_INSTRUMENT
   #include "interpolation.hpp"
   #include "logger.hpp"
   #include "module_base.hpp"
   #include "triangulation.hpp"
   
   #include "math/coordinates.hpp"
   #include "math/LinearAlgebra.hpp"
   
   #include <physics/PhysConst.h>
   #include "physics/Atmosphere.h"
   
   #include <meteoio/MeteoIO.h>
   
   #include <cmath>
   #include <gsl/gsl_cdf.h>
   #include <gsl/gsl_integration.h>
   #include <gsl/gsl_math.h>
   #include <gsl/gsl_randist.h>
   #include <gsl/gsl_sf_lambert.h>
   #include <vector>
   
   #include <armadillo>
   
   #include <cstdlib>
   #include <string>
   //#define _USE_MATH_DEFINES
   //#include <math.h>
   
   #include <boost/math/tools/roots.hpp>
   #include <boost/math/tools/tuple.hpp>
   
   
   class PBSM3D : public module_base
   {
       REGISTER_MODULE_HPP(PBSM3D);
   
     public:
       PBSM3D(config_file cfg);
       ~PBSM3D();
       void run(mesh& domain);
       void init(mesh& domain);
   
       double nLayer;
       double susp_depth;
       double v_edge_height;
   
       // Beta * K, this is beta and scales the eddy diffusivity
       double snow_diffusion_const;
       double l__max;                // vertical mixing length (m)
       bool rouault_diffusion_coeff; // use the spatially variable diffusivity coefficient of Rouault 1991
   
       bool do_fixed_settling;   // should we have a constant settling velocity?
                                 // true: constant settling velocity = settling_velocity (see below)
                                 // false: use the parameterization of Pomeroy et al. (1993) and Pomeroy and Gray (1995):
                                 //        In this case, the settling velocity decreases with height above the snow surface
                                 //        due to a decrease with height in the mean particle size.
       double settling_velocity; // Variable used if do_fixed_settling = true
       double n_non_edge_tri;
       double eps; //lapacian smoothing epilson.
       bool do_sublimation; // should we have a sink sublimation term?
       bool do_lateral_diff; // should have lateral diffusion
       bool enable_veg;      // should we consider vegetation ?
   
       bool use_PomLi_probability; // Use areal Pomeroy Li 2000 probability function.
       bool use_exp_fetch;         // Enable the exp Liston 2006 fetch
       bool use_tanh_fetch;        // Enable the tanh Pomeroy and Male 1986 fetch
   
       bool use_subgrid_topo;    // Enable effect of subgrid topography on snow transport
       bool use_subgrid_topo_V2; // Enable effect of subgrid topography on snow transport
   
       bool iterative_subl; // if True, enables the iterative sublimation calculation as per Pomeroy and Li 2000
       bool use_R94_lambda; // use the Raupach 1990 lambda expression using LAI/2 instead of pomeroy stalk density
   
   
   
       // this is the suspension transport matrix
       double nnz; // number non zero
   
       // this is the drift matrix
       double nnz_drift; // number non zero
   
       bool debug_output;
       double cutoff; // cutoff veg-snow diff (m) that we inhibit saltation entirely
   
       // don't allow transport if below this threshold.
       // This gives models like snobal a chance to build up their snowpack and avoid convergence issues with thin snowcovers
       double min_sd_trans;
   
       // Couple the calculation of u* and z0 via the z0 value from
       // Li and Pomeroy 2000, eqn 5.
       // to modify the u* estimation instead of using a snow z0 for u* estimation
       bool z0_ustar_coupling;
   
       class data : public face_info
       {
         public:
           // edge unit normals
           arma::vec m[5];
   
           // prism areas
           double A[5];
   
           // face neighbors
           bool face_neigh[3];
   
           std::vector<double> u_z_susp; // suspension layer windspeeds
           size_t cell_local_id;
   
           double CanopyHeight;
           double LAI;
   
           // used for the pomeroy stalk formulation
           double N;  // vegetation number density
           double dv; // stalk diameter
   
           // saltation height
           double hs;
   
           bool is_edge;
   
           // used to flag the large vegetation areas or other via landcover types to not do any saltation at this point.
           bool saltation;
   
           double z0;
   
           double sum_drift;
           double sum_subl;
           std::vector<double> csubl; //vertical col of sublimation coeffs
   
           gsl_function F_fill;
           gsl_function F_fill2;
           gsl_function F_fill3;
           // gsl_function F_roots;
           // struct my_fill_topo_params params = { 1.1, 0.3, 0.6 , 0.4};
       };
   
       void checkpoint(mesh& domain,  netcdf& chkpt);
       void load_checkpoint(mesh& domain,  netcdf& chkpt);
   
   private:
   
     // For detecting if there is suspension and/or saltation
     bool suspension_present, deposition_present;
     constexpr static double suspension_present_threshold=1e-12;
     constexpr static double deposition_present_threshold=1e-12;
   
     std::unique_ptr<math::LinearAlgebra::NearestNeighborProblem> deposition_NNP;
     std::unique_ptr<math::LinearAlgebra::NearestNeighborProblem> suspension_NNP;
   
   };