Skip to content

Commit

Permalink
version
Browse files Browse the repository at this point in the history
  • Loading branch information
bergm committed May 29, 2024
1 parent 6005c2d commit 11c0bb2
Show file tree
Hide file tree
Showing 10 changed files with 71 additions and 102 deletions.
3 changes: 3 additions & 0 deletions .idea/misc.xml

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

95 changes: 35 additions & 60 deletions src/core/crop-module.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,7 @@ Copyright (C) Leibniz Centre for Agricultural Landscape Research (ZALF)
#include "voc-common.h"
#include "photosynthesis-FvCB.h"
#include "O3-impact.h"
#include "io/output.h"

const double PI = 3.14159265358979323;

Expand Down Expand Up @@ -492,7 +493,6 @@ void CropModule::deserialize(mas::schema::model::monica::CropModuleState::Reader
dyingOut = reader.getDyingOut();
vc_AccumulatedETa = reader.getAccumulatedETa();
vc_AccumulatedTranspiration = reader.getAccumulatedTranspiration();
vc_AccumulatedPrimaryCropYield = reader.getAccumulatedPrimaryCropYield();
vc_sumExportedCutBiomass = reader.getSumExportedCutBiomass();
vc_exportedCutBiomass = reader.getExportedCutBiomass();
vc_sumResidueCutBiomass = reader.getSumResidueCutBiomass();
Expand Down Expand Up @@ -758,7 +758,6 @@ void CropModule::serialize(mas::schema::model::monica::CropModuleState::Builder
builder.setDyingOut(dyingOut);
builder.setAccumulatedETa(vc_AccumulatedETa);
builder.setAccumulatedTranspiration(vc_AccumulatedTranspiration);
builder.setAccumulatedPrimaryCropYield(vc_AccumulatedPrimaryCropYield);
builder.setSumExportedCutBiomass(vc_sumExportedCutBiomass);
builder.setExportedCutBiomass(vc_exportedCutBiomass);
builder.setSumResidueCutBiomass(vc_sumResidueCutBiomass);
Expand Down Expand Up @@ -1176,32 +1175,27 @@ pair<double, double> CropModule::fc_VernalisationFactor(double vw_MeanAirTempera
* @author Claas Nendel
*/
double CropModule::fc_OxygenDeficiency(double d_CriticalOxygenContent) {
double vc_AirFilledPoreVolume = 0.0;
double vc_MaxOxygenDeficit = 0.0;

int timeUnderAnoxiaThresholdAtStage = vc_DevelopmentalStage < vc_TimeUnderAnoxiaThreshold.size()
? vc_TimeUnderAnoxiaThreshold.at(vc_DevelopmentalStage) : TimeUnderAnoxiaThresholdDefault;
// Reduktion bei Luftmangel Stauwasser berücksichtigen!!!!
vc_AirFilledPoreVolume =
((soilColumn[0].vs_Saturation() + soilColumn[1].vs_Saturation() + soilColumn[2].vs_Saturation()) -
(soilColumn[0].get_Vs_SoilMoisture_m3() + soilColumn[1].get_Vs_SoilMoisture_m3() +
soilColumn[2].get_Vs_SoilMoisture_m3())) / 3.0;
if (vc_AirFilledPoreVolume < d_CriticalOxygenContent) {
vc_TimeUnderAnoxia += int(vc_TimeStep);
if (vc_TimeUnderAnoxia > vc_TimeUnderAnoxiaThreshold) {
vc_TimeUnderAnoxia = vc_TimeUnderAnoxiaThreshold;
}
if (vc_AirFilledPoreVolume < 0.0) {
vc_AirFilledPoreVolume = 0.0;
}
vc_MaxOxygenDeficit = vc_AirFilledPoreVolume / d_CriticalOxygenContent;
vc_OxygenDeficit = 1.0 - double(vc_TimeUnderAnoxia / vc_TimeUnderAnoxiaThreshold) * (1.0 - vc_MaxOxygenDeficit);
double sumSaturation = 0, sumSoilMoisture = 0;
int sumLayers = 0;
auto nols = std::min(std::max(3UL, vc_RootingDepth), soilColumn.vs_NumberOfLayers());
for(size_t i = 0; i < nols; i++) {
sumSaturation += soilColumn[i].vs_Saturation();
sumSoilMoisture += soilColumn[i].get_Vs_SoilMoisture_m3();
sumLayers++;
}
double avgAirFilledPoreVolume = (sumSaturation - sumSoilMoisture) / sumLayers;
if (avgAirFilledPoreVolume < d_CriticalOxygenContent) {
vc_TimeUnderAnoxia = std::min(vc_TimeUnderAnoxia + int(vc_TimeStep), timeUnderAnoxiaThresholdAtStage);
avgAirFilledPoreVolume = std::max(0.0, avgAirFilledPoreVolume);
double vc_MaxOxygenDeficit = avgAirFilledPoreVolume / d_CriticalOxygenContent;
vc_OxygenDeficit = 1.0 - double(vc_TimeUnderAnoxia / double(timeUnderAnoxiaThresholdAtStage)) * (1.0 - vc_MaxOxygenDeficit);
} else {
vc_TimeUnderAnoxia = 0;
vc_OxygenDeficit = 1.0;
}
if (vc_OxygenDeficit > 1.0) {
vc_OxygenDeficit = 1.0;
}

return vc_OxygenDeficit;
}

Expand Down Expand Up @@ -2242,7 +2236,7 @@ void CropModule::fc_CropPhotosynthesis(double vw_MeanAirTemperature,
Voc::SpeciesData species;
// species.id = 0; // right now we just have one crop at a time, so no need to distinguish multiple crops
species.lai = LAI;
species.mFol = get_OrganGreenBiomass(LEAF) / (100. * 100.); // kg/ha -> kg/m2
species.mFol = get_OrganGreenBiomass(OId::LEAF) / (100. * 100.); // kg/ha -> kg/m2
species.sla =
species.mFol > 0 ? species.lai / species.mFol : pc_SpecificLeafArea[vc_DevelopmentalStage] * 100. *
100.; // ha/kg -> m2/kg
Expand Down Expand Up @@ -2290,7 +2284,7 @@ void CropModule::fc_CropPhotosynthesis(double vw_MeanAirTemperature,
// JJV
for (const auto &lf: {FvCB_res.sunlit, FvCB_res.shaded}) {
species.lai = lf.LAI;
species.mFol = get_OrganGreenBiomass(LEAF) / (100. * 100.) * lf.LAI /
species.mFol = get_OrganGreenBiomass(OId::LEAF) / (100. * 100.) * lf.LAI /
(sun_LAI + sh_LAI); // kg/ha -> kg/m2
species.sla =
species.mFol > 0 ? species.lai / species.mFol : pc_SpecificLeafArea[vc_DevelopmentalStage] * 100. *
Expand Down Expand Up @@ -2885,7 +2879,7 @@ void CropModule::fc_CropDryMatter(double vw_MeanAirTemperature) {
//! TODO: hard coded organ ids; must be more generalized because in database organ_ids can be mixed
// vc_OrganBiomass[i_Organ];

if (i_Organ == LEAF) { // leaf
if (i_Organ == OId::LEAF) { // leaf

double incr = assimilate_partition_leaf * vc_NetPhotosynthesis;
if (fabs(incr) <= vc_OrganBiomass[i_Organ]) {
Expand All @@ -2908,7 +2902,7 @@ void CropModule::fc_CropDryMatter(double vw_MeanAirTemperature) {
// vc_OrganGrowthIncrement[i_Organ] = incr;
// debug() << "LEAF - Reducing organ by " << incr << " (" << vc_OrganBiomass[i_Organ] + vc_OrganGrowthIncrement[i_Organ] << ")"<< endl;
}
} else if (i_Organ == SHOOT) { // shoot
} else if (i_Organ == OId::SHOOT) { // shoot

double incr = assimilate_partition_leaf * vc_NetPhotosynthesis; // should be negative

Expand Down Expand Up @@ -2996,10 +2990,10 @@ void CropModule::fc_CropDryMatter(double vw_MeanAirTemperature) {

// update the root biomass and dead root biomass vars
// root dead biomass will be transfered to proper AOM pools
if (i_Organ == ROOT) {
vc_OrganBiomass[ROOT] -= dailyDeadBiomassIncrement[ROOT];
vc_OrganDeadBiomass[ROOT] -= dailyDeadBiomassIncrement[ROOT];
vc_TotalBiomassNContent -= dailyDeadBiomassIncrement[ROOT] * vc_NConcentrationRoot;
if (i_Organ == OId::ROOT) {
vc_OrganBiomass[OId::ROOT] -= dailyDeadBiomassIncrement[OId::ROOT];
vc_OrganDeadBiomass[OId::ROOT] -= dailyDeadBiomassIncrement[OId::ROOT];
vc_TotalBiomassNContent -= dailyDeadBiomassIncrement[OId::ROOT] * vc_NConcentrationRoot;
}
} else {
vc_OrganBiomass[i_Organ] += (vc_OrganGrowthIncrement[i_Organ] * vc_TimeStep); // [kg CH2O ha-1]
Expand Down Expand Up @@ -3100,10 +3094,11 @@ void CropModule::fc_CropDryMatter(double vw_MeanAirTemperature) {

// Determining root penetration rate according to soil clay content [m °C-1 d-1]
double vc_RootPenetrationRate = 0.0; // [m °C-1 d-1]
if (soilColumn[vc_RootingDepth].vs_SoilClayContent() <= 0.02) {
auto layerIndexBelowRootingDepth = std::min(vc_RootingDepth, nols - 1);
if (soilColumn[layerIndexBelowRootingDepth].vs_SoilClayContent() <= 0.02) {
vc_RootPenetrationRate = 0.5 * pc_RootPenetrationRate;
} else if (soilColumn[vc_RootingDepth].vs_SoilClayContent() <= 0.08) {
vc_RootPenetrationRate = ((1.0 / 3.0) + (0.5 / 0.06 * soilColumn[vc_RootingDepth].vs_SoilClayContent())) *
} else if (soilColumn[layerIndexBelowRootingDepth].vs_SoilClayContent() <= 0.08) {
vc_RootPenetrationRate = ((1.0 / 3.0) + (0.5 / 0.06 * soilColumn[layerIndexBelowRootingDepth].vs_SoilClayContent())) *
pc_RootPenetrationRate; // [m °C-1 d-1]
} else {
vc_RootPenetrationRate = pc_RootPenetrationRate; // [m °C-1 d-1]
Expand All @@ -3121,28 +3116,12 @@ void CropModule::fc_CropDryMatter(double vw_MeanAirTemperature) {
vc_RootingDepth_m += (vc_DailyTemperatureRoot * vc_RootPenetrationRate); // [m]
}

if (vc_RootingDepth_m <= pc_InitialRootingDepth) {
vc_RootingDepth_m = pc_InitialRootingDepth;
}

if (vc_RootingDepth_m > vc_MaxRootingDepth) {
vc_RootingDepth_m = vc_MaxRootingDepth;
} // [m]

if (vc_RootingDepth_m > vs_MaxEffectiveRootingDepth) {
vc_RootingDepth_m = vs_MaxEffectiveRootingDepth;
}

// Calculating rooting depth layer []
vc_RootingDepth = int(std::floor(0.5 + (vc_RootingDepth_m / layerThickness))); // []
if (vc_RootingDepth > nols) {
vc_RootingDepth = nols;
}
if (vc_RootingDepth_m <= pc_InitialRootingDepth) vc_RootingDepth_m = pc_InitialRootingDepth;
if (vc_RootingDepth_m > vc_MaxRootingDepth) vc_RootingDepth_m = vc_MaxRootingDepth; // [m]
if (vc_RootingDepth_m > vs_MaxEffectiveRootingDepth) vc_RootingDepth_m = vs_MaxEffectiveRootingDepth;

vc_RootingZone = int(std::floor(0.5 + ((1.3 * vc_RootingDepth_m) / layerThickness))); // []
if (vc_RootingZone > nols) {
vc_RootingZone = nols;
}
vc_RootingDepth = std::min(int(std::round(vc_RootingDepth_m/layerThickness)), int(nols)); // layer no
vc_RootingZone = std::min(int(std::round(1.3 * vc_RootingDepth_m / layerThickness)), int(nols)); // layer no

vc_TotalRootLength = vc_RootBiomass * pc_SpecificRootLength; //[m m-2]

Expand Down Expand Up @@ -3949,7 +3928,7 @@ void CropModule::calculateVOCEmissions(const Voc::MicroClimateData &mcd) {
Voc::SpeciesData species;
// species.id = 0; // right now we just have one crop at a time, so no need to distinguish multiple crops
species.lai = get_LeafAreaIndex();
species.mFol = get_OrganBiomass(LEAF) / (100. * 100.); // kg/ha -> kg/m2
species.mFol = get_OrganBiomass(OId::LEAF) / (100. * 100.); // kg/ha -> kg/m2
species.sla = pc_SpecificLeafArea[vc_DevelopmentalStage] * 100. * 100.; // ha/kg -> m2/kg

species.EF_MONO = speciesPs.EF_MONO;
Expand Down Expand Up @@ -4655,10 +4634,6 @@ double CropModule::get_AccumulatedTranspiration() const {
return vc_AccumulatedTranspiration;
}

double CropModule::get_AccumulatedPrimaryCropYield() const {
return vc_AccumulatedPrimaryCropYield;
}

/**
* Returns the depth of the maximum active and effective root.
* [m]
Expand Down
11 changes: 2 additions & 9 deletions src/core/crop-module.h
Original file line number Diff line number Diff line change
Expand Up @@ -40,13 +40,6 @@ Copyright (C) Leibniz Centre for Agricultural Landscape Research (ZALF)

namespace monica {

enum {
ROOT = 0,
LEAF = 1,
SHOOT = 2,
STORAGE_ORGAN = 3
};

/*
* @brief Crop part of model
*
Expand Down Expand Up @@ -656,7 +649,8 @@ class CropModule {
double vc_TargetNConcentration{0.0}; //! old GEHMAX
double vc_TimeStep{1.0}; //! old dt
int vc_TimeUnderAnoxia{0};
int vc_TimeUnderAnoxiaThreshold{4};
int TimeUnderAnoxiaThresholdDefault = 4;
std::vector<int> vc_TimeUnderAnoxiaThreshold;
double vs_Tortuosity{}; //! old AD
double vc_TotalBiomass{0.0};
double vc_TotalBiomassNContent{0.0}; //! old PESUM
Expand All @@ -680,7 +674,6 @@ class CropModule {
bool dyingOut{false};
double vc_AccumulatedETa{0.0};
double vc_AccumulatedTranspiration{0.0};
double vc_AccumulatedPrimaryCropYield{0.0};
double vc_sumExportedCutBiomass{0.0};
double vc_exportedCutBiomass{0.0};
double vc_sumResidueCutBiomass{0.0};
Expand Down
17 changes: 8 additions & 9 deletions src/core/monica-model.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -321,7 +321,7 @@ void MonicaModel::seedCrop(Crop *crop) {
* Deletes the current crop.
*/
void
MonicaModel::harvestCurrentCrop(bool exported, Harvest::Spec spec, Harvest::OptCarbonManagementData optCarbMgmtData) {
MonicaModel::harvestCurrentCrop(bool exported, const Harvest::Spec& spec, Harvest::OptCarbonManagementData optCarbMgmtData) {
if (_currentCropModule) {
// prepare to add root and crop residues to soilorganic (AOMs)
// dead root biomass has already been added daily, so just living root biomass is left
Expand Down Expand Up @@ -404,22 +404,21 @@ MonicaModel::harvestCurrentCrop(bool exported, Harvest::Spec spec, Harvest::OptC
auto organIdsForPrimaryYield = _currentCropModule->organIdsForPrimaryYield();
for (const auto &p: spec.organ2specVal) {
// ignore root, is probably an error, when the user specified the root organ (0) as something to harvest
if (p.first == 0)
continue;
if (p.first == 0) continue;
auto organBiomass = _currentCropModule->get_OrganBiomass(p.first);
auto organYield = organBiomass * p.second.exportPercentage / 100.0;
cropYield += organYield;
if (organIdsForPrimaryYield.find(p.first + 1) != organIdsForPrimaryYield.end())
if (organIdsForPrimaryYield.find(p.first + 1) != organIdsForPrimaryYield.end()) {
primaryCropYield += organYield;
if (p.second.incorporate)
sumOrganResidueBiomassToIncorporate += organBiomass - organYield;
else
sumOrganResidueBiomassAsOverlay += organBiomass - organYield;
}
if (p.second.incorporate) sumOrganResidueBiomassToIncorporate += organBiomass - organYield;
else sumOrganResidueBiomassAsOverlay += organBiomass - organYield;
}
auto totalResidueBiomass = _currentCropModule->get_ResidueBiomass(false, cropYield);
auto totalResidueBiomassToIncorporate = totalResidueBiomass - sumOrganResidueBiomassAsOverlay;
auto residuesNConcentration = _currentCropModule->get_ResiduesNConcentration(primaryCropYield);
_soilOrganic->addOrganicMatter(_currentCropModule->residueParameters(), totalResidueBiomassToIncorporate,
_soilOrganic->addOrganicMatter(_currentCropModule->residueParameters(),
totalResidueBiomassToIncorporate,
residuesNConcentration);

debug()
Expand Down
2 changes: 1 addition & 1 deletion src/core/monica-model.h
Original file line number Diff line number Diff line change
Expand Up @@ -81,7 +81,7 @@ class MonicaModel {

bool isCropPlanted() const { return _currentCropModule; }

void harvestCurrentCrop(bool exported, Harvest::Spec spec, Harvest::OptCarbonManagementData optCarbMgmtData = Harvest::OptCarbonManagementData());
void harvestCurrentCrop(bool exported, const Harvest::Spec& spec, Harvest::OptCarbonManagementData optCarbMgmtData = Harvest::OptCarbonManagementData());

void incorporateCurrentCrop();

Expand Down
7 changes: 6 additions & 1 deletion src/core/monica-parameters.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -1456,7 +1456,12 @@ Errors CropModuleParameters::merge(json11::Json j) {
set_double_value(pc_GrowthRespirationParameter2, j, "GrowthRespirationParameter2");
set_double_value(pc_Tortuosity, j, "Tortuosity");
set_bool_value(pc_AdjustRootDepthForSoilProps, j, "AdjustRootDepthForSoilProps");
set_int_value(pc_TimeUnderAnoxiaThreshold, j, "TimeUnderAnoxiaThreshold");
if (j["TimeUnderAnoxiaThreshold"].is_number()) {
std::fill(pc_TimeUnderAnoxiaThreshold.begin(), pc_TimeUnderAnoxiaThreshold.end(),
j["TimeUnderAnoxiaThreshold"].number_value());
} else if(j["TimeUnderAnoxiaThreshold"].is_array()) {
set_int_vector(pc_TimeUnderAnoxiaThreshold, j, "TimeUnderAnoxiaThreshold");
}

set_bool_value(__enable_Photosynthesis_WangEngelTemperatureResponse__, j,
"__enable_Photosynthesis_WangEngelTemperatureResponse__");
Expand Down
2 changes: 1 addition & 1 deletion src/core/monica-parameters.h
Original file line number Diff line number Diff line change
Expand Up @@ -710,7 +710,7 @@ struct DLL_API CropModuleParameters : public Tools::Json11Serializable {
double pc_GrowthRespirationParameter2{ 0.0 };
double pc_Tortuosity{ 0.0 };
bool pc_AdjustRootDepthForSoilProps{ true };
int pc_TimeUnderAnoxiaThreshold{ 4 };
std::vector<int> pc_TimeUnderAnoxiaThreshold{ 4, 4, 4, 4, 4, 4, 4 };

bool __enable_Phenology_WangEngelTemperatureResponse__{ false };
bool __enable_Photosynthesis_WangEngelTemperatureResponse__{ false };
Expand Down
24 changes: 9 additions & 15 deletions src/io/build-output.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -546,11 +546,10 @@ BOTRes& monica::buildOutputTable()
return monica.cropGrowth() ? round(monica.cropGrowth()->get_PrimaryCropYield(), 1) : 0.0;
});

build({ id++, "SumYield", "kgDM ha-1", "get_AccumulatedPrimaryCropYield" },
[](const MonicaModel& monica, OId oid)
{
return monica.cropGrowth() ? round(monica.cropGrowth()->get_AccumulatedPrimaryCropYield(), 1) : 0.0;
});
build({id++, "SecondaryYield", "kgDM ha-1", "get_SecondaryCropYield"},
[](const MonicaModel &monica, OId oid) {
return monica.cropGrowth() ? round(monica.cropGrowth()->get_SecondaryCropYield(), 3) : 0.0;
});

build({ id++, "sumExportedCutBiomass", "kgDM ha-1", "return sum (across cuts) of exported cut biomass for current crop" },
[](const MonicaModel& monica, OId oid)
Expand Down Expand Up @@ -594,11 +593,6 @@ BOTRes& monica::buildOutputTable()
return round(monica.humusBalanceCarryOver(), 1);
});

build({ id++, "SecondaryYield", "kgDM ha-1", "get_SecondaryCropYield" },
[](const MonicaModel& monica, OId oid)
{
return monica.cropGrowth() ? round(monica.cropGrowth()->get_SecondaryCropYield(), 1) : 0.0;
});

build({ id++, "GroPhot", "kgCH2O ha-1", "GrossPhotosynthesisHaRate" },
[](const MonicaModel& monica, OId oid)
Expand Down Expand Up @@ -729,11 +723,11 @@ BOTRes& monica::buildOutputTable()
return Nstress;
});

build({ id++, "YieldNc", "kgN ha-1", "PrimaryYieldNConcentration" },
[](const MonicaModel& monica, OId oid)
{
return monica.cropGrowth() ? round(monica.cropGrowth()->get_PrimaryYieldNConcentration(), 3) : 0.0;
});
build({id++, "YieldNc", "kgN ha-1", "PrimaryYieldNConcentration"},
[](const MonicaModel &monica, OId oid) {
return monica.cropGrowth() ? round(monica.cropGrowth()->get_PrimaryYieldNConcentration(), 3) : 0.0;
});


build({ id++, "YieldN", "kgN ha-1", "PrimaryYieldNContent" },
[](const MonicaModel& monica, OId oid)
Expand Down
2 changes: 1 addition & 1 deletion src/resource/version.h
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@

#define VERSION_MAJOR 3
#define VERSION_MINOR 6
#define VERSION_REVISION 17
#define VERSION_REVISION 18
#define VERSION_BUILD 214

#define VER_FILE_VERSION VERSION_MAJOR, VERSION_MINOR, VERSION_REVISION, VERSION_BUILD
Expand Down
10 changes: 5 additions & 5 deletions src/run/cultivation-method.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -464,7 +464,7 @@ Errors Harvest::merge(json11::Json j) {
set_double_value(_optCarbMgmtData.maxResidueRecoverFraction, j, "max-residue-recover-fraction");

for (const string &organName: {"leaf", "shoot", "fruit", "struct", "sugar"}) {
for (auto kv: j.object_items()) {
for (const auto& kv: j.object_items()) {
if (toLower(kv.first) == organName && kv.second.is_object()) {
Spec::Value sv;
set_double_value(sv.exportPercentage, kv.second, "export");
Expand Down Expand Up @@ -493,8 +493,9 @@ json11::Json Harvest::to_json(bool includeFullCropParameters) const {
};

for (const auto &p: _spec.organ2specVal) {
jo[organNameFromId(p.first)] = J11Object{{"export", J11Array{p.second.exportPercentage, "%"}},
{"incorporate", p.second.incorporate}};
jo[organNameFromId(p.first)] = J11Object
{{"export", J11Array{p.second.exportPercentage, "%"}},
{"incorporate", p.second.incorporate}};
}

return jo;
Expand All @@ -505,8 +506,7 @@ bool Harvest::apply(MonicaModel *model) {

if (model->cropGrowth()) {
model->harvestCurrentCrop(_exported, _spec, _optCarbMgmtData);
if (_sowing)
debug() << "harvesting crop: " << _sowing->crop()->toString() << " at: " << date().toString() << endl;
if (_sowing) debug() << "harvesting crop: " << _sowing->crop()->toString() << " at: " << date().toString() << endl;
model->addEvent("Harvest");
}

Expand Down

0 comments on commit 11c0bb2

Please sign in to comment.