fixes for a segfault in test/modules/PenmanMonteith, updated (b1f7b460) · Commits · rootmodels / OpenSimRoot

OpenSimRoot/src/modules/Atmosphere/DiffusionResistance.cpp

+16 −7

Original line number	Diff line number	Diff line
		@@ -132,8 +132,11 @@ AerodynamicResistance::AerodynamicResistance(SimulaDynamic* const pSV) :
		//check if unit given in input file agrees with this function
		pSD->checkUnit("s/m");
		cropHeight_ = pSD->existingPath("/plants/canopyHeight", "cm");
		if (!cropHeight_) msg::warning("AerodynamicResistance: /plants/canopyHeight not found. Assuming 0.5 m");
		if (!cropHeight_) {
		msg::warning("AerodynamicResistance: /plants/canopyHeight not found. Assuming 0.5 m");
		}else{
		cropHeight_->checkUnit("cm");
		}
		windSpeed_ = pSD->existingPath("/environment/atmosphere/windSpeed"); //, "m/s"
		windSpeed_->checkUnit("m/s");
		}
		@@ -176,19 +179,25 @@ void AerodynamicResistance::calculate(const Time &t, double& r_a) {
		// Allen, Richard G., et al. "Crop evapotranspiration-Guidelines for computing crop water requirements-FAO Irrigation and drainage paper 56." Fao, Rome 300.9 (1998): D05109.
		// See also http://www.fao.org/3/X0490E/x0490e06.htm
		const double vonKarmanConstant = 0.41;
		double measurementHeight = 2;
		double cropHeight = 0.50; //m
		if (cropHeight_) {
		cropHeight_->get(t, cropHeight);
		cropHeight/=100; //cm to m
		cropHeight/=100.0; //cm to m
		if(cropHeight<0.1)cropHeight=0.1; //asume crop height is at least 10 cm, given that the soil will be uneven.
		}
		measurementHeight=std::max(measurementHeight,cropHeight+1);// make sure the measurement height is at least 1 m above the crop.
		const double measurementHeight=std::max(2.0,cropHeight+1.0);// make sure the measurement height is at least 1 m above the crop.

		const double zeroPlaneDisplacementHeight = 2./3. * cropHeight; // Convert from cm to m
		const double roughnessLengthMomentumTransfer = 0.123*cropHeight;
		const double roughnessLengthHeatAndVapour = 0.10.123cropHeight;
		r_a = log((measurementHeight - zeroPlaneDisplacementHeight)/roughnessLengthMomentumTransfer)log((measurementHeight - zeroPlaneDisplacementHeight)/roughnessLengthHeatAndVapour)/(vonKarmanConstantvonKarmanConstant*windSpeed);
		const double roughnessLengthHeatAndVapour = 0.1*roughnessLengthMomentumTransfer;

		const double zd=measurementHeight - zeroPlaneDisplacementHeight;
		const double zom=zd/roughnessLengthMomentumTransfer;
		const double zoh=zd/roughnessLengthHeatAndVapour;

		r_a = log(zom)log(zoh)/(vonKarmanConstantvonKarmanConstant*windSpeed);


		if (std::isnan(r_a)) msg::error("AerodynamicResistance: NaN value");
		if (std::isinf(r_a)) msg::error("AerodynamicResistance: inf value");
		if (r_a <= 0) msg::error("AerodynamicResistance: Aerodynamic resistance smaller than or equal to 0. Check if crop height is less than 2, if not, adjust measurement height in the code.");

OpenSimRoot/tests/modules/PenmanMonteith.xml

+786 −82

Original line number	Diff line number	Diff line
		@@ -30,10 +30,8 @@
		<SimulaBase name="plants">
		<SimulaTable name_column1="time" unit_column1="day"
		name="meanLeafAreaIndex" unit="cm2/cm2">
		0 0.
		10 5.
		20 5.
		30 0.
		0 5.
		3000 5.
		</SimulaTable>
		<SimulaConstant name="meanExtinctionCoefficient" unit="noUnit">
		0.6
		@@ -77,110 +75,816 @@
		</SimulaConstant>
		</SimulaDerivative>


		<SimulaDerivative name="interception" unit="cm/day"
		function="Interception" />

		<SimulaDerivative name="priestleyTaylor" unit="cm/day" function="PriestleyTaylor"
		<SimulaDerivative name="ET0_priestleyTaylor" unit="cm/day" function="PriestleyTaylor"
		/>
		<SimulaDerivative name="ET0_stanghellini" unit="cm/day" function="Stanghellini"
		/>
		<SimulaDerivative name="stanghellini" unit="cm/day" function="Stanghellini"
		<SimulaDerivative name="ET0_grassStandardReference" unit="cm/day" function="Grass_reference_evapotranspiration"
		/>
		<SimulaDerivative name="grassStandardReference" unit="cm/day" function="Grass_reference_evapotranspiration"
		<SimulaDerivative name="ET0_tallStandardReference" unit="cm/day" function="Tall_reference_Crop"
		/>
		<SimulaDerivative name="tallStandardReference" unit="cm/day" function="Tall_reference_Crop"
		<SimulaDerivative name="ET0_penmanEQ" unit="cm/day" function="penmanEQ"
		/>
		<SimulaDerivative name="penmanEQ" unit="cm/day" function="penmanEQ"
		<SimulaDerivative name="ET0_monteithEQ" unit="cm/day" function="monteithEQ"
		/>
		<SimulaDerivative name="monteithEQ" unit="cm/day" function="monteithEQ"
		<SimulaDerivative name="ET0_fromFile" unit="cm/day" function="useET0fromTable"
		/>

		</SimulaBase>


		<SimulaBase name="environment">
		<SimulaConstant name="startDay" unit="noUnit" type="integer"> 1
		<SimulaConstant
		name="altitude"
		unit="m"> 64
		</SimulaConstant>
		<SimulaConstant
		name="startDay"
		unit="noUnit"
		type="integer"> 25
		</SimulaConstant>
		<SimulaConstant
		name="startMonth"
		unit="noUnit"
		type="integer"> 3
		</SimulaConstant>
		<SimulaConstant
		name="latitude"
		unit="noUnit"> 50.624071
		</SimulaConstant>
		<SimulaConstant name="startMonth" unit="noUnit" type="integer"> 1
		<SimulaConstant
		name="longitude"
		unit="noUnit"> 6.985767
		</SimulaConstant>

		<SimulaBase name="atmosphere">
		<SimulaTable name_column1="time" name_column2="irradiation"
		unit_column1="day" unit_column2="umol/cm2/day">
		0 3000
		100 3000
		<SimulaTable
		name_column1="time"
		name_column2="shortwaveRadiation"
		unit_column1="day"
		unit_column2="MJ/m2/day">
		0 13.59
		1 8.77
		2 18.42
		3 18.22
		4 17.9
		5 8.96
		6 10.22
		7 15.49
		8 17.51
		9 19.56
		10 20.11
		11 19.21
		12 21.2
		13 21.03
		14 20.95
		15 20.59
		16 12
		17 21.6
		18 21.76
		19 14.8
		20 17.54
		21 14.76
		22 5.01
		23 4.76
		24 8.5
		25 18.75
		26 18.55
		27 18.71
		28 13.62
		29 13.95
		30 5.98
		31 14.87
		32 23.59
		33 24.17
		34 24.9
		35 25.26
		36 24.72
		37 25.21
		38 22.69
		39 14.14
		40 15.72
		41 16.01
		42 18.47
		43 13.09
		44 22.3
		45 24.25
		46 26.79
		47 26.82
		48 26.49
		49 27.35
		50 26.24
		51 21.98
		52 25.42
		53 16.95
		54 13.68
		55 25.8
		56 21.98
		57 20.7
		58 21.94
		59 20.37
		60 8.99
		61 14.12
		62 20.53
		63 13.95
		64 18.5
		65 10.22
		66 25.53
		67 18.64
		68 18.63
		69 22.75
		70 20.83
		71 19.65
		72 16.11
		73 20.3
		74 19.9
		75 24.45
		76 19.34
		77 16.39
		78 27.6
		79 29.4
		80 27.8
		81 22.43
		82 16.06
		83 21.44
		84 28.45
		85 27.6
		86 29.06
		87 29.74
		88 29.42
		89 26.1
		90 19.68
		91 24.78
		92 23.93
		93 20.02
		94 14.4
		95 25.06
		96 26.28
		97 28.65
		98 26.94
		</SimulaTable>
		<SimulaTable name_column2="precipitation" name_column1="time"
		unit_column1="day" unit_column2="cm/day" interpolationMethod="step">
		<SimulaTable
		name_column2="precipitation"
		name_column1="time"
		unit_column1="day"
		unit_column2="cm/day"
		interpolationMethod="step">
		0 0
		1 0
		1 0.31
		2 0
		3 0.056666667
		4 0.26
		5 0.053333333
		6 0.203333333
		7 0.046666667
		8 0.19
		3 0
		4 0.07
		5 0.02
		6 0.08
		7 0
		8 0
		9 0
		10 0.19
		10 0
		11 0
		30 0
		</SimulaTable>
		<!-- for pennman monteith -->
		<SimulaConstant name="netRadiationSoil" unit="W/m2">
		12 0
		13 0
		14 0
		15 0
		</SimulaConstant>
		<SimulaConstant name="netRadiation" unit="W/m2">
		300
		</SimulaConstant>
		<SimulaConstant name="windSpeed" unit="m/s">
		1.2
		</SimulaConstant>
		<SimulaConstant name="relativeHumidity" unit="m/s"> 65
		</SimulaConstant>
		<SimulaTable name_column2="averageDailyTemperature"
		name_column1="time" unit_column1="day" unit_column2="degreesC"
		16 0
		17 0
		18 0
		19 0.52
		20 0.01
		21 0.06
		22 0.41
		23 0.86
		24 0.67
		25 0.5
		26 0.58
		27 0.04
		28 0.74
		29 1.05
		30 0.82
		31 0.25
		32 0
		33 0
		34 0
		35 0
		36 0
		37 0
		38 0.08
		39 0.99
		40 0.09
		41 0.03
		42 0.03
		43 0
		44 0
		45 0
		46 0
		47 0
		48 0
		49 0
		50 0
		51 0.01
		52 0
		53 0.05
		54 0.1
		55 0
		56 0
		57 0
		58 0.19
		59 0.09
		60 0.72
		61 0.42
		62 0.02
		63 0.13
		64 0.62
		65 0.6
		66 0
		67 0.48
		68 0.22
		69 0.07
		70 0
		71 0.1
		72 0.28
		73 0.53
		74 0.17
		75 0.85
		76 0.52
		77 0.02
		78 0
		79 0
		80 0
		81 0.54
		82 0.77
		83 0.5
		84 0
		85 0
		86 0
		87 0
		88 0
		89 0
		90 0.03
		91 0
		92 0
		93 0.94
		94 0.05
		95 0
		96 0
		97 0
		98 0

		</SimulaTable>
		<SimulaTable
		name="windSpeed"
		unit="m/s"
		name_column1="time"
		unit_column1="day">
		0 15.8
		1 15.5
		2 4.7
		3 15.3
		4 17.9
		5 29.3
		6 19.5
		7 23.4
		8 25.4
		9 17.8
		10 7.5
		11 16
		12 19.1
		13 9.5
		14 9.2
		15 13.6
		16 13
		17 9.6
		18 16.6
		19 21
		20 10.9
		21 15
		22 14
		23 22.2
		24 16
		25 11.9
		26 16.4
		27 15.5
		28 11.4
		29 15.2
		30 11.3
		31 13.7
		32 12.8
		33 7.8
		34 8.7
		35 6.8
		36 5.4
		37 11.3
		38 14.7
		39 18.4
		40 15.1
		41 15.2
		42 16.3
		43 12.6
		44 15.1
		45 16.2
		46 11.1
		47 14.3
		48 13.9
		49 11.9
		50 14.9
		51 18.6
		52 18.9
		53 15.2
		54 16.5
		55 8
		56 11.5
		57 15.4
		58 22.5
		59 19.8
		60 20.7
		61 21.8
		62 18
		63 25.8
		64 25.6
		65 16.9
		66 14.8
		67 15
		68 16.5
		69 16.3
		70 18.2
		71 19.8
		72 24.2
		73 27.3
		74 21.7
		75 29
		76 12.7
		77 25.2
		78 10.9
		79 19.1
		80 14.8
		81 20
		82 17.7
		83 12.5
		84 6.5
		85 14.8
		86 13.9
		87 10
		88 12.1
		89 16.1
		90 29
		91 24.6
		92 12.8
		93 22.7
		94 20.8
		95 17
		96 15.2
		97 10.8
		98 8.7

		</SimulaTable>
		<SimulaTable
		name="relativeHumidity"
		unit="%"
		name_column1="time"
		unit_column1="day"
		>
		0 82
		1 82
		2 77
		3 70
		4 72
		5 67
		6 70
		7 68
		8 64
		9 63
		10 65
		11 65
		12 52
		13 60
		14 64
		15 68
		16 77
		17 76
		18 62
		19 73
		20 71
		21 77
		22 92
		23 92
		24 84
		25 75
		26 79
		27 77
		28 86
		29 86
		30 93
		31 82
		32 72
		33 68
		34 72
		35 65
		36 67
		37 65
		38 63
		39 87
		40 75
		41 70
		42 72
		43 72
		44 61
		45 58
		46 56
		47 55
		48 48
		49 49
		50 52
		51 56
		52 54
		53 61
		54 71
		55 60
		56 60
		57 56
		58 61
		59 59
		60 68
		61 82
		62 65
		63 66
		64 80
		65 82
		66 71
		67 79
		68 77
		69 67
		70 65
		71 67
		72 75
		73 73
		74 72
		75 71
		76 66
		77 69
		78 69
		79 54
		80 47
		81 56
		82 72
		83 64
		84 61
		85 58
		86 57
		87 54
		88 45
		89 45
		90 54
		91 54
		92 56
		93 71
		94 70
		95 65
		96 67
		97 59
		98 43

		</SimulaTable>
		<SimulaTable
		name_column2="averageDailyTemperature"
		name_column1="time"
		unit_column1="day"
		unit_column2="degreesC"
		interpolationMethod="linear">
		0 5.4
		1 5.8
		2 6.5
		3 7.3
		0 8.5
		1 8.3
		2 7.7
		3 9.5
		4 8.4
		5 10.9
		6 7.4
		7 9
		8 5.8
		9
		5.4
		10 5.9
		30 30
		5 8.4
		6 7.9
		7 6.9
		8 9.8
		9 12.3
		10 12.4
		11 11.3
		12 6.4
		13 5.7
		14 7.7
		15 8.7
		16 7.3
		17 9.1
		18 13
		19 14.5
		20 13
		21 13.6
		22 10.4
		23 8
		24 7.8
		25 9.9
		26 11.5
		27 11
		28 11.1
		29 11.5
		30 10.2
		31 11
		32 11.5
		33 13
		34 13.6
		35 14.4
		36 16.3
		37 17.9
		38 20
		39 15.4
		40 10
		41 9.4
		42 10.3
		43 10.7
		44 10.6
		45 11
		46 12.9
		47 15.3
		48 16.8
		49 15.7
		50 16
		51 13.3
		52 12.7
		53 13.6
		54 11.9
		55 14.4
		56 15.9
		57 14.4
		58 11.6
		59 9.8
		60 10.5
		61 15.4
		62 14.5
		63 15.4
		64 15.6
		65 14.7
		66 19.9
		67 20.6
		68 19.1
		69 16.3
		70 17.1
		71 18.2
		72 16.1
		73 16.7
		74 15.7
		75 13.7
		76 13.9
		77 16.2
		78 15.4
		79 19.9
		80 24.9
		81 25.2
		82 19.2
		83 18.4
		84 20.4
		85 21.3
		86 20.7
		87 19.6
		88 22.2
		89 25.6
		90 20.5
		91 19.8
		92 24.4
		93 21.9
		94 19.6
		95 23.2
		96 23.6
		97 25.2
		98 28.4


		</SimulaTable>

		<!-- these are necessary to compute incoming radiation by the radiation
		module -->
		<SimulaConstant name="albedoSoil" unit="noUnit">
		0.17
		</SimulaConstant>
		<SimulaConstant name="albedoCrop" unit="noUnit">
		0.23
		</SimulaConstant>
		<SimulaConstant name="altitude" unit="m">
		64
		</SimulaConstant>
		<SimulaConstant name="startDay" unit="noUnit" type="integer">
		25
		<SimulaConstant
		name="albedoSoil"
		unit="noUnit"> 0.17
		</SimulaConstant>
		<SimulaConstant name="startMonth" unit="noUnit" type="integer">
		3
		<SimulaConstant
		name="albedoCrop"
		unit="noUnit"> 0.23
		</SimulaConstant>
		<SimulaConstant name="latitude" unit="noUnit">
		50.8
		</SimulaConstant>
		<SimulaTable name_column2="actualDurationofSunshine"
		name_column1="time" unit_column1="day" unit_column2="hour"
		<!-- 0.23 grass 0.17 bare soil 0.05 open water 0.4 desert sand -->
		<SimulaTable
		name_column2="actualDurationofSunshine"
		name_column1="time"
		unit_column1="day"
		unit_column2="hour"
		interpolationMethod="step">
		<!-- time of the day that there are no clouds -->
		0 7.78043333333333
		1 2.42387777777778
		2 12.0575
		3 12.1541944444444
		4 10.9080861111111
		5 5.62731944444444
		6 6.46415277777778
		7 12.7703
		8 12.5184777777778
		9 12.8552194444444
		10 12.8918972222222
		11 12.2940583333333
		12 12.9651138888889
		13 13
		14 13
		15 13
		16 7.143825
		17 13
		18 13
		19 10.6974027777778
		20 11.3215888888889
		21 8.34925
		22 0
		23 0
		24 1.81589166666667
		25 11.6806888888889
		26 10.1578916666667
		27 11.3243277777778
		28 8.20366666666667
		29 6.06355555555556
		30 0
		31 9.81924444444445
		32 13.415175
		33 13.6029861111111
		34 14.0061861111111
		35 13.9823361111111
		36 14.2820527777778
		37 14.571325
		38 14.2381194444444
		39 7.12981388888889
		40 12.1244111111111
		41 10.4050111111111
		42 11.9504888888889
		43 5.14893333333333
		44 12.5982055555556
		45 14.9864805555556
		46 15
		47 15
		48 15
		49 15
		50 15
		51 11
		52 13.5908277777778
		53 11.7929361111111
		54 7.00061388888889
		55 15
		56 15
		57 12.2591111111111
		58 13.7334361111111
		59 14
		60 3.3232
		61 6.85596111111111
		62 14.6261944444444
		63 5.13385
		64 10.2347333333333
		65 2.87796111111111
		66 15
		67 9.18990555555556
		68 10.8661861111111
		69 15
		70 14.33815
		71 11.7950388888889
		72 9.51720555555556
		73 11.4338083333333
		74 11
		75 14.8591888888889
		76 10.8227527777778
		77 6.45055555555556
		78 15.0605888888889
		79 15.1460083333333
		80 15.0514194444444
		81 13
		82 7.05885833333333
		83 13.1488694444444
		84 15.1620972222222
		85 15.1594333333333
		86 15.1570361111111
		87 15.2457111111111
		88 14.3951916666667
		89 13.1222777777778
		90 10.9601
		91 12.8373333333333
		92 13.7843833333333
		93 11.9497555555556
		94 5.81592222222222
		95 15.1494555555556
		96 15.1502833333333
		97 15.1514611111111
		98 15.1570583333333

		</SimulaTable>
		<SimulaTable
		name_column2="ET0"
		unit_column2="cm/day"
		name_column1="time"
		unit_column1="day"
		interpolationMethod="step">
		0 7.66
		129 7.89

		0 0.176
		1 0.129
		2 0.235
		3 0.274
		4 0.241
		5 0.198
		6 0.177
		7 0.229
		8 0.295
		9 0.33
		10 0.323
		11 0.31
		12 0.297
		13 0.27
		14 0.286
		15 0.301
		16 0.17
		17 0.294
		18 0.378
		19 0.274
		20 0.298
		21 0.252
		22 0.081
		23 0.077
		24 0.124
		25 0.269
		26 0.287
		27 0.28
		28 0.204
		29 0.202
		30 0.097
		31 0.22
		32 0.361
		33 0.371
		34 0.391
		35 0.408
		36 0.415
		37 0.47
		38 0.468
		39 0.228
		40 0.227
		41 0.245
		42 0.275
		43 0.216
		44 0.349
		45 0.382
		46 0.442
		47 0.494
		48 0.532
		49 0.49
		50 0.497
		51 0.412
		52 0.431
		53 0.327
		54 0.241
		55 0.426
		56 0.397
		57 0.394
		58 0.356
		59 0.332
		60 0.209
		61 0.251
		62 0.366
		63 0.328
		64 0.306
		65 0.197
		66 0.476
		67 0.361
		68 0.341
		69 0.397
		70 0.426
		71 0.408
		72 0.306
		73 0.373
		74 0.351
		75 0.368
		76 0.331
		77 0.353
		78 0.449
		79 0.624
		80 0.693
		81 0.551
		82 0.337
		83 0.415
		84 0.527
		85 0.546
		86 0.58
		87 0.561
		88 0.641
		89 0.703
		90 0.501
		91 0.549
		92 0.531
		93 0.402
		94 0.323
		95 0.552
		96 0.543
		97 0.617
		98 0.626

		</SimulaTable>
		</SimulaBase>

		</SimulaBase>

		<SimulaBase name="rootTypeParameters">
		@@ -197,7 +901,7 @@
		0.
		</SimulaConstant>
		<SimulaConstant name="endTime" type="time">
		30.
		98
		</SimulaConstant>
		<SimulaConstant name="timeInterval" type="time">
		1.