Compute moisture rmf as in official RothC26C3

[hbessler]

This commit fixes a discrepancy between the official RothC and SoilR implementation in the calculation of the soil moisture - rate modifyer (rmf) under bare soil.

In the official RothC implementation (https://github.com/Rothamsted-Models/RothC_Code):

• maximum topsoil moisture deficit (Max.TSMD) is calculated from clay content and thickness of the topsoil,

• accumulated soil moisture deficit (Acc.TSMD) is limited to the bare soil moisture deficit (BareSMD = Max.TSMD / 1.8), and

• rate modifyer (beta) is calculated from Max.TSMD and Acc.TSMD.

In the SoilR implementation:

• Max.TSMD is calculated from clay content and thickness of the topsoil as in the official RothC implementation but set to Max.TSMD / 1.8 (= BareSMD),

• Acc.TSMD is limited to Max.TSMD, and - beta is calculated from Max.TSMD and Acc.TSMD,

which leads to wrong results.

In this commit:

• Max.TSMD and b are calculated as in the official RothC implementation, and

• Acc.TSMD is limited to TSMD.limit, which is Max.TSMD / 1.8 for bare and Max.TSMD for vegetated soil.

[crlsierra]

Thanks for pointing this potential issue. I checked the code, and I don't see a major difference with the previous code, except that the bare soil moisture deficit modifier is computed separately and not as part of the Acc.TSMD as I have done previously.

Can you provide a test example that shows how the previous implementation fails and produce results from the original RothC model? The function currently implemented in the package reproduces the Table provide in the original documentation, and this test is presented in the example of the SoilR function. Can you provide two separate tests, one that shows how the SoilR implementation fails to reproduce results from the original RothC model, and another test that shows that your function produces the same results as the Table in the RothC documentation?

[hbessler]

As a test, we can use the data from the original RothC documentation and calculate Acc.TSMD and the rate modifying factor (beta) for bareground (bare = TRUE) using [1] original RothC code, [2] SoilR, [3] provided fix for SoilR.

The results are given in the table below. They show:

1. Minor differences in Acc.TSMD between [1] and [2],[3] from Jun to Aug, when the max. bare soil moisture deficit (BareSMD) is reached. This is due to rounding error because BareSMD is calculated as Max.TSMD / 1.8 in [2] and [3], which is consistent with the official RothC documentation, but as 0.556 x Max.TSMD in [1]. This discrepancy is not part of the commited fix.
2. Major differences in beta between [1] and [2] from Jun to Sep. Note that the results of [2] are not consistent with those in Figure 3 in the original RothC documentation. This discrepancy is resolved with the commited fix [3].

Month	P	E	Acc.TSMD[1]	Acc.TSMD[2]	Acc.TSMD[3]	beta[1]	beta [2]	beta [3]
Jan	74	8	0	0	0	1.00	1.00	1.00
Feb	59	10	0	0	0	1.00	1.00	1.00
Mar	62	27	0	0	0	1.00	1.00	1.00
Apr	51	49	0	0	0	1.00	1.00	1.00
May	52	83	-10.25	-10.25	-10.25	1.00	1.00	1.00
Jun	57	99	-24.99	-24.97	-24.97	0.84	0.20	0.84
Jul	34	103	-24.99	-24.97	-24.97	0.84	0.20	0.84
Aug	55	91	-24.99	-24.97	-24.97	0.84	0.20	0.84
Sep	58	69	-18.74	-18.74	-18.74	1.00	0.56	1.00
Oct	56	34	0	0	0	1.00	1.00	1.00
Nov	75	16	0	0	0	1.00	1.00	1.00
Dec	71	8	0	0	0	1.00	1.00	1.00

P...Precipitation, E...Open pan evaporation

[crlsierra]

OK, great. This looks convincing. I will merge the pull request and the changes will be available in the latest SoilR version in GitHub.