Subversion Repositories public iLand

/********************************************************************************************

**    iLand - an individual based forest landscape and disturbance model

**    http://iland.boku.ac.at

**    Copyright (C) 2009-  Werner Rammer, Rupert Seidl

**

**    This program 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.

**

**    This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.

********************************************************************************************/

#ifndef WINDMODULE_H

#define WINDMODULE_H

#include "grid.h"

#include "layeredgrid.h"

#include "expression.h"

#include <QHash>

class Tree; // forward

class Species; // forward

class ResourceUnit; // forward

/** data structure for a single wind cell (usually 10x10m).

  @ingroup windmodule

*/

class WindCell {

public:

    WindCell() { topex = 0; n_affected=0; sum_volume_killed=0.; edge_age=0; clear(); }

    void clear() {height = edge = 0.f; n_trees=0.f; tree=0; n_killed = 0; basal_area_killed = 0.f; cws_uproot = 0.; cws_break= crown_windspeed= 0.; n_iteration = 0;}

    bool isValid() const { return height<9999.f; } ///< returns true if the pixel is on the valid project area

    float topex; ///< topographic modifier for wind speed (-)

    float height; ///< top height (m).

    float n_trees; ///< number of trees on pixel

    const Tree *tree; ///< pointer to the tallest tree on the pixel (if already populated)

    float edge; ///< maximum difference to neighboring cells (m)

    // statistics

    int n_iteration; ///< number of iteration this pixel is processed (and trees are killed)

    int n_killed; ///< number of trees killed on the pixel

    int edge_age; ///< age of an edge (consecutive number of years of being an edge)

    double basal_area_killed; ///< basal area of trees that died (m2) (during an event)

    double cws_uproot; ///< critital wind speed for uprooting (m/s)

    double cws_break; ///< critical wind speed for tree breakage (m/s)

    double crown_windspeed; ///< wind speed (m/s) on the cecll

    int n_affected; ///< number of storm that killed trees on that pixel.

    double sum_volume_killed; ///< running sum of killed tree volume on that pixel (m3)

};

// data structure for a resource unit

class WindRUCell {

public:

    WindRUCell(): flag(0), soilIsFrozen(false) {}

    int flag; // flag to indicate that the trees of the current resource are already processed

    bool soilIsFrozen;

};

/** Helper class manage and visualize data layers related to fire.

  @ingroup windmodule

*/

class WindLayers: public LayeredGrid<WindCell> {

  public:

    void setGrid(const Grid<WindCell> &grid) { mGrid = &grid; }

    double value(const WindCell& data, const int index) const;

    const QVector<LayerElement> &names();

    // specifics for wind layers

    void setRUGrid(const Grid<WindRUCell> *grid) { mRUGrid = grid; }

private:

    QVector<LayerElement> mNames;

    double ruValueAt(const WindCell *cell, const int what) const; // get resource-unit level factor at cell "cell"

    const Grid<WindRUCell> *mRUGrid;

};

/** Species parameters that are specific to the wind module

  */

struct WindSpeciesParameters

{

    WindSpeciesParameters():  crown_area_factor(0.5), crown_length(0.5), Creg(111), MOR(30.6), wet_biomass_factor(1.86) {}

    double crown_area_factor; // empirical factor related to the crown shape (fraction of crown shape compared to rectangle)

    double crown_length; // crown length of the tree (fraction of tree height)

    double Creg; // Nm/kg, critical turning coefficient from tree pulling

    double MOR; // MPa, modulus of rupture

    double wet_biomass_factor; // conversion factor between dry and wet biomass (wet = dry*factor)

};

/** @class WindModule

    @ingroup windmodule

    The WindModule is the disturbance module for simulation wind and windthrow in iLand.  */

class WindModule

{

public:

    WindModule();

    static double cellsize() { return 10.; }

    /// the general setup routine after starting iland

    void setup();

    void setupResourceUnit(const ResourceUnit* ru);

    WindLayers &layers() { return mWindLayers; }

    /// main function of the disturbance module

    void run(const int iteration=-1, const bool execute_from_script=false);

    // test functions

    void setWindProperties(const double direction_rad, const double speed_ms) { mWindDirection = direction_rad; mWindSpeed = speed_ms; }

    void setSimulationMode(const bool mode) { mSimulationMode = mode; }

    bool simulationMode() const { return mSimulationMode; }

    void setMaximumIterations(const double maxit) { mMaxIteration = maxit; }

    void testFetch(double degree_direction);

    void testEffect();

private:

    // main functions

    bool eventTriggered(); ///< determine details of this years' wind event (and return false if no event happens)

    void initWindGrid(); ///< load state from iland main module

    void detectEdges(bool at_startup=false); ///< detect all pixels that are higher than the surrounding and therefore are likely candidates for damage

    void calculateFetch(); ///< calculate maximum gap sizes in upwind direction

    int calculateWindImpact(); ///< do one round of wind effect calculations

    // details

    /// find distance to the next pixels that give shelter

    bool checkFetch(const int startx, const int starty, const double direction, const double max_distance, const double threshold) ;

    /// perform the wind effect calculations for a given grid cell

    bool windImpactOnPixel(const QPoint position, WindCell *cell);

    ///

    double calculateCrownWindSpeed(const Tree *tree, const WindSpeciesParameters &params, const double n_trees, const double wind_speed_10);

    double calculateCrititalWindSpeed(const Tree *tree, const WindSpeciesParameters &params, const double gap_length, double &rCWS_uproot, double &rCWS_break);

    /// check if the soil on current resource unit is frozen or not

    bool isSoilFrozen(const ResourceUnit *ru, const int day_of_year) const;

    // helping functions

    void scanResourceUnitTrees(const QPoint &position);

    void loadSpeciesParameter(const QString &table_name);

    const WindSpeciesParameters &speciesParameter(const Species *s);

    void initializeEdgeAge(const int years);

    void incrementEdgeAge();

    // variables

    bool mTopexFromGrid; ///< indicates if topex grid is a real grid or by ressource unit

    double mWindDirection; ///< direction of the current wind event (rad)

    double mWindDirectionVariation; ///< random variation in wind direction

    double mWindSpeed; ///< wind speed (TODO: per resource unit!)

    double mEdgeDetectionThreshold; ///< minimum height difference of height-grid pixels to be detected as edges (default is 10m)

    double mFactorEdge; ///< constant ratio between the maximujm turning moments at the stand edge and conditions well inside the forest (default: 5)

    int mWindDayOfYear; ///< day of year of the wind event (0..365)

    bool mSimulationMode; ///< if true, no trees are removed (test mode)

    int mCurrentIteration; ///< current iteration (1..n)

    int mMaxIteration; ///< maximum number of iterations

    double mGustModifier; ///< Modification range accounting for differences in wind speed between iterations (0..1)

    double mCurrentGustFactor; ///< gustfactor of the current year (multiplier)

    enum ETopexFactorModificationType {gfMultiply, gfAdd} mTopexFactorModificationType; ///< determines if topo-modifier is added multiplicatively or additively.

    double mIterationsPerMinute; ///< number of iterations per minute of the events' duration

    int mEdgeAgeBaseValue; ///< initial value for the age of edges in the initial state of the forest

    Expression mEdgeProbability; ///< function that determines the probability that an edge is calculated (old edges are less likely to being tested)

    double mEdgeBackgroundProbability; ///< probability that a 10x10m pixel is marked as edge during a wind event (randomly)

    // some statistics

    int mPixelAffected; ///< total number of pixels that are impacted

    int mTreesKilled; ///< total number of killed trees

    double mTotalKilledBasalArea; ///< total basal area of killed trees

    double mTotalKilledVolume; ///< total tree volume (m3) killed

    enum ESoilFreezeMode {esfFrozen, esfNotFrozen, esfAuto, esfInvalid} mSoilFreezeMode; ///< if "esfAuto", soil-freeze-state is derived from climate

    Grid<WindCell> mGrid; ///< wind grid (10x10m)

    Grid<WindRUCell> mRUGrid; ///< grid for resource unit data

    WindLayers mWindLayers; ///< helping structure

    // species parameters for the wind module

    QHash<const Species*, WindSpeciesParameters> mSpeciesParameters;

    /// formula for the transfer function LRI

    Expression mLRITransferFunction;

    QString mAfterExecEvent;

    friend class WindScript;

    friend class WindOut;

    friend void nc_calculateFetch(WindCell *begin, WindCell *end);

    friend void nc_calculateWindImpact(ResourceUnit *unit);

};

#endif // WINDMODULE_H