Program Listing for File PBSM3D.hpp

↰ Return to documentation for file (src/modules/PBSM3D.hpp)

//
// 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;

};