Commit f45307a2 authored by Ronan Trepos's avatar Ronan Trepos
Browse files

Pea: compute coef_fixmax3

- from weather stats (ahead of simulation time)
- Integration of works of M. Chanis on Pea
parent 97b50548
Package: AZODYN
Version: 0.1.0
Depends:
Build-Depends: vle.discrete-time
Build-Depends: vle.discrete-time,vle.reader
Conflicts:
Maintainer: Ronan Trépos <ronan.trepos@inra.fr>
Description: Still in progress
......
......@@ -41,7 +41,7 @@
<port name="TX"/>
</out>
</model>
<model observables="oNitrogenFixationPea" conditions="cBegin,cSoil,cPlant,cNitrogenFixationPea,cPlantPea,cCrop" width="100" dynamics="dNitrogenFixationPea" height="180" x="988" y="1140" name="NitrogenFixationPea" type="atomic">
<model observables="oNitrogenFixationPea" conditions="cBegin,cSoil,cPlant,cNitrogenFixationPea,cPlantPea,cCrop" width="100" dynamics="dNitrogenFixationPea" height="220" x="988" y="1140" name="NitrogenFixationPea" type="atomic">
<in>
<port name="MS"/>
<port name="ProfRacmax"/>
......@@ -52,7 +52,8 @@
<port name="Stade"/>
<port name="dQN_besoin"/>
<port name="date_DRG"/>
</in>
<port name="coef_fixmax3"/>
</in>
<out>
<port name="NJstressH_nod"/>
<port name="Ndfa"/>
......@@ -317,7 +318,13 @@
<port name="I"/>
</out>
</model>
</submodels>
<model observables="" conditions="cBegin,condClimate,cPlant,cPlantPea" width="50" dynamics="dCoef_fixmax3" height="50" y="1290" x="299" name="Coef_fixmax3" type="atomic">
<in/>
<out>
<port name="coef_fixmax3"/>
</out>
</model>
</submodels>
<connections>
<connection type="internal">
<origin port="stressGEL" model="Frost"/>
......@@ -655,6 +662,10 @@
<origin port="MS_MP" model="PlantGrowthPea"/>
<destination port="MS_MP" model="YieldPea"/>
</connection>
<connection type="internal">
<origin model="Coef_fixmax3" port="coef_fixmax3"/>
<destination model="NitrogenFixationPea" port="coef_fixmax3"/>
</connection>
</connections>
</model>
</structures>
......@@ -672,7 +683,8 @@
<dynamic library="StressT" package="AZODYN" name="dStressT"/>
<dynamic library="WaterSoil" package="AZODYN" name="dWaterSoil"/>
<dynamic library="YieldPea" package="AZODYN" name="dYieldPea"/>
</dynamics>
<dynamic library="Coef_fixmax3" package="AZODYN" name="dCoef_fixmax3"/>
</dynamics>
<experiment name="AzodynPea">
<conditions>
<condition name="simulation_engine">
......@@ -701,9 +713,7 @@
<port name="GELmax_seuil">
<double>5</double>
</port>
<port name="GEL_activation">
<double>1.000000000000000</double>
</port>
<port name="date_fin_GEL">
<string>2015-03-31</string>
</port>
......@@ -931,9 +941,7 @@
<port name="coef_fixmax2">
<double>0.011900000000000</double>
</port>
<port name="coef_fixmax3">
<double>0.06888791</double>
</port>
<port name="coef_prot_gr">
<double>6.250000000000000</double>
</port>
......
// @@tagdynamic@@
// @@tagdepends: vle.reader, vle.discrete-time @@endtagdepends
/*
* Copyright (c) 2014-2017 INRA http://www.inra.fr
*
* Permission is hereby granted, free of charge, to any person obtaining a copy
* of this software and associated documentation files (the "Software"), to
* deal in the Software without restriction, including without limitation the
* rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
* sell copies of the Software, and to permit persons to whom the Software is
* furnished to do so, subject to the following conditions:
*
* The above copyright notice and this permission notice shall be included in
* all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
* IN THE SOFTWARE.
*/
/*
* Very close to record.meteo/src/MetoeReader.cpp
*/
#include <iostream>
#include <sstream>
#include <vle/utils/Tools.hpp>
#include <vle/utils/Context.hpp>
#include <vle/utils/Package.hpp>
#include <vle/utils/DateTime.hpp>
#include <vle/DiscreteTime.hpp>
#include <vle/reader/table_file_reader.hpp>
namespace AZODYN {
namespace vv = vle::value;
namespace vu = vle::utils;
using namespace vle::discrete_time;
enum METEO_TYPE {GENERIC_WITH_HEADER, AGROCLIM, DRIAS_ASCII};
class Coef_fixmax3 : public DiscreteTimeDyn
{
public:
Coef_fixmax3(const vle::devs::DynamicsInit& init,
const vle::devs::InitEventList& events)
: DiscreteTimeDyn(init, events), m_table_reader(),
meteo_type(GENERIC_WITH_HEADER), line_read(), vars(), year_i(-1),
month_i(-1),day_i(-1), nb_compute(0),
begin_date(std::numeric_limits<double>::infinity())
{
int indexVar_TMC = 9999999;
if (!events.exist("meteo_type")) {
throw vu::ArgError(vu::format("[%s] meteo_type parameter is missing!" ,
getModelName().c_str()));
} else {
const std::string& type = events.getString("meteo_type");
if (type == "generic_with_header") {
meteo_type = GENERIC_WITH_HEADER;
} else if (type == "agroclim") {
meteo_type = AGROCLIM;
} else if (type == "drias_ascii") {
meteo_type = DRIAS_ASCII;
} else {
throw vu::ArgError(vu::format("[%s] Unexpected meteo_type `%s`." ,
getModelName().c_str(), type.c_str()));
}
}
if (events.exist("column_separator")) {
m_table_reader.getParams().separator =
events.getString("column_separator");
} else {
switch (meteo_type) {
case GENERIC_WITH_HEADER: {
m_table_reader.getParams().separator = "\t";
break;
} case AGROCLIM: {
m_table_reader.getParams().separator = ";";
break;
} case DRIAS_ASCII: {
m_table_reader.getParams().separator = ";";
break;
} default:
throw vu::ArgError(vu::format("[%s]: column_separator parameter have no default value defined for meteo_type : `%s`.",
getModelName().c_str(), enum2str(meteo_type)));
}
}
std::string meteo_file_path;
if (not events.exist("meteo_file")) {
throw vu::ArgError(vu::format("[%s]: meteo_file parameter is missing!",
getModelName().c_str()));
}
if (events.exist("PkgName") and events.getString("PkgName") != "") {
vu::Package pkg(context(), events.getString("PkgName"));
meteo_file_path = pkg.getDataFile(events.getString("meteo_file"));
} else {
meteo_file_path = events.getString("meteo_file");
}
m_table_reader.setFilePath(meteo_file_path);
//configure table file reader
switch (meteo_type) {
case GENERIC_WITH_HEADER: {
m_table_reader.readLine(line_read);
vv::Set::iterator itb = line_read.begin();
vv::Set::iterator ite = line_read.end();
for (; itb!=ite; itb++) {
vars.push_back(Var());
vars[vars.size()-1].init(this, (*itb)->toString().value(), events);
}
for (unsigned int i=0; i<vars.size(); i++) {
m_table_reader.getParams().col_types.push_back(vv::Value::DOUBLE);
}
break;
} case AGROCLIM: {
std::string line;
bool skip = true;
while (skip) {
m_table_reader.readLine(line);
if (line.substr(0, 9) == "NUM_POSTE") {
std::vector<std::string> strs;
vle::utils::tokenize(line, strs, ";", false);
for (unsigned int i=0; i<strs.size();i++) {
vars.push_back(Var());
vars[vars.size()-1].init(this, strs[i], events);
m_table_reader.getParams().col_types.push_back(
vv::Value::DOUBLE);
}
skip = false;
}
}
break;
} case DRIAS_ASCII: {
std::string line;
bool skip = true;
while (skip) {
m_table_reader.readLine(line);
if (line.substr(0, 22) == "# Format de la ligne :") {
skip = false;
}
}
m_table_reader.readLine(line);
line = line.substr(11, line.size() - 12);
line = "AN MOIS JOUR " + line;
std::vector<std::string> strs;
vle::utils::tokenize(line, strs, " ", false);
for (unsigned int i=0; i<strs.size();i++) {
vars.push_back(Var());
vars[vars.size()-1].init(this, strs[i], events);
m_table_reader.getParams().col_types.push_back(
vv::Value::DOUBLE);
}
skip = true;
while (skip) {
m_table_reader.readLine(line);
if (line.substr(0, 57) == "# http://www.drias-climat.fr/accompagnement/conditionsEDC") {
skip = false;
m_table_reader.readLine(line);
}
}
break;
} default:
throw vu::ArgError(vu::format("[%s] Unexpected meteo_type `%s`" ,
getModelName().c_str(), enum2str(meteo_type)));
}
//jump to the correct line
if (events.exist("begin_date")) {
if (events.get("begin_date")->isDouble()) {
begin_date = events.getDouble("begin_date");
} else if (events.get("begin_date")->isString()) {
begin_date = vu::DateTime::toJulianDay(
events.getString("begin_date"));
}
std::string year_col ="AN";//default agroclim
std::string month_col ="MOIS";//default agroclim
std::string day_col = "JOUR";//default agroclim
if (events.exist("year_column")) {
year_col.assign(events.getString("year_column"));
}
if (events.exist("month_column")) {
month_col.assign(events.getString("month_column"));
}
if (events.exist("day_column")) {
day_col.assign(events.getString("day_column"));
}
std::vector<Var>::const_iterator itb = vars.begin();
std::vector<Var>::const_iterator ite = vars.end();
int i =0;
for (; itb != ite; itb++) {
if (itb->name == year_col) {
year_i = i;
} else if (itb->name == month_col) {
month_i = i;
} else if (itb->name == day_col) {
day_i = i;
}
if (itb->name == "TMC") {
indexVar_TMC = i;
}
i++;
}
if (year_i == -1) {
throw vu::ArgError(vu::format("[%s] Index of column `%s` not found for year_column",
getModelName().c_str(), year_col.c_str()));
}
if (month_i == -1) {
throw vu::ArgError(vu::format("[%s] Index of column `%s` not found for month_column",
getModelName().c_str(), month_col.c_str()));
}
if (day_i == -1) {
throw vu::ArgError(vu::format("[%s] Index of column `%s` not found for day_column",
getModelName().c_str(), day_col.c_str()));
}
//read first line to get beginning date of file
vle::value::Set lineToFill;
res = m_table_reader.readLine(line_read);
if (line_read.empty()) {
throw vu::ArgError(vu::format("[%s] No data found in file '%s'",
getModelName().c_str(),
m_table_reader.getFilePath().c_str()));
}
std::stringstream ss;
ss << line_read.getDouble(year_i) << "-" << line_read.getDouble(month_i) << "-" << line_read.getDouble(day_i);
int begin_date_file = vu::DateTime::toJulianDayNumber(ss.str());
if (begin_date_file > int(begin_date + 1)) {
throw vu::ArgError(vu::format("[%s] Cannot find data in file '%s' for begin_date equals: `%s` ",
getModelName().c_str(),
m_table_reader.getFilePath().c_str(),
vu::DateTime::toJulianDay(begin_date).c_str()));
}
for (int toskip = int(begin_date + 1) - begin_date_file;
toskip > 0; toskip--) {
res = m_table_reader.readLine(line_read);
}
} else {
res = m_table_reader.readLine(line_read);
}
//*** compute coef_fixmax3
//jump to sowing
coef_fixmax3.init(this, "coef_fixmax3", events);
int current_time = int(begin_date);
int sowing_time = vu::DateTime::toJulianDayNumber(events.getString("date_S"));
for (; current_time < sowing_time; current_time++) {}
//sum temperature betwenn sowing and crop emergence
double sumT_sowing_emergence = 0;
double ST_S_LEVEE = events.getDouble("ST_S_LEVEE");
for (; sumT_sowing_emergence < ST_S_LEVEE ; current_time++) {
res = m_table_reader.readLine(line_read);
sumT_sowing_emergence += line_read.getDouble(indexVar_TMC);
}
std::string date_levee = vu::DateTime::toJulianDay(current_time);
//sum temperature between crop emergence and flowering
double sumT_emergence_flowering = 0;
int flowering_time = vu::DateTime::toJulianDayNumber(events.getString("date_FLO"));
for (; current_time <= flowering_time ; current_time++) {
res = m_table_reader.readLine(line_read);
sumT_emergence_flowering += line_read.getDouble(indexVar_TMC);
}
//regression equation
double coef_fixmax1 = events.getDouble("coef_fixmax1");
double coef_fixmax2 = events.getDouble("coef_fixmax2");
double ST_DF_DRG = events.getDouble("ST_DF_DRG");
coef_fixmax3.init_value(
coef_fixmax2 - (coef_fixmax1*(sumT_emergence_flowering + ST_DF_DRG)));
}
virtual ~Coef_fixmax3()
{
}
void compute(const vle::devs::Time& /*t*/) override
{
}
vle::reader::TableFileReader m_table_reader;
METEO_TYPE meteo_type;
vv::Set line_read;
std::vector<Var> vars;
Var coef_fixmax3;
int year_i;
int month_i;
int day_i;
double nb_compute;
double begin_date;
bool res;
private:
char* enum2str(METEO_TYPE u) {
switch (u) {
case GENERIC_WITH_HEADER:
return (char*)&"generic_with_header";
case AGROCLIM:
return (char*)&"agroclim";
case DRIAS_ASCII:
return (char*)&"drias_ascii";
default:
return (char*)&"INVALID ENUM";
}
}
};
}//namespaces
DECLARE_DYNAMICS(AZODYN::Coef_fixmax3)
......@@ -38,6 +38,7 @@ public:
STapLEVEE.init(this, "STapLEVEE", events);
date_DRG.init(this, "date_DRG", events);
MS.init(this, "MS", events);
coef_fixmax3.init(this, "coef_fixmax3", events);
Profnod.init(this, "Profnod", events);
RU_nod.init(this, "RU_nod", events);
......@@ -116,7 +117,7 @@ public:
dQNfixmax_tot = 0.0 ;
} else if (Stade() < DRG_FSLA){
dQNfixmax_tot = std::max(0.0,pp.coef_tot_a * (MS()-MS(-1))*10.0*
(pp.coef_fixmax1*STapLEVEE() + pp.coef_fixmax3));
(pp.coef_fixmax1*STapLEVEE() + coef_fixmax3()));
} else if (Stade() < FSLA_MP){
dQNfixmax_tot =std::max(0.0,(MS()-MS(-1))*10.0*pp.coef_fixmax2
*pp.coef_tot_a);
......@@ -166,6 +167,7 @@ private:
/*Sync*/ Var STapLEVEE ;
/*Sync*/ Var date_DRG;
/*Sync*/ Var MS;
/*Sync*/ Var coef_fixmax3;;
Var Profnod;
Var RU_nod ;
......
......@@ -78,7 +78,6 @@ public:
//sert pour calculer QNfixmax
double coef_fixmax1;
double coef_fixmax2;
double coef_fixmax3;
/*** Plant growth **/
......@@ -167,7 +166,6 @@ public:
STfixmax_seuil = Utils::extractDouble(events, "STfixmax_seuil");
coef_fixmax1 = Utils::extractDouble(events, "coef_fixmax1");
coef_fixmax2 = Utils::extractDouble(events, "coef_fixmax2");
coef_fixmax3 = Utils::extractDouble(events, "coef_fixmax3");
Ebv = Utils::extractDouble(events, "Ebv");
coef_MS_init = Utils::extractDouble(events, "coef_MS_init");
vRac_RG = Utils::extractDouble(events, "vRac_RG");
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment