30 day months are assumed for f_frost

[twest820]

At the moment, f_frost is calculated as

! frost modifier
f_frost(:,i) = 1.d0 - kF(i) * ( frost_days(:) / 30.d0)

which I think matches equation A6 of the 3-PGmix user manual (the manual defines fF in terms of a month's mean frost days but I presume the "mean" bit is a typo). The 3-PGmix manual cites Landsberg and Waring 1997 for this but the corresponding equation in Landsberg and Waring is

fT = 1 - (frost days per month) / (number of days per month)

Under what seems to be the default assumption of kF = 1, treating all months as having 30 days appears problematic, particularly on colder sites.

1. f_frost is biased 3.3% low in months with 31 days and 6.7% or 3.3% high in February.
2. f_frost cannot go to zero for February, even if the whole month is below freezing.
3. f_frost can become negative in 31 day months.

Since kF is a scalar, no value of kF exists which avoids both of these problems. Additionally, choosing kF = 30/31 avoid negative f_frost in 31 day months expands the February problem to all 30 day months.

This makes me suspect both the 3-PGmix manual and Fortran code have a bug and that Landsberg and Waring's use of the number of days in a given month would be preferred. It also seems like, particularly in cases where kF ≠ 1, f_frost should be not be allowed to move out of the [0, 1] range.

Am I missing something?

[trotsiuk]

Hi @twest820

Thank you for finding this interesting issue. Indeed the manual of the 3-PG is declaring the number of days per month.

However, if you look at the code itself (also of 3-PGpjs which can be downloaded here https://3pg.forestry.ubc.ca/software/) it will be recognised that the same formula as in r3PG or 3PGmix is used for the frost days, see below:

'calculate frost modifier
 fFrost = 1 - kF * (FrostDays / 30)

Since we were translating the code and aimed to test compatibility it was done the same way as in the original model.

I fully support your suggestions for improvements. Please feel free to create a new branch an make a push request to incorporate the changes.

Unfortunately from my side there is very little time to allocate for it.

3PGpjs27.zip

[twest820]

Yeah, I know how that goes. Thanks Volodymyr, I'll address this in my branch but, as I'm not coding in Fortran, a pull request is unlikely to materialize. A daysInMonth(month) lookup, as done elsewhere in the existing code, would address this. Modulo the leap year issue I'm about to capture, anyways.

[trotsiuk]

I'll address this in my branch but, as I'm not coding in Fortran

So is your intension to write model fully in R or other language? Would be great if you implement all your detailed comments and suggestions into the code!

[twest820]

I think it's most accurate to say I'm looking into an experimental version for evaluating certain possibilities. Definitely not funded for anything like an R package though it's possible that, if things go well, the code might eventually become of wider interest. If that happens I'll follow up with you (most likely in email) but I anticipate it will be at least several months to find out, possibly well into next year.

[trotsiuk]

I have now implemented the number of days per month to calculate the f_frost

r3PG/pkg/src/md_3PG.f95

Line 88 in 40289a7

f_frost(:,i) = 1.d0 - kF(i) * ( frost_days(:) / daysInMonth(month_vector(:)) )

Additionally add the check to the prepare_climate function to limit the frost_day to the maximum possible

r3PG/pkg/R/prepare_climate.R

Line 100 in 40289a7

climate$frost_days <- pmin( climate$frost_days, daysInMonth[climate$month])

@DavidForrester this slightly however impacted the output of the model. Particularly it has deviated the biomass outputs by third digit for the stem biomass. The general pattern remain unchanged.