diff --git a/docs/source/_toc.yml b/docs/source/_toc.yml
index 1c587e0e..73c15779 100644
--- a/docs/source/_toc.yml
+++ b/docs/source/_toc.yml
@@ -28,6 +28,7 @@ subtrees:
         - file: users/pmodel/c3c4model.md
         - file: users/pmodel/isotopic_discrimination.md
   - file: users/tmodel/tmodel.md
+  - file: users/tmodel/canopy.md
   - file: users/splash.md
   - file: users/constants.md
   - file: users/hygro.md 
diff --git a/docs/source/api/demography_api.md b/docs/source/api/demography_api.md
new file mode 100644
index 00000000..1244feff
--- /dev/null
+++ b/docs/source/api/demography_api.md
@@ -0,0 +1,28 @@
+---
+jupytext:
+  formats: md:myst
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 0.13
+kernelspec:
+  display_name: Python 3
+  language: python
+  name: python3
+---
+
+# The {mod}`~pyrealm.demography` module
+
+```{eval-rst}
+.. automodule:: pyrealm.demography
+    :autosummary:
+    :members:
+```
+
+## The {mod}`~pyrealm.demography.flora` module
+
+```{eval-rst}
+.. automodule:: pyrealm.demography.flora
+    :autosummary:
+    :members:
+```
diff --git a/docs/source/refs.bib b/docs/source/refs.bib
index 700ff609..0dcb8e4a 100644
--- a/docs/source/refs.bib
+++ b/docs/source/refs.bib
@@ -82,6 +82,22 @@ @article{BerberanSantos:2009bk
   file = {/Users/dorme/Zotero/storage/EN33QA2H/J Math Chem 2009 Berberan-Santos.pdf}
 }
 
+@article{berger:1993a,
+  title = {Insolation and {{Earth}}'s Orbital Periods},
+  author = {Berger, Andr{\'e} and Loutre, Marie-France and Tricot, Christian},
+  year = {1993},
+  journal = {Journal of Geophysical Research: Atmospheres},
+  volume = {98},
+  number = {D6},
+  pages = {10341--10362},
+  issn = {2156-2202},
+  doi = {10.1029/93JD00222},
+  urldate = {2024-08-02},
+  abstract = {Solar irradiance received on a horizontal surface depends on the solar output, the semimajor axis of the elliptical orbit of the Earth around the sun (a), the distance from the Earth to the sun (r), and the zenith distance (z). The spectrum of the distance, r, for a given value of the true longitude, {$\lambda$}, displays mainly the precessional periods and, with much less power, half precession periods, eccentricity periods, and some combination tones. The zenith distance or its equivalent, the elevation angle (E), is only a function of obliquity ({$\epsilon$}) for a given latitude, {$\phi$}, true longitude, and hour angle, H. Therefore the insolation at a given constant value of z is only a function of precession and eccentricity. On the other hand, the value of the hour angle, H, corresponding to this fixed value of z varies with {$\varepsilon$}, except for the equinoxes, where H corresponding to a constant z also remains constant through time. Three kinds of insolation have been computed both analytically and numerically: the instantaneous insolation (irradiance) at noon, the daily irradiation, and the irradiations received during particular time intervals of the day defined by two constant values of the zenith distance (diurnal irradiations). Mean irradiances (irradiations divided by the length of the time interval over which they are calculated) are also computed for different time intervals, like the interval between sunrise and sunset, in particular. Examples of these insolations are given in this paper for the equinoxes and the solstices. At the equinoxes, for each latitude, all insolations are only a function of precession (this invalidates the results obtained by Cerveny [1991]). At the solstices, both precession and obliquity are present, although precession dominates for most of the latitudes. Because the lengths of the astronomical seasons are secularly variable (in terms of precession only), a particular calendar day does not always correspond to the same position relative to the sun through geological time. Similarly, a given longitude of the Sun on its orbit does not correspond to the same calendar day. For example, 103 kyr ago, assuming arbitrarily that the spring equinox is always on March 21, autumn began on September 13, and 114 kyr ago, it began on September 27, the length of the summer season being 85 and 98 calendar days, respectively, at these remote times in the past.},
+  langid = {english},
+  file = {/Users/dorme/Zotero/storage/ALTXDVSM/Berger et al. - 1993 - Insolation and Earth's orbital periods.pdf;/Users/dorme/Zotero/storage/HWQ4XSZN/93JD00222.html}
+}
+
 @article{berger:1978a,
   title = {Long-{{Term Variations}} of {{Daily Insolation}} and {{Quaternary Climatic Changes}}},
   author = {Berger, Andr{\'e}L},
@@ -101,20 +117,25 @@ @article{berger:1978a
   file = {/Users/dorme/Zotero/storage/DPM5WN32/Berger - 1978 - Long-Term Variations of Daily Insolation and Quate.pdf}
 }
 
-@article{berger:1993a,
-  title = {Insolation and {{Earth}}'s Orbital Periods},
-  author = {Berger, Andr{\'e} and Loutre, Marie-France and Tricot, Christian},
-  year = {1993},
-  journal = {Journal of Geophysical Research: Atmospheres},
-  volume = {98},
-  number = {D6},
-  pages = {10341--10362},
-  issn = {2156-2202},
-  doi = {10.1029/93JD00222},
-  urldate = {2024-08-02},
-  abstract = {Solar irradiance received on a horizontal surface depends on the solar output, the semimajor axis of the elliptical orbit of the Earth around the sun (a), the distance from the Earth to the sun (r), and the zenith distance (z). The spectrum of the distance, r, for a given value of the true longitude, {$\lambda$}, displays mainly the precessional periods and, with much less power, half precession periods, eccentricity periods, and some combination tones. The zenith distance or its equivalent, the elevation angle (E), is only a function of obliquity ({$\epsilon$}) for a given latitude, {$\phi$}, true longitude, and hour angle, H. Therefore the insolation at a given constant value of z is only a function of precession and eccentricity. On the other hand, the value of the hour angle, H, corresponding to this fixed value of z varies with {$\varepsilon$}, except for the equinoxes, where H corresponding to a constant z also remains constant through time. Three kinds of insolation have been computed both analytically and numerically: the instantaneous insolation (irradiance) at noon, the daily irradiation, and the irradiations received during particular time intervals of the day defined by two constant values of the zenith distance (diurnal irradiations). Mean irradiances (irradiations divided by the length of the time interval over which they are calculated) are also computed for different time intervals, like the interval between sunrise and sunset, in particular. Examples of these insolations are given in this paper for the equinoxes and the solstices. At the equinoxes, for each latitude, all insolations are only a function of precession (this invalidates the results obtained by Cerveny [1991]). At the solstices, both precession and obliquity are present, although precession dominates for most of the latitudes. Because the lengths of the astronomical seasons are secularly variable (in terms of precession only), a particular calendar day does not always correspond to the same position relative to the sun through geological time. Similarly, a given longitude of the Sun on its orbit does not correspond to the same calendar day. For example, 103 kyr ago, assuming arbitrarily that the spring equinox is always on March 21, autumn began on September 13, and 114 kyr ago, it began on September 27, the length of the summer season being 85 and 98 calendar days, respectively, at these remote times in the past.},
+@article{Bernacchi:2003dc,
+  title = {In Vivo Temperature Response Functions of Parameters Required to Model {{RuBP-limited}} Photosynthesis},
+  author = {Bernacchi, C J and Pimentel, C and Long, S P},
+  year = {2003},
+  month = sep,
+  journal = {Plant, Cell \& Environment},
+  volume = {26},
+  number = {9},
+  pages = {1419--1430},
+  publisher = {John Wiley \& Sons, Ltd},
+  doi = {10.1046/j.0016-8025.2003.01050.x},
+  abstract = {The leaf model of C3 photosynthesis of Farquhar, von Caemmerer \& Berry (Planta 149, 78--90, 1980) provides the basis for scaling carbon exchange from leaf to canopy and Earth-System models, and is wid...},
+  date-added = {2020-11-30T14:15:56GMT},
+  date-modified = {2020-11-30T14:42:52GMT},
   langid = {english},
-  file = {/Users/dorme/Zotero/storage/ALTXDVSM/Berger et al. - 1993 - Insolation and Earth's orbital periods.pdf;/Users/dorme/Zotero/storage/HWQ4XSZN/93JD00222.html}
+  local-url = {file://localhost/Users/dorme/References/Library.papers3/Articles/2003/Bernacchi/Plant\%20Cell\%20\&\%20Environment\%202003\%20Bernacchi.pdf},
+  rating = {0},
+  uri = {papers3://publication/doi/10.1046/j.0016-8025.2003.01050.x},
+  file = {/Users/dorme/Zotero/storage/7J2I98DM/Plant Cell & Environment 2003 Bernacchi.pdf}
 }
 
 @article{Bernacchi:2001kg,
@@ -136,27 +157,6 @@ @article{Bernacchi:2001kg
   file = {/Users/dorme/Zotero/storage/FI9562Z3/Plant Cell & Environment 2001 Bernacchi.pdf}
 }
 
-@article{Bernacchi:2003dc,
-  title = {In Vivo Temperature Response Functions of Parameters Required to Model {{RuBP-limited}} Photosynthesis},
-  author = {Bernacchi, C J and Pimentel, C and Long, S P},
-  year = {2003},
-  month = sep,
-  journal = {Plant, Cell \& Environment},
-  volume = {26},
-  number = {9},
-  pages = {1419--1430},
-  publisher = {John Wiley \& Sons, Ltd},
-  doi = {10.1046/j.0016-8025.2003.01050.x},
-  abstract = {The leaf model of C3 photosynthesis of Farquhar, von Caemmerer \& Berry (Planta 149, 78--90, 1980) provides the basis for scaling carbon exchange from leaf to canopy and Earth-System models, and is wid...},
-  date-added = {2020-11-30T14:15:56GMT},
-  date-modified = {2020-11-30T14:42:52GMT},
-  langid = {english},
-  local-url = {file://localhost/Users/dorme/References/Library.papers3/Articles/2003/Bernacchi/Plant\%20Cell\%20\&\%20Environment\%202003\%20Bernacchi.pdf},
-  rating = {0},
-  uri = {papers3://publication/doi/10.1046/j.0016-8025.2003.01050.x},
-  file = {/Users/dorme/Zotero/storage/7J2I98DM/Plant Cell & Environment 2003 Bernacchi.pdf}
-}
-
 @article{boyd:2015a,
   title = {Temperature Response of {{C4}} Photosynthesis: {{Biochemical}} Analysis of {{Rubisco}}, {{Phosphoenolpyruvate Carboxylase}} and {{Carbonic Anhydrase}} in {{Setaria}} Viridis.},
   shorttitle = {Temperature Response of {{C4}} Photosynthesis},
@@ -267,19 +267,6 @@ @article{Diaz:2016by
   file = {/Users/dorme/Zotero/storage/RGG46LIG/Nature 2016 Díaz.pdf}
 }
 
-@book{duffie:2013a,
-  title = {Solar {{Engineering}} of {{Thermal Processes}}},
-  author = {Duffie, John A. and Beckman, William A.},
-  year = {2013},
-  month = apr,
-  publisher = {John Wiley \& Sons},
-  abstract = {The updated fourth edition of the "bible" of solar energy theory and applications Over several editions, Solar Engineering of Thermal Processes has become a classic solar engineering text and reference. This revised Fourth Edition offers current coverage of solar energy theory, systems design, and applications in different market sectors along with an emphasis on solar system design and analysis using simulations to help readers translate theory into practice. An important resource for students of solar engineering, solar energy, and alternative energy as well as professionals working in the power and energy industry or related fields, Solar Engineering of Thermal Processes, Fourth Edition features:  Increased coverage of leading-edge topics such as photovoltaics and the design of solar cells and heaters A brand-new chapter on applying CombiSys (a readymade TRNSYS simulation program available for free download) to simulate a solar heated house with solar- heated domestic hot water Additional simulation problems available through a companion website An extensive array of homework problems and exercises},
-  googlebooks = {5uDdUfMgXYQC},
-  isbn = {978-1-118-41541-2},
-  langid = {english},
-  keywords = {Technology & Engineering / Electronics / General,Technology & Engineering / Mechanical,Technology & Engineering / Power Resources / General}
-}
-
 @article{Farquhar:1980ft,
   title = {A Biochemical Model of Photosynthetic {{CO2}} Assimilation in Leaves of {{C3}} Species.},
   author = {Farquhar, G D and {von Caemmerer}, S and Berry, J A},
@@ -502,6 +489,14 @@ @phdthesis{Lancelot:2020tm
   file = {/Users/dorme/Zotero/storage/86NBGKCA/2020 Lancelot.pdf}
 }
 
+@article{lavergne:2022a,
+  title = {A Semi-Empirical Model for Primary Production, Isotopic Discrimination and Competition of {{C3}} and {{C4}} Plants},
+  author = {Lavergne, Ali{\'e}nor and Harrison, Sandy P. and Atsawawaranunt, Kamolphat and Dong, Ning and Prentice, Iain Colin},
+  year = {2022},
+  journal = {Global Ecology and Biogeography (submitted)},
+  file = {/Users/dorme/Zotero/storage/I797PRT2/Lavergneetal_GEB.docx}
+}
+
 @article{lavergne:2020a,
   title = {Impacts of Soil Water Stress on the Acclimated Stomatal Limitation of Photosynthesis: {{Insights}} from Stable Carbon Isotope Data},
   shorttitle = {Impacts of Soil Water Stress on the Acclimated Stomatal Limitation of Photosynthesis},
@@ -520,14 +515,6 @@ @article{lavergne:2020a
   file = {/Users/dorme/Zotero/storage/ECFRUWX5/Lavergne et al. - 2020 - Impacts of soil water stress on the acclimated sto.pdf}
 }
 
-@article{lavergne:2022a,
-  title = {A Semi-Empirical Model for Primary Production, Isotopic Discrimination and Competition of {{C3}} and {{C4}} Plants},
-  author = {Lavergne, Ali{\'e}nor and Harrison, Sandy P. and Atsawawaranunt, Kamolphat and Dong, Ning and Prentice, Iain Colin},
-  year = {2022},
-  journal = {Global Ecology and Biogeography (submitted)},
-  file = {/Users/dorme/Zotero/storage/I797PRT2/Lavergneetal_GEB.docx}
-}
-
 @article{Li:2014bc,
   title = {Simulation of Tree-Ring Widths with a Model for Primary Production, Carbon Allocation, and Growth},
   author = {Li, G and Harrison, S P and Prentice, I. C. and Falster, D},
@@ -564,22 +551,6 @@ @article{Lin:2015wh
   file = {/Users/dorme/Zotero/storage/J7WANJNI/2015 Lin-2.pdf}
 }
 
-@article{linacre:1968a,
-  title = {Estimating the Net-Radiation Flux},
-  author = {Linacre, E. T.},
-  year = {1968},
-  month = jan,
-  journal = {Agricultural Meteorology},
-  volume = {5},
-  number = {1},
-  pages = {49--63},
-  issn = {0002-1571},
-  doi = {10.1016/0002-1571(68)90022-8},
-  urldate = {2024-08-02},
-  abstract = {A major influence controlling the water loss from irrigated crops is the net-radiation intensity Qn, but measurements of this are not normally available, and so attempts are often made to deduce it from other climatic data, such as the solar-radiation. Here it is shown that the relationship between net and solar-radiation intensities depends on the degree of cloudiness and the ambient temperature. By making appropriate assumptions, a series of expressions for Qn is derived, with decreasing accuracy but increasing simplicity of estimation. It appears that clouds lower the net-radiation intensity when it exceeds a critical value in the region of 0.02 cal./cm2 min, but increase it when the intensity is lower.},
-  file = {/Users/dorme/Zotero/storage/UY4NRV4W/0002157168900228.html}
-}
-
 @article{long:1993a,
   title = {Quantum Yields for Uptake of Carbon Dioxide in {{C3}} Vascular Plants of Contrasting Habitats and Taxonomic Groupings},
   author = {Long, S. P. and Postl, W. F. and {Bolh{\'a}r-Nordenkampf}, H. R.},
@@ -705,15 +676,6 @@ @article{Prentice:2014bc
   file = {/Users/dorme/Zotero/storage/MDRBDVTF/Ecol Lett 2014 Prentice.pdf}
 }
 
-@article{sandoval:in_prep,
-  title = {Aridity and Growth Temperature Effects on Phio Manuscript},
-  author = {Sandoval, David},
-  year = {in\_prep},
-  journal = {Placeholder},
-  volume = {?},
-  pages = {??-??}
-}
-
 @article{Smith:2019dv,
   title = {Global Photosynthetic Capacity Is Optimized to the Environment},
   author = {Smith, Nicholas G and Keenan, Trevor F and Colin Prentice, I and Wang, Han and Wright, Ian J. and Niinemets, {\"U}lo and Crous, Kristine Y and Domingues, Tomas F and Guerrieri, Rossella and Yoko Ishida, F and Kattge, Jens and Kruger, Eric L and Maire, Vincent and Rogers, Alistair and Serbin, Shawn P and Tarvainen, Lasse and Togashi, Henrique F and Townsend, Philip A and Wang, Meng and Weerasinghe, Lasantha K and Zhou, Shuang Xi},
@@ -734,6 +696,25 @@ @article{Smith:2019dv
   file = {/Users/dorme/Zotero/storage/EIWDBL6T/Ecol Lett 2019 Smith.pdf}
 }
 
+@article{Stocker:2020dh,
+  title = {P-Model v1.0: An Optimality-Based Light Use Efficiency Model for Simulating Ecosystem Gross Primary Production},
+  author = {Stocker, Benjamin D and Wang, Han and Smith, Nicholas G and Harrison, Sandy P and Keenan, Trevor F and Sandoval, David and Davis, Tyler and Prentice, I. Colin},
+  year = {2020},
+  journal = {Geoscientific Model Development},
+  volume = {13},
+  number = {3},
+  pages = {1545--1581},
+  doi = {10.5194/gmd-13-1545-2020},
+  date-added = {2020-11-30T12:24:06GMT},
+  date-modified = {2021-01-28T11:55:51GMT},
+  langid = {english},
+  local-url = {file://localhost/Users/dorme/References/Library.papers3/Articles/2020/Stocker/Geoscientific\%20Model\%20Development\%202020\%20Stocker.pdf},
+  rating = {0},
+  read = {Yes},
+  uri = {papers3://publication/doi/10.5194/gmd-13-1545-2020},
+  file = {/Users/dorme/Zotero/storage/J4YM5X2R/Geoscientific Model Development 2020 Stocker.pdf}
+}
+
 @article{Stocker:2018be,
   title = {Quantifying Soil Moisture Impacts on Light Use Efficiency across Biomes},
   author = {Stocker, Benjamin D and Zscheischler, Jakob and Keenan, Trevor F and Prentice, I. Colin and Penuelas, Josep and Seneviratne, Sonia I},
@@ -754,25 +735,6 @@ @article{Stocker:2018be
   file = {/Users/dorme/Zotero/storage/E6PPD38Z/New Phytol. 2018 Stocker.pdf}
 }
 
-@article{Stocker:2020dh,
-  title = {P-Model v1.0: An Optimality-Based Light Use Efficiency Model for Simulating Ecosystem Gross Primary Production},
-  author = {Stocker, Benjamin D and Wang, Han and Smith, Nicholas G and Harrison, Sandy P and Keenan, Trevor F and Sandoval, David and Davis, Tyler and Prentice, I. Colin},
-  year = {2020},
-  journal = {Geoscientific Model Development},
-  volume = {13},
-  number = {3},
-  pages = {1545--1581},
-  doi = {10.5194/gmd-13-1545-2020},
-  date-added = {2020-11-30T12:24:06GMT},
-  date-modified = {2021-01-28T11:55:51GMT},
-  langid = {english},
-  local-url = {file://localhost/Users/dorme/References/Library.papers3/Articles/2020/Stocker/Geoscientific\%20Model\%20Development\%202020\%20Stocker.pdf},
-  rating = {0},
-  read = {Yes},
-  uri = {papers3://publication/doi/10.5194/gmd-13-1545-2020},
-  file = {/Users/dorme/Zotero/storage/J4YM5X2R/Geoscientific Model Development 2020 Stocker.pdf}
-}
-
 @article{Togashi:2018es,
   title = {Functional Trait Variation Related to Gap Dynamics in Tropical Moist forests{{{\textsubscript{A}}}} Vegetation Modelling Perspective},
   author = {Togashi, Henrique F{\"u}rstenau and Atkin, Owen K and Bloomfield, Keith J and Bradford, Matt and Cao, Kunfang and Dong, Ning and Evans, Bradley J and Fan, Zexin and Harrison, Sandy P and Hua, Zhu and Liddell, Michael J and Lloyd, Jon and Ni, Jian and Wang, Han and Weerasinghe, Lasantha K and Prentice, Iain Colin},
@@ -847,25 +809,6 @@ @article{Walker:2014ce
   file = {/Users/dorme/Zotero/storage/N8LKZ4EP/Ecology and Evolution 2014 Walker.pdf}
 }
 
-@article{Wang:2017go,
-  title = {Towards a Universal Model for Carbon Dioxide Uptake by Plants},
-  author = {Wang, Han and Prentice, I. Colin and Keenan, Trevor F and Davis, Tyler W and Wright, Ian J. and Cornwell, William K and Evans, Bradley J and Peng, Changhui},
-  year = {2017},
-  month = sep,
-  journal = {Nature Plants},
-  pages = {1--8},
-  publisher = {Springer US},
-  doi = {10.1038/s41477-017-0006-8},
-  abstract = {Nature Plants, doi:10.1038/s41477-017-0006-8},
-  date-added = {2020-11-30T12:27:00GMT},
-  date-modified = {2021-01-28T09:14:19GMT},
-  local-url = {file://localhost/Users/dorme/References/Library.papers3/Articles/2017/Wang/Nature\%20Plants\%202017\%20Wang.pdf},
-  rating = {0},
-  read = {Yes},
-  uri = {papers3://publication/doi/10.1038/s41477-017-0006-8},
-  file = {/Users/dorme/Zotero/storage/D4RYI3NP/Nature Plants 2017 Wang.pdf}
-}
-
 @article{Wang:2020ik,
   title = {Acclimation of Leaf Respiration Consistent with Optimal Photosynthetic Capacity},
   author = {Wang, Han and Atkin, Owen K and Keenan, Trevor F and Smith, Nicholas G and Wright, Ian J. and Bloomfield, Keith J and Kattge, Jens and Reich, Peter B and Prentice, I. Colin},
@@ -886,6 +829,25 @@ @article{Wang:2020ik
   file = {/Users/dorme/Zotero/storage/2XYW4BF4/gcb15314-sup-0004-Supinfo.pdf;/Users/dorme/Zotero/storage/RAGU56IZ/Global Change Biol 2020 Wang.pdf}
 }
 
+@article{Wang:2017go,
+  title = {Towards a Universal Model for Carbon Dioxide Uptake by Plants},
+  author = {Wang, Han and Prentice, I. Colin and Keenan, Trevor F and Davis, Tyler W and Wright, Ian J. and Cornwell, William K and Evans, Bradley J and Peng, Changhui},
+  year = {2017},
+  month = sep,
+  journal = {Nature Plants},
+  pages = {1--8},
+  publisher = {Springer US},
+  doi = {10.1038/s41477-017-0006-8},
+  abstract = {Nature Plants, doi:10.1038/s41477-017-0006-8},
+  date-added = {2020-11-30T12:27:00GMT},
+  date-modified = {2021-01-28T09:14:19GMT},
+  local-url = {file://localhost/Users/dorme/References/Library.papers3/Articles/2017/Wang/Nature\%20Plants\%202017\%20Wang.pdf},
+  rating = {0},
+  read = {Yes},
+  uri = {papers3://publication/doi/10.1038/s41477-017-0006-8},
+  file = {/Users/dorme/Zotero/storage/D4RYI3NP/Nature Plants 2017 Wang.pdf}
+}
+
 @book{woolf:1968a,
   title = {On the {{Computation}} of {{Solar Elevation Angles}} and the {{Determination}} of {{Sunrise}} and {{Sunset Times}}},
   author = {Woolf, Harold M.},
@@ -893,3 +855,54 @@ @book{woolf:1968a
   publisher = {{National Aeronautics and Space Administration}},
   langid = {english}
 }
+
+@book{duffie:2013a,
+  title = {Solar {{Engineering}} of {{Thermal Processes}}},
+  author = {Duffie, John A. and Beckman, William A.},
+  year = {2013},
+  month = apr,
+  publisher = {John Wiley \& Sons},
+  abstract = {The updated fourth edition of the "bible" of solar energy theory and applications Over several editions, Solar Engineering of Thermal Processes has become a classic solar engineering text and reference. This revised Fourth Edition offers current coverage of solar energy theory, systems design, and applications in different market sectors along with an emphasis on solar system design and analysis using simulations to help readers translate theory into practice. An important resource for students of solar engineering, solar energy, and alternative energy as well as professionals working in the power and energy industry or related fields, Solar Engineering of Thermal Processes, Fourth Edition features:  Increased coverage of leading-edge topics such as photovoltaics and the design of solar cells and heaters A brand-new chapter on applying CombiSys (a readymade TRNSYS simulation program available for free download) to simulate a solar heated house with solar- heated domestic hot water Additional simulation problems available through a companion website An extensive array of homework problems and exercises},
+  googlebooks = {5uDdUfMgXYQC},
+  isbn = {978-1-118-41541-2},
+  langid = {english},
+  keywords = {Technology & Engineering / Electronics / General,Technology & Engineering / Mechanical,Technology & Engineering / Power Resources / General}
+}
+
+@article{linacre:1968a,
+  title = {Estimating the Net-Radiation Flux},
+  author = {Linacre, E. T.},
+  year = {1968},
+  month = jan,
+  journal = {Agricultural Meteorology},
+  volume = {5},
+  number = {1},
+  pages = {49--63},
+  issn = {0002-1571},
+  doi = {10.1016/0002-1571(68)90022-8},
+  urldate = {2024-08-02},
+  abstract = {A major influence controlling the water loss from irrigated crops is the net-radiation intensity Qn, but measurements of this are not normally available, and so attempts are often made to deduce it from other climatic data, such as the solar-radiation. Here it is shown that the relationship between net and solar-radiation intensities depends on the degree of cloudiness and the ambient temperature. By making appropriate assumptions, a series of expressions for Qn is derived, with decreasing accuracy but increasing simplicity of estimation. It appears that clouds lower the net-radiation intensity when it exceeds a critical value in the region of 0.02 cal./cm2 min, but increase it when the intensity is lower.},
+  file = {/Users/dorme/Zotero/storage/UY4NRV4W/0002157168900228.html}
+}
+
+@article{sandoval:in_prep,
+  title = {Aridity and Growth Temperature Effects on Phio Manuscript},
+  author = {Sandoval, David},
+  year = {in\_prep},
+  journal = {Placeholder},
+  volume = {?},
+  pages = {??-??}
+}
+
+@techreport{joshi:2022a,
+  title = {Plant-{{FATE}}: {{Predicting}} the Adaptive Responses of Biodiverse Plant Communities Using Functional-Trait Evolution},
+  author = {Joshi, Jaideep and Prentice, Iain Colin and Br{\"a}nnstr{\"o}m, {\AA}ke and Singh, Shipra and Hofhansl, Florian and Dieckmann, Ulf},
+  year = {2022},
+  month = mar,
+  number = {EGU22-9994},
+  institution = {Copernicus Meetings},
+  doi = {10.5194/egusphere-egu22-9994},
+  urldate = {2024-09-05},
+  langid = {english},
+  file = {/Users/dorme/Zotero/storage/47RVAZGK/EGU22-9994.html}
+}
diff --git a/docs/source/users/tmodel/canopy.md b/docs/source/users/tmodel/canopy.md
new file mode 100644
index 00000000..b0911b87
--- /dev/null
+++ b/docs/source/users/tmodel/canopy.md
@@ -0,0 +1,677 @@
+---
+jupytext:
+  formats: md:myst
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 0.13
+kernelspec:
+  display_name: python3
+  language: python
+  name: python3
+---
+
+# Canopy model
+
+This notebook walks through the steps in generating the canopy model as (hopefully) used
+in Plant-FATE.
+
+## The T Model
+
+The T Model provides a numerical description of how tree geometry scales with stem
+diameter, and an allocation model of how GPP predicts changes in stem diameter.
+
+The implementation in `pyrealm` provides a class representing a particular plant
+functional type, using a set of traits. The code below creates a PFT with the default
+set of trait values.
+
+```{warning}
+This sketch:
+
+* Assumes a single individual of each stem diameter, but in practice
+  we are going to want to include a number of individuals to capture cohorts.
+* Assumes a single PFT, where we will need to provide a mixed community.
+* Consequently handles forest inventory properties in a muddled way: we will
+  likely package all of the stem data into a single class, probably a community
+  object.
+```
+
+```{code-cell}
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.optimize import root_scalar
+
+from pyrealm.tmodel import TTree
+
+np.set_printoptions(precision=3)
+
+# Plant functional type with default parameterization
+pft = TTree(diameters=np.array([0.1, 0.15, 0.2, 0.25, 0.38, 0.4, 1.0]))
+pft.traits
+```
+
+The scaling of a set of trees is automatically calculated using the initial diameters to
+the `TTree` instance. This automatically calculates the other dimensions, such as
+height, using the underlying scaling equations of the T Model.
+
+```{code-cell}
+pft.height
+```
+
+```{code-cell}
+:lines_to_next_cell: 2
+
+pft.crown_area
+```
+
++++ {"lines_to_next_cell": 2}
+
+### Crown shape
+
+Jaideep's extension of the T Model adds a crown shape model, driven by two parameters
+($m$ and $n$) that provide a very flexible description of the vertical crown profile.
+These are used to derive two further invariant parameters: $q_m$ scales the canopy
+radius to match the predictions of crown area under the T model and $p_{zm}$ is the
+proportion of the total stem height at which the maximum crown radius occurs. Both these
+values are constant for a plant functional type.
+
+$$
+\begin{align}
+q_m &= m n \left(\dfrac{n-1}{m n -1}\right)^ {1 - \tfrac{1}{n}}
+    \left(\dfrac{\left(m-1\right) n}{m n -1}\right)^ {m-1} \\[8pt]
+p_{zm} &= \left(\dfrac{n-1}{m n -1}\right)^ {\tfrac{1}{n}}
+\end{align}
+$$
+
+For individual stems, with expected height $H$ and crown area $A_c$, we can then
+estimate:
+
+* the height $z_m$ at which the maximum crown radius $r_m$ is found, and
+* a slope $r_0$ that scales the relative canopy radius so that the $r_m$ matches
+  the allometric prediction of $A_c$ from the T Model.
+
+$$
+\begin{align}
+z_m &= H p_{zm}\\[8pt]
+r_0 &= \frac{1}{q_m}\sqrt{\frac{A_c}{\pi}}
+\end{align}
+$$
+
+```{code-cell}
+:lines_to_next_cell: 2
+
+def calculate_qm(m, n):
+
+    # Constant q_m
+    return (
+        m
+        * n
+        * ((n - 1) / (m * n - 1)) ** (1 - 1 / n)
+        * (((m - 1) * n) / (m * n - 1)) ** (m - 1)
+    )
+
+
+def calculate_stem_canopy_factors(pft, m, n):
+
+    # Height of maximum crown radius
+    zm = pft.height * ((n - 1) / (m * n - 1)) ** (1 / n)
+
+    # Slope to give Ac at zm
+    r0 = 1 / qm * np.sqrt(pft.crown_area / np.pi)
+
+    return zm, r0
+
+
+# Shape parameters for a fairly top heavy crown profile
+m = 2
+n = 5
+qm = calculate_qm(m=m, n=n)
+zm, r0 = calculate_stem_canopy_factors(pft=pft, m=m, n=n)
+
+print("qm = ", np.round(qm, 4))
+print("zm = ", zm)
+print("r0 = ", r0)
+```
+
++++ {"lines_to_next_cell": 2}
+
+The following functions then provide the value at height $z$ of relative $q(z)$ and
+actual $r(z)$ canopy radius:
+
+$$
+\begin{align}
+q(z) &= m n \left(\dfrac{z}{H}\right) ^ {n -1}
+    \left( 1 - \left(\dfrac{z}{H}\right) ^ n \right)^{m-1}\\[8pt]
+r(z) &= r_0 \; q(z)
+\end{align}
+$$
+
+```{code-cell}
+def calculate_relative_canopy_radius_at_z(z, H, m, n):
+    """Calculate q(z)"""
+
+    z_over_H = z / H
+
+    return m * n * z_over_H ** (n - 1) * (1 - z_over_H**n) ** (m - 1)
+```
+
+```{code-cell}
+# Validate that zm and r0 generate the predicted maximum crown area
+q_zm = calculate_relative_canopy_radius_at_z(zm, pft.height, m, n)
+rm = r0 * q_zm
+print("rm = ", rm)
+```
+
+```{code-cell}
+np.allclose(rm**2 * np.pi, pft.crown_area)
+```
+
+Vertical crown radius profiles can now be calculated for each stem:
+
+```{code-cell}
+# Create an interpolation from ground to maximum stem height, with 5 cm resolution.
+# Also append a set of values _fractionally_ less than the exact height  of stems
+# so that the height at the top of each stem is included but to avoid floating
+# point issues with exact heights.
+
+zres = 0.05
+z = np.arange(0, pft.height.max() + 1, zres)
+z = np.sort(np.concatenate([z, pft.height - 0.00001]))
+
+# Convert the heights into a column matrix to broadcast against the stems
+# and then calculate r(z) = r0 * q(z)
+rz = r0 * calculate_relative_canopy_radius_at_z(z[:, None], pft.height, m, n)
+
+# When z > H, rz < 0, so set radius to 0 where rz < 0
+rz[np.where(rz < 0)] = 0
+
+np.cumsum(np.convolve(rm, np.ones(2), "valid") + 0.1)
+```
+
+Those can be plotted out to show the vertical crown radius profiles
+
+```{code-cell}
+# Separate the stems along the x axis for plotting
+stem_x = np.concatenate(
+    [np.array([0]), np.cumsum(np.convolve(rm, np.ones(2), "valid") + 0.4)]
+)
+
+# Get the canopy sections above and below zm
+rz_above_zm = np.where(np.logical_and(rz > 0, np.greater.outer(z, zm)), rz, np.nan)
+rz_below_zm = np.where(np.logical_and(rz > 0, np.less_equal.outer(z, zm)), rz, np.nan)
+
+# Plot the canopy parts
+plt.plot(stem_x + rz_below_zm, z, color="khaki")
+plt.plot(stem_x - rz_below_zm, z, color="khaki")
+plt.plot(stem_x + rz_above_zm, z, color="forestgreen")
+plt.plot(stem_x - rz_above_zm, z, color="forestgreen")
+
+# Add the maximum radius
+plt.plot(np.vstack((stem_x - rm, stem_x + rm)), np.vstack((zm, zm)), color="firebrick")
+
+# Plot the stem centre lines
+plt.vlines(stem_x, 0, pft.height, linestyles="-", color="grey")
+
+plt.gca().set_aspect("equal")
+```
+
+## Canopy structure
+
+The canopy structure model uses the perfect plasticity approximation (PPA), which
+assumes that plants can arrange their canopies to fill the available space $A$.
+It takes the **projected area of stems** $Ap(z)$ within the canopy and finds the heights
+at which each canopy layer closes ($z^*_l$ for $l = 1, 2, 3 ...$) where the total projected
+area of the canopy equals $lA$.
+
+### Canopy projected area
+
+The projected area $A_p$ for a stem with height $H$, a maximum crown area $A_c$ at a
+height $z_m$ and $m$, $n$ and $q_m$ for the associated plant functional type is
+
+$$
+A_p(z)=
+\begin{cases}
+A_c, & z \le z_m \\
+A_c \left(\dfrac{q(z)}{q_m}\right)^2, & H > z > z_m \\
+0, & z > H
+\end{cases}
+$$
+
+```{code-cell}
+Stems = float | np.ndarray
+
+
+def calculate_projected_area(
+    z: float,
+    pft,
+    m: Stems,
+    n: Stems,
+    qm: Stems,
+    zm: Stems,
+) -> np.ndarray:
+    """Calculate projected crown area above a given height.
+
+    This function takes PFT specific parameters (shape parameters) and stem specific
+    sizes and estimates the projected crown area above a given height $z$. The inputs
+    can either be scalars describing a single stem or arrays representing a community
+    of stems. If only a single PFT is being modelled then `m`, `n`, `qm` and `fg` can
+    be scalars with arrays `H`, `Ac` and `zm` giving the sizes of stems within that
+    PFT.
+
+    Args:
+        z: Canopy height
+        m, n, qm : PFT specific shape parameters
+        pft, qm, zm: stem data
+
+    """
+
+    # Calculate q(z)
+    qz = calculate_relative_canopy_radius_at_z(z, pft.height, m, n)
+
+    # Calculate Ap given z > zm
+    Ap = pft.crown_area * (qz / qm) ** 2
+    # Set Ap = Ac where z <= zm
+    Ap = np.where(z <= zm, pft.crown_area, Ap)
+    # Set Ap = 0 where z > H
+    Ap = np.where(z > pft.height, 0, Ap)
+
+    return Ap
+```
+
+The code below calculates the projected crown area for each stem and then plots the
+vertical profile for individual stems and across the community.
+
+```{code-cell}
+:lines_to_next_cell: 2
+
+# Calculate the projected area for each stem
+Ap_z = calculate_projected_area(z=z[:, None], pft=pft, m=m, n=n, qm=qm, zm=zm)
+
+# Plot the calculated values for each stem and across the community
+fig, (ax1, ax2) = plt.subplots(1, 2, sharey=True, figsize=(10, 5))
+
+# Plot the individual stem projected area
+ax1.set_ylabel("Height above ground ($z$, m)")
+ax1.plot(Ap_z, z)
+ax1.set_xlabel("Stem $A_p(z)$ (m2)")
+
+# Plot the individual stem projected area
+ax2.plot(np.nansum(Ap_z, axis=1), z)
+ax2.set_xlabel("Total community $A_p(z)$ (m2)")
+
+plt.tight_layout()
+```
+
++++ {"lines_to_next_cell": 2}
+
+### Canopy closure and canopy gap fraction
+
+The total cumulative projected area shown above is modified by a community-level
+**canopy gap fraction** ($f_G$) that captures the overall proportion of the canopy area
+that is left unfilled by canopy. This gap fraction, capturing processes such as crown
+shyness, describes the proportion of open sky visible from the forest floor.
+
+The definition of the height of canopy layer closure ($z^*_l$) for a given canopy
+layer $l = 1, ..., l_m$ is then:
+
+$$
+\sum_1^{N_s}{ A_p(z^*_l)} = l A(1 - f_G)
+$$
+
+This can be found numerically using a root solver as:
+
+$$
+\sum_1^{N_s}{ A_p(z^*_l)} - l A(1 - f_G) = 0
+$$
+
+The total number of layers $l_m$ in a canopy, where the final layer may not be fully closed,
+can be found given the total crown area across stems as:
+
+$$
+l_m = \left\lceil \frac{\sum_1^{N_s}{ A_c}}{ A(1 - f_G)}\right\rceil
+$$
+
+```{code-cell}
+def solve_canopy_closure_height(
+    z: float,
+    l: int,
+    A: float,
+    fG: float,
+    m: Stems,
+    n: Stems,
+    qm: Stems,
+    pft: Stems,
+    zm: Stems,
+) -> np.ndarray:
+    """Solver function for canopy closure height.
+
+    This function returns the difference between the total community projected area
+    at a height $z$ and the total available canopy space for canopy layer $l$, given
+    the community gap fraction for a given height. It is used with a root solver to
+    find canopy layer closure heights $z^*_l* for a community.
+
+    Args:
+        m, n, qm : PFT specific shape parameters
+        H, Ac, zm: stem specific sizes
+        A, l: cell area and layer index
+        fG: community gap fraction
+    """
+
+    # Calculate Ap(z)
+    Ap_z = calculate_projected_area(z=z, pft=pft, m=m, n=n, qm=qm, zm=zm)
+
+    # Return the difference between the projected area and the available space
+    return Ap_z.sum() - (A * l) * (1 - fG)
+
+
+def calculate_canopy_heights(
+    A: float,
+    fG: float,
+    m: Stems,
+    n: Stems,
+    qm: Stems,
+    pft,
+    zm: Stems,
+):
+
+    # Calculate the number of layers
+    total_community_ca = pft.crown_area.sum()
+    n_layers = int(np.ceil(total_community_ca / (A * (1 - fG))))
+
+    # Data store for z*
+    z_star = np.zeros(n_layers)
+
+    # Loop over the layers TODO - edge case of completely filled final layer
+    for lyr in np.arange(n_layers - 1):
+        z_star[lyr] = root_scalar(
+            solve_canopy_closure_height,
+            args=(lyr + 1, A, fG, m, n, qm, pft, zm),
+            bracket=(0, pft.height.max()),
+        ).root
+
+    return z_star
+```
+
+The example below calculates the projected crown area above ground level for the example
+stems. These should be identical to the crown area of the stems.
+
+```{code-cell}
+# Set the total available canopy space and community gap fraction
+canopy_area = 32
+community_gap_fraction = 2 / 32
+
+z_star = calculate_canopy_heights(
+    A=canopy_area, fG=community_gap_fraction, m=m, n=n, qm=qm, pft=pft, zm=zm
+)
+
+print("z_star = ", z_star)
+```
+
+We can now plot the canopy stems alongside the community $A_p(z)$ profile, and
+superimpose the calculated $z^*_l$ values and the cumulative canopy area for each layer
+to confirm that the calculated values coincide with the profile. Note here that the
+total area at each closed layer height is omitting the community gap fraction.
+
+```{code-cell}
+community_Ap_z = np.nansum(Ap_z, axis=1)
+
+fig, (ax1, ax2) = plt.subplots(1, 2, sharey=True, figsize=(10, 5))
+
+zcol = "red"
+
+# Plot the canopy parts
+ax1.plot(stem_x + rz_below_zm, z, color="khaki")
+ax1.plot(stem_x - rz_below_zm, z, color="khaki")
+ax1.plot(stem_x + rz_above_zm, z, color="forestgreen")
+ax1.plot(stem_x - rz_above_zm, z, color="forestgreen")
+
+# Add the maximum radius
+ax1.plot(np.vstack((stem_x - rm, stem_x + rm)), np.vstack((zm, zm)), color="firebrick")
+
+# Plot the stem centre lines
+ax1.vlines(stem_x, 0, pft.height, linestyles="-", color="grey")
+
+ax1.set_ylabel("Height above ground ($z$, m)")
+ax1.set_xlabel("Arbitrary stem position")
+ax1.hlines(z_star, 0, (stem_x + rm).max(), color=zcol, linewidth=0.5)
+
+# Plot the projected community crown area by height along with the heights
+# at which different canopy layers close
+ax2.hlines(z_star, 0, community_Ap_z.max(), color=zcol, linewidth=0.5)
+ax2.vlines(
+    canopy_area * np.arange(1, len(z_star)) * (1 - community_gap_fraction),
+    0,
+    pft.height.max(),
+    color=zcol,
+    linewidth=0.5,
+)
+
+ax2.plot(community_Ap_z, z)
+ax2.set_xlabel("Projected crown area above height $z$ ($A_p(z)$, m2)")
+
+# Add z* values on the righthand axis
+ax3 = ax2.twinx()
+
+
+def z_star_labels(X):
+    return [f"$z^*_{l + 1}$ = {z:.2f}" for l, z in enumerate(X)]
+
+
+ax3.set_ylim(ax2.get_ylim())
+ax3.set_yticks(z_star)
+ax3.set_yticklabels(z_star_labels(z_star))
+
+ax4 = ax2.twiny()
+
+# Add canopy layer closure areas on top axis
+cum_area = np.arange(1, len(z_star)) * canopy_area * (1 - community_gap_fraction)
+
+
+def cum_area_labels(X):
+    return [f"$A_{l + 1}$ = {z:.1f}" for l, z in enumerate(X)]
+
+
+ax4.set_xlim(ax2.get_xlim())
+ax4.set_xticks(cum_area)
+ax4.set_xticklabels(cum_area_labels(cum_area))
+
+plt.tight_layout()
+```
+
+The projected area from individual stems to each canopy layer can then be calculated at
+$z^*_l$ and hence the projected area of canopy **within each layer**.
+
+```{code-cell}
+# Calculate the canopy area above z_star for each stem
+Ap_z_star = calculate_projected_area(z=z_star[:, None], pft=pft, m=m, n=n, qm=qm, zm=zm)
+
+print(Ap_z_star)
+```
+
+```{code-cell}
+:lines_to_next_cell: 2
+
+# Calculate the contribution _within_ each layer per stem
+Ap_within_layer = np.diff(Ap_z_star, axis=0, prepend=0)
+
+print(Ap_within_layer)
+```
+
++++ {"lines_to_next_cell": 2}
+
+### Leaf area within canopy layers
+
+The projected area occupied by leaves at a given height $\tilde{A}_{cp}(z)$ is
+needed to calculate light transmission through those layers. This differs from the
+projected area $A_p(z)$ because, although a tree occupies an area in the canopy
+following the PPA, a **crown gap fraction** ($f_g$) reduces the actual leaf area
+at a given height $z$.
+
+The crown gap fraction does not affect the overall projected canopy area at ground
+level or the community gap fraction: the amount of clear sky at ground level is
+governed purely by $f_G$. Instead it models how leaf gaps in the upper canopy are
+filled by leaf area at lower heights. It captures the vertical distribution of
+leaf area within the canopy: a higher $f_g$ will give fewer leaves at the top of
+the canopy and more leaves further down within the canopy.
+
+The calculation of $\tilde{A}_{cp}(z)$ is defined as:
+
+$$
+\tilde{A}_{cp}(z)=
+\begin{cases}
+0, & z \gt H \\
+A_c \left(\dfrac{q(z)}{q_m}\right)^2 \left(1 - f_g\right), & H \gt z \gt z_m \\
+Ac - A_c \left(\dfrac{q(z)}{q_m}\right)^2 f_g, & zm \gt z
+\end{cases}
+$$
+
+The function below calculates $\tilde{A}_{cp}(z)$.
+
+```{code-cell}
+def calculate_leaf_area(
+    z: float,
+    fg: float,
+    pft,
+    m: Stems,
+    n: Stems,
+    qm: Stems,
+    zm: Stems,
+) -> np.ndarray:
+    """Calculate leaf area above a given height.
+
+    This function takes PFT specific parameters (shape parameters) and stem specific
+    sizes and estimates the projected crown area above a given height $z$. The inputs
+    can either be scalars describing a single stem or arrays representing a community
+    of stems. If only a single PFT is being modelled then `m`, `n`, `qm` and `fg` can
+    be scalars with arrays `H`, `Ac` and `zm` giving the sizes of stems within that
+    PFT.
+
+    Args:
+        z: Canopy height
+        fg: crown gap fraction
+        m, n, qm : PFT specific shape parameters
+        pft, qm, zm: stem data
+    """
+
+    # Calculate q(z)
+    qz = calculate_relative_canopy_radius_at_z(z, pft.height, m, n)
+
+    # Calculate Ac term
+    Ac_term = pft.crown_area * (qz / qm) ** 2
+    # Set Acp either side of zm
+    Acp = np.where(z <= zm, pft.crown_area - Ac_term * fg, Ac_term * (1 - fg))
+    # Set Ap = 0 where z > H
+    Acp = np.where(z > pft.height, 0, Acp)
+
+    return Acp
+```
+
+The plot below shows how the vertical leaf area profile for the community changes for
+different values of $f_g$. When $f_g = 0$, then $A_cp(z) = A_p(z)$ (red line) because
+there are no crown gaps and hence all of the leaf area is within the crown surface. As
+$f_g \to 1$, more of the leaf area is displaced deeper into the canopy, leaves in the
+lower crown intercepting light coming through holes in the upper canopy.
+
+```{code-cell}
+fig, ax1 = plt.subplots(1, 1, figsize=(6, 5))
+
+for fg in np.arange(0, 1.01, 0.05):
+
+    if fg == 0:
+        color = "red"
+        label = "$f_g = 0$"
+        lwd = 0.5
+    elif fg == 1:
+        color = "blue"
+        label = "$f_g = 1$"
+        lwd = 0.5
+    else:
+        color = "black"
+        label = None
+        lwd = 0.25
+
+    Acp_z = calculate_leaf_area(z=z[:, None], fg=fg, pft=pft, m=m, n=n, qm=qm, zm=zm)
+    ax1.plot(np.nansum(Acp_z, axis=1), z, color=color, linewidth=lwd, label=label)
+
+ax1.set_xlabel(r"Projected leaf area above height $z$ ($\tilde{A}_{cp}(z)$, m2)")
+ax1.legend(frameon=False)
+```
+
+We can now calculate the crown area occupied by leaves above the height of each closed
+layer $z^*_l$:
+
+```{code-cell}
+# Calculate the leaf area above z_star for each stem
+crown_gap_fraction = 0.05
+Acp_z_star = calculate_leaf_area(
+    z=z_star[:, None], fg=crown_gap_fraction, pft=pft, m=m, n=n, qm=qm, zm=zm
+)
+
+print(Acp_z_star)
+```
+
+And from that, the area occupied by leaves **within each layer**. These values are
+similar to the projected crown area within layers (`Ap_within_layer`, above) but
+leaf area is displaced into lower layers because $f_g > 0$.
+
+```{code-cell}
+# Calculate the contribution _within_ each layer per stem
+Acp_within_layer = np.diff(Acp_z_star, axis=0, prepend=0)
+
+print(Acp_within_layer)
+```
+
+### Light transmission through the canopy
+
+Now we can use the leaf area by layer and the Beer-Lambert equation to calculate light
+attenuation through the canopy layers.
+
+$f_{abs} = 1 - e ^ {-kL}$,
+
+where $k$ is the light extinction coefficient ($k$) and $L$ is the leaf area index
+(LAI). The LAI can be calculated for each stem and layer:
+
+```{code-cell}
+LAI = Acp_within_layer / canopy_area
+print(LAI)
+```
+
+This can be used to calculate the LAI of individual stems but also the LAI of each layer
+in the canopy:
+
+```{code-cell}
+LAI_stem = LAI.sum(axis=0)
+LAI_layer = LAI.sum(axis=1)
+
+print("LAI stem = ", LAI_stem)
+print("LAI layer = ", LAI_layer)
+```
+
+The layer LAI values can now be used to calculate the light transmission of each layer and
+hence the cumulative light extinction profile through the canopy.
+
+```{code-cell}
+f_abs = 1 - np.exp(-pft.traits.par_ext * LAI_layer)
+ext = np.cumprod(f_abs)
+
+print("f_abs = ", f_abs)
+print("extinction = ", ext)
+```
+
+One issue that needs to be resolved is that the T Model implementation in `pyrealm`
+follows the original implementation of the T Model in having LAI as a fixed trait of
+a given plant functional type, so is constant for all stems of that PFT.
+
+```{code-cell}
+print("f_abs = ", (1 - np.exp(-pft.traits.par_ext * pft.traits.lai)))
+```
+
+## Things to worry about later
+
+Herbivory - leaf fall (transmission increases, truncate at 0, replace from NSCs) vs leaf
+turnover (transmission constant, increased GPP penalty)
+
+Leaf area dynamics in PlantFATE - acclimation to herbivory and transitory decreases in
+transimission, need non-structural carbohydrates  to recover from total defoliation.
+
+Leaf economics.
diff --git a/poetry.lock b/poetry.lock
index f21ca1e6..aca832f9 100644
--- a/poetry.lock
+++ b/poetry.lock
@@ -1794,6 +1794,48 @@ files = [
     {file = "MarkupSafe-2.1.5.tar.gz", hash = "sha256:d283d37a890ba4c1ae73ffadf8046435c76e7bc2247bbb63c00bd1a709c6544b"},
 ]
 
+[[package]]
+name = "marshmallow"
+version = "3.22.0"
+description = "A lightweight library for converting complex datatypes to and from native Python datatypes."
+optional = false
+python-versions = ">=3.8"
+files = [
+    {file = "marshmallow-3.22.0-py3-none-any.whl", hash = "sha256:71a2dce49ef901c3f97ed296ae5051135fd3febd2bf43afe0ae9a82143a494d9"},
+    {file = "marshmallow-3.22.0.tar.gz", hash = "sha256:4972f529104a220bb8637d595aa4c9762afbe7f7a77d82dc58c1615d70c5823e"},
+]
+
+[package.dependencies]
+packaging = ">=17.0"
+
+[package.extras]
+dev = ["marshmallow[tests]", "pre-commit (>=3.5,<4.0)", "tox"]
+docs = ["alabaster (==1.0.0)", "autodocsumm (==0.2.13)", "sphinx (==8.0.2)", "sphinx-issues (==4.1.0)", "sphinx-version-warning (==1.1.2)"]
+tests = ["pytest", "pytz", "simplejson"]
+
+[[package]]
+name = "marshmallow-dataclass"
+version = "8.7.0"
+description = "Python library to convert dataclasses into marshmallow schemas."
+optional = false
+python-versions = ">=3.8"
+files = [
+    {file = "marshmallow_dataclass-8.7.0-py3-none-any.whl", hash = "sha256:9e528d72b83f2b6b0f60cb29fd38781a6f7ce2155295adb1ed33289826a93c4b"},
+    {file = "marshmallow_dataclass-8.7.0.tar.gz", hash = "sha256:0218008fec3fd4b5f739b2a0c6d7593bcc403308f6da953e341e4e359e268849"},
+]
+
+[package.dependencies]
+marshmallow = ">=3.18.0"
+typeguard = ">=4.0.0,<4.1.0"
+typing-extensions = {version = ">=4.2.0", markers = "python_version < \"3.11\""}
+typing-inspect = ">=0.9.0,<0.10.0"
+
+[package.extras]
+dev = ["pre-commit (>=2.17,<3.0)", "pytest (>=5.4)", "pytest-mypy-plugins (>=1.2.0)", "sphinx"]
+docs = ["sphinx"]
+lint = ["pre-commit (>=2.17,<3.0)"]
+tests = ["pytest (>=5.4)", "pytest-mypy-plugins (>=1.2.0)"]
+
 [[package]]
 name = "matplotlib"
 version = "3.9.1"
@@ -3758,6 +3800,24 @@ files = [
 docs = ["myst-parser", "pydata-sphinx-theme", "sphinx"]
 test = ["argcomplete (>=3.0.3)", "mypy (>=1.7.0)", "pre-commit", "pytest (>=7.0,<8.2)", "pytest-mock", "pytest-mypy-testing"]
 
+[[package]]
+name = "typeguard"
+version = "4.0.1"
+description = "Run-time type checker for Python"
+optional = false
+python-versions = ">=3.7.4"
+files = [
+    {file = "typeguard-4.0.1-py3-none-any.whl", hash = "sha256:43f55cc9953f26dae362adb973b6c9ad6b97bfffcc6757277912eddd5cfa345b"},
+    {file = "typeguard-4.0.1.tar.gz", hash = "sha256:db35142d1f92fc8c1b954e5cc03b57810428f9cd4e4604647bdf5764fc5bbba9"},
+]
+
+[package.dependencies]
+typing-extensions = {version = ">=4.7.0", markers = "python_version < \"3.12\""}
+
+[package.extras]
+doc = ["Sphinx (<7)", "packaging", "sphinx-autodoc-typehints (>=1.2.0)", "sphinx-rtd-theme"]
+test = ["mypy (>=1.2.0)", "pytest (>=7)"]
+
 [[package]]
 name = "types-python-dateutil"
 version = "2.9.0.20240316"
@@ -3802,6 +3862,21 @@ files = [
     {file = "typing_extensions-4.12.2.tar.gz", hash = "sha256:1a7ead55c7e559dd4dee8856e3a88b41225abfe1ce8df57b7c13915fe121ffb8"},
 ]
 
+[[package]]
+name = "typing-inspect"
+version = "0.9.0"
+description = "Runtime inspection utilities for typing module."
+optional = false
+python-versions = "*"
+files = [
+    {file = "typing_inspect-0.9.0-py3-none-any.whl", hash = "sha256:9ee6fc59062311ef8547596ab6b955e1b8aa46242d854bfc78f4f6b0eff35f9f"},
+    {file = "typing_inspect-0.9.0.tar.gz", hash = "sha256:b23fc42ff6f6ef6954e4852c1fb512cdd18dbea03134f91f856a95ccc9461f78"},
+]
+
+[package.dependencies]
+mypy-extensions = ">=0.3.0"
+typing-extensions = ">=3.7.4"
+
 [[package]]
 name = "tzdata"
 version = "2024.1"
@@ -3959,4 +4034,4 @@ test = ["big-O", "importlib-resources", "jaraco.functools", "jaraco.itertools",
 [metadata]
 lock-version = "2.0"
 python-versions = ">=3.10"
-content-hash = "47bf9006da7dea679c3ebd0892d08041d477bbcc3d898326280382a6e48777c8"
+content-hash = "aef31655540ecd90f7e317fe133563bcfa68b2dadbd711cb61389563fd57e555"
diff --git a/pyproject.toml b/pyproject.toml
index 7110af96..a170c766 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -37,6 +37,9 @@ python = ">=3.10"
 scipy = "^1.7.3" 
 tabulate = "^0.8.10" 
 
+marshmallow = "^3.22.0"
+pandas = "^2.2.2"
+marshmallow-dataclass = "^8.7.0"
 [tool.poetry.group.types.dependencies]
 pandas-stubs = "^2.2.0.240218"
 types-tabulate = "^0.9.0.0"
diff --git a/pyrealm/demography/flora.py b/pyrealm/demography/flora.py
new file mode 100644
index 00000000..6c754a7d
--- /dev/null
+++ b/pyrealm/demography/flora.py
@@ -0,0 +1,312 @@
+"""The flora module implements definitions of:
+
+* The ``PlantFunctionalType`` and ``PlantFunctionalTypeStrict`` dataclasses, which are
+  used to parameterise the traits of different plant functional types. The
+  ``PlantFunctionalType`` dataclass is a subclass of ``PlantFunctionalTypeStrict`` that
+  simply adds default values to the attributes.
+* The ``PlantFunctionalTypeStrict`` dataclass is used as the basis of a ``marshmallow``
+  schema for validating the creation of plant functional types from data files. This
+  intentionally enforces a complete description of the traits in the input data. The
+  ``PlantFunctionalType`` is provided as a more convenient API for programmatic use.
+* The Flora class, which is simply a dictionary of named plant functional types for use
+  in describing a plant community in a simulation. The Flora class also defines factory
+  methods to create instances from plant functional type data stored in JSON, TOML or
+  CSV formats.
+"""  # noqa: D415
+
+from __future__ import annotations
+
+import json
+import sys
+from collections import Counter
+from collections.abc import Sequence
+from dataclasses import dataclass, field
+from pathlib import Path
+
+import marshmallow_dataclass
+import pandas as pd
+from marshmallow.exceptions import ValidationError
+
+if sys.version_info[:2] >= (3, 11):
+    import tomllib
+    from tomllib import TOMLDecodeError
+else:
+    import tomli as tomllib
+    from tomli import TOMLDecodeError
+
+
+@dataclass(frozen=True)
+class PlantFunctionalTypeStrict:
+    """The PlantFunctionalTypeStrict dataclass.
+
+    This dataclass implements the set of traits required to define a plant functional
+    type for use in ``pyrealm``.
+
+    * Most traits are taken from the definition of the T Model of plant growth and GPP
+      allocation :cite:`Li:2014bc`.
+    * The foliage maintenance respiration fraction was not explicitly included in
+      :cite:t:`Li:2014bc` - there was assumed to be a 10% penalty on GPP before
+      calculating the other component - but has been explicitly included here.
+    * This implementation adds two further canopy shape parameters (``m`` and ``n``),
+      which are then used to calculate two derived attributes (``q_m`` and
+      ``z_max_ratio``). These are used to define the vertical distribution of leaves
+      around a stem and follow the implementation developed in the PlantFATE model
+      :cite:`joshi:2022a`.
+
+    See also :class:`~pyrealm.demography.flora.PlantFunctionalType` for the default
+    values implemented in that subclass.
+    """
+
+    name: str
+    r"""The name of the plant functional type."""
+    a_hd: float
+    r"""Initial slope of height-diameter relationship (:math:`a`, -)"""
+    ca_ratio: float
+    r"""Initial ratio of crown area to stem cross-sectional area
+    (:math:`c`, -)"""
+    h_max: float
+    r"""Maximum tree height (:math:`H_m`, m)"""
+    rho_s: float
+    r"""Sapwood density (:math:`\rho_s`, kg Cm-3)"""
+    lai: float
+    """Leaf area index within the crown (:math:`L`,  -)"""
+    sla: float
+    r"""Specific leaf area (:math:`\sigma`,  m2 kg-1 C)"""
+    tau_f: float
+    r"""Foliage turnover time (:math:`\tau_f`,years)"""
+    tau_r: float
+    r"""Fine-root turnover time (:math:`\tau_r`,  years)"""
+    par_ext: float
+    r"""Extinction coefficient of photosynthetically active radiation (PAR) (:math:`k`,
+     -)"""
+    yld: float
+    r"""Yield factor (:math:`y`,  -)"""
+    zeta: float
+    r"""Ratio of fine-root mass to foliage area (:math:`\zeta`, kg C m-2)"""
+    resp_r: float
+    r"""Fine-root specific respiration rate (:math:`r_r`, year-1)"""
+    resp_s: float
+    r"""Sapwood-specific respiration rate (:math:`r_s`,  year-1)"""
+    resp_f: float
+    r"""Foliage maintenance respiration fraction (:math:`r_f`,  -)"""
+    m: float
+    r"""Canopy shape parameter (:math:`m`, -)"""
+    n: float
+    r"""Canopy shape parameter (:math:`n`, -)"""
+
+    q_m: float = field(init=False)
+    """Scaling factor to derive maximum crown radius from crown area."""
+    z_max_prop: float = field(init=False)
+    """Proportion of stem height at which maximum crown radius is found."""
+
+    def __post_init__(self) -> None:
+        """Populate derived attributes.
+
+        This method populates the ``q_m`` and ``z_max_ratio`` attributes from the
+        provided values of ``m`` and ``n``.
+        """
+
+        # Calculate q_m and z_max proportion. Need to use __setattr__ because the
+        # dataclass is frozen.
+        object.__setattr__(self, "q_m", calculate_q_m(m=self.m, n=self.n))
+        object.__setattr__(
+            self, "z_max_prop", calculate_z_max_proportion(m=self.m, n=self.n)
+        )
+
+
+@dataclass(frozen=True)
+class PlantFunctionalType(PlantFunctionalTypeStrict):
+    r"""The PlantFunctionalType dataclass.
+
+    This dataclass is a subclass of
+    :class:`~pyrealm.demography.flora.PlantFunctionalTypeStrict` that implements exactly
+    the same set of traits but provides default values. This class is intended as a
+    convenience API for programmatic use, where the parent provides a strict schema for
+    generating plant functional type instances from data.
+
+    The table below lists the attributes and default values taken from Table 1 of
+    :cite:t:`Li:2014bc`
+
+    .. csv-table::
+        :header: "Attribute", "Default", "Unit"
+        :widths: 15, 10, 30
+
+        a_hd, 116.0, -
+        ca_ratio, 390.43, -
+        h_max,  25.33, m
+        rho_s, 200.0, kg Cm-3
+        lai,  1.8, -
+        sla,  14.0, m2 kg-1 C
+        tau_f, 4.0, years
+        tau_r,  1.04, years
+        par_ext, 0.5, -
+        yld, 0.17, -
+        zeta, 0.17, kg C m-2
+        resp_r,  0.913, year-1
+        resp_s,  0.044, year-1
+        resp_f,  0.1, -
+        m, 2, -
+        n, 5, -
+    """
+
+    a_hd: float = 116.0
+    ca_ratio: float = 390.43
+    h_max: float = 25.33
+    rho_s: float = 200.0
+    lai: float = 1.8
+    sla: float = 14.0
+    tau_f: float = 4.0
+    tau_r: float = 1.04
+    par_ext: float = 0.5
+    yld: float = 0.17
+    zeta: float = 0.17
+    resp_r: float = 0.913
+    resp_s: float = 0.044
+    resp_f: float = 0.1
+    m: float = 2
+    n: float = 5
+
+
+PlantFunctionalTypeSchema = marshmallow_dataclass.class_schema(
+    PlantFunctionalTypeStrict
+)
+"""Marshmallow validation schema class for validating PlantFunctionalType data.
+
+This schema explicitly uses the strict version of the dataclass, which enforces complete
+descriptions of plant functional type data rather than allowing partial data and filling
+in gaps from the default values.
+"""
+
+
+def calculate_q_m(m: float, n: float) -> float:
+    """Calculate a q_m value.
+
+    The value of q_m is a constant canopy scaling parameter derived from the ``m`` and
+    ``n`` attributes defined for a plant functional type.
+
+    Args:
+        m: Canopy shape parameter
+        n: Canopy shape parameter
+    """
+    return (
+        m
+        * n
+        * ((n - 1) / (m * n - 1)) ** (1 - 1 / n)
+        * (((m - 1) * n) / (m * n - 1)) ** (m - 1)
+    )
+
+
+def calculate_z_max_proportion(m: float, n: float) -> float:
+    """Calculate the z_m proportion.
+
+    The z_m proportion is the constant proportion of stem height at which the maximum
+    crown radius is found for a given plant functional type.
+
+    Args:
+        m: Canopy shape parameter
+        n: Canopy shape parameter
+    """
+
+    return ((n - 1) / (m * n - 1)) ** (1 / n)
+
+
+class Flora(dict[str, type[PlantFunctionalTypeStrict]]):
+    """Defines the flora used in a ``virtual_ecosystem`` model.
+
+    The flora is the set of plant functional types used within a particular simulation
+    and this class provides dictionary-like access to a defined set of
+    :class:`~pyrealm.demography.flora.PlantFunctionalType` or
+    :class:`~pyrealm.demography.flora.PlantFunctionalTypeStrict` instances.
+
+    Instances of this class should not be altered during model fitting, at least until
+    the point where plant evolution is included in the modelling process.
+
+    Args:
+        pfts: A sequence of ``PlantFunctionalType`` or ``PlantFunctionalTypeStrict``
+            instances, which must not have duplicated
+            :attr:`~pyrealm.demography.flora.PlantFunctionalTypeStrict.name` attributes.
+    """
+
+    def __init__(self, pfts: Sequence[type[PlantFunctionalTypeStrict]]) -> None:
+        # Initialise the dict superclass to implement dict like behaviour
+        super().__init__()
+
+        # Check the PFT data
+        if (not isinstance(pfts, Sequence)) or (
+            not all([isinstance(v, PlantFunctionalTypeStrict) for v in pfts])
+        ):
+            raise ValueError(
+                "The pfts argument must be a sequence of PlantFunctionalType instances"
+            )
+
+        # Validate the PFT instances - check there are no duplicate PFT names.
+        pft_names = Counter([p.name for p in pfts])
+        duplicates = [k for k, v in pft_names.items() if v > 1]
+
+        if duplicates:
+            raise ValueError(
+                f"Duplicated plant functional type names: {','.join(duplicates)}"
+            )
+
+        # Populate the dictionary using the PFT name as key
+        for name, pft in zip(pft_names, pfts):
+            self[name] = pft
+
+    @classmethod
+    def _from_file_data(cls, file_data: dict) -> Flora:
+        """Create a Flora object from a JSON string.
+
+        Args:
+            file_data: The payload from a data file defining plant functional types.
+        """
+        try:
+            pfts = PlantFunctionalTypeSchema().load(file_data["pft"], many=True)  # type: ignore[attr-defined]
+        except ValidationError as excep:
+            raise excep
+
+        return cls(pfts=pfts)
+
+    @classmethod
+    def from_json(cls, path: Path) -> Flora:
+        """Create a Flora object from a JSON file.
+
+        Args:
+            path: A path to a JSON file of plant functional type definitions.
+        """
+
+        try:
+            file_data = json.load(open(path))
+        except (FileNotFoundError, json.JSONDecodeError) as excep:
+            raise excep
+
+        return cls._from_file_data(file_data=file_data)
+
+    @classmethod
+    def from_toml(cls, path: Path) -> Flora:
+        """Create a Flora object from a TOML file.
+
+        Args:
+            path: A path to a TOML file of plant functional type definitions.
+        """
+
+        try:
+            file_data = tomllib.load(open(path, "rb"))
+        except (FileNotFoundError, TOMLDecodeError) as excep:
+            raise excep
+
+        return cls._from_file_data(file_data)
+
+    @classmethod
+    def from_csv(cls, path: Path) -> Flora:
+        """Create a Flora object from a CSV file.
+
+        Args:
+            path: A path to a CSV file of plant functional type definitions.
+        """
+
+        try:
+            data = pd.read_csv(path)
+        except (FileNotFoundError, pd.errors.ParserError) as excep:
+            raise excep
+
+        return cls._from_file_data({"pft": data.to_dict(orient="records")})
diff --git a/pyrealm_build_data/community/pfts.csv b/pyrealm_build_data/community/pfts.csv
new file mode 100644
index 00000000..1f0c0aab
--- /dev/null
+++ b/pyrealm_build_data/community/pfts.csv
@@ -0,0 +1,3 @@
+name,a_hd,ca_ratio,h_max,lai,par_ext,resp_f,resp_r,resp_s,rho_s,sla,tau_f,tau_r,yld,zeta,m,n
+test1,116.0,390.43,25.33,1.8,0.5,0.1,0.913,0.044,200.0,14.0,4.0,1.04,0.17,0.17,2,5
+test2,116.0,390.43,15.33,1.8,0.5,0.1,0.913,0.044,200.0,14.0,4.0,1.04,0.17,0.17,2,5
\ No newline at end of file
diff --git a/pyrealm_build_data/community/pfts.json b/pyrealm_build_data/community/pfts.json
new file mode 100644
index 00000000..5cafa9d2
--- /dev/null
+++ b/pyrealm_build_data/community/pfts.json
@@ -0,0 +1,41 @@
+{"pft": [
+  {
+    "name": "test1",
+    "a_hd": 116.0,
+    "ca_ratio": 390.43,
+    "h_max": 25.33,
+    "lai": 1.8,
+    "par_ext": 0.5,
+    "resp_f": 0.1,
+    "resp_r": 0.913,
+    "resp_s": 0.044,
+    "rho_s": 200.0,
+    "sla": 14.0,
+    "tau_f": 4.0,
+    "tau_r": 1.04,
+    "yld": 0.17,
+    "zeta": 0.17,
+    "m": 2,
+    "n": 5
+  },
+  {
+    "name": "test2",
+    "a_hd": 116.0,
+    "ca_ratio": 390.43,
+    "h_max": 15.33,
+    "lai": 1.8,
+    "par_ext": 0.5,
+    "resp_f": 0.1,
+    "resp_r": 0.913,
+    "resp_s": 0.044,
+    "rho_s": 200.0,
+    "sla": 14.0,
+    "tau_f": 4.0,
+    "tau_r": 1.04,
+    "yld": 0.17,
+    "zeta": 0.17,
+    "m": 2,
+    "n": 5
+  }
+]
+}
\ No newline at end of file
diff --git a/pyrealm_build_data/community/pfts.toml b/pyrealm_build_data/community/pfts.toml
new file mode 100644
index 00000000..4e31c832
--- /dev/null
+++ b/pyrealm_build_data/community/pfts.toml
@@ -0,0 +1,37 @@
+[[pft]]
+a_hd = 116.0
+ca_ratio = 390.43
+h_max = 25.33
+lai = 1.8
+name = 'test1'
+par_ext = 0.5
+resp_f = 0.1
+resp_r = 0.913
+resp_s = 0.044
+rho_s = 200.0
+sla = 14.0
+tau_f = 4.0
+tau_r = 1.04
+yld = 0.17
+zeta = 0.17
+m = 2
+n = 5
+
+[[pft]]
+a_hd = 116.0
+ca_ratio = 390.43
+h_max = 15.33
+lai = 1.8
+name = 'test2'
+par_ext = 0.5
+resp_f = 0.1
+resp_r = 0.913
+resp_s = 0.044
+rho_s = 200.0
+sla = 14.0
+tau_f = 4.0
+tau_r = 1.04
+yld = 0.17
+zeta = 0.17
+m = 2
+n = 5
diff --git a/pyrealm_build_data/community/pfts_invalid.csv b/pyrealm_build_data/community/pfts_invalid.csv
new file mode 100644
index 00000000..e1e94908
--- /dev/null
+++ b/pyrealm_build_data/community/pfts_invalid.csv
@@ -0,0 +1,3 @@
+name,a_hd,ca_ratio,h_max,laindex,par_ext,resp_f,resp_r,resp_s,rho_s,sla,tau_f,tau_r,yld,zeta,m,n
+test1,116.0,390.43,25.33,1.8,0.5,0.1,0.913,0.044,200.0,14.0,4.0,1.04,0.17,0.17,a,5
+test2,116.0,390.43,15.33,1.8,0.5,0.1,0.913,0.044,200.0,14.0,4.0,1.04,0.17,0.17,a,5
\ No newline at end of file
diff --git a/pyrealm_build_data/community/pfts_invalid.json b/pyrealm_build_data/community/pfts_invalid.json
new file mode 100644
index 00000000..9bececdc
--- /dev/null
+++ b/pyrealm_build_data/community/pfts_invalid.json
@@ -0,0 +1,41 @@
+{"pft": [
+  {
+    "name": "test1",
+    "a_hd": 116.0,
+    "ca_ratio": 390.43,
+    "h_max": 25.33,
+    "laindex": 1.8,
+    "par_ext": 0.5,
+    "resp_f": 0.1,
+    "resp_r": 0.913,
+    "resp_s": 0.044,
+    "rho_s": 200.0,
+    "sla": 14.0,
+    "tau_f": 4.0,
+    "tau_r": 1.04,
+    "yld": 0.17,
+    "zeta": 0.17,
+    "m": "a",
+    "n": 5
+  },
+  {
+    "name": "test2",
+    "a_hd": 116.0,
+    "ca_ratio": 390.43,
+    "h_max": 15.33,
+    "laindex": 1.8,
+    "par_ext": 0.5,
+    "resp_f": 0.1,
+    "resp_r": 0.913,
+    "resp_s": 0.044,
+    "rho_s": 200.0,
+    "sla": 14.0,
+    "tau_f": 4.0,
+    "tau_r": 1.04,
+    "yld": 0.17,
+    "zeta": 0.17,
+    "m": "a",
+    "n": 5
+  }
+]
+}
\ No newline at end of file
diff --git a/pyrealm_build_data/community/pfts_invalid.toml b/pyrealm_build_data/community/pfts_invalid.toml
new file mode 100644
index 00000000..1993518f
--- /dev/null
+++ b/pyrealm_build_data/community/pfts_invalid.toml
@@ -0,0 +1,37 @@
+[[pft]]
+a_hd = 116.0
+ca_ratio = 390.43
+h_max = 25.33
+laindex = 1.8
+name = 'test1'
+par_ext = 0.5
+resp_f = 0.1
+resp_r = 0.913
+resp_s = 0.044
+rho_s = 200.0
+sla = 14.0
+tau_f = 4.0
+tau_r = 1.04
+yld = 0.17
+zeta = 0.17
+m = 'a'
+n = 5
+
+[[pft]]
+a_hd = 116.0
+ca_ratio = 390.43
+h_max = 15.33
+laindex = 1.8
+name = 'test2'
+par_ext = 0.5
+resp_f = 0.1
+resp_r = 0.913
+resp_s = 0.044
+rho_s = 200.0
+sla = 14.0
+tau_f = 4.0
+tau_r = 1.04
+yld = 0.17
+zeta = 0.17
+m = 'a'
+n = 5
diff --git a/pyrealm_build_data/community/pfts_partial.csv b/pyrealm_build_data/community/pfts_partial.csv
new file mode 100644
index 00000000..2b16910b
--- /dev/null
+++ b/pyrealm_build_data/community/pfts_partial.csv
@@ -0,0 +1,3 @@
+name,a_hd,ca_ratio,h_max,par_ext,resp_f,resp_r,resp_s,rho_s,sla,tau_f,tau_r,yld,zeta,m
+test1,116.0,390.43,25.33,0.5,0.1,0.913,0.044,200.0,14.0,4.0,1.04,0.17,0.17,2
+test2,116.0,390.43,15.33,0.5,0.1,0.913,0.044,200.0,14.0,4.0,1.04,0.17,0.17,2
\ No newline at end of file
diff --git a/pyrealm_build_data/community/pfts_partial.json b/pyrealm_build_data/community/pfts_partial.json
new file mode 100644
index 00000000..a51371f6
--- /dev/null
+++ b/pyrealm_build_data/community/pfts_partial.json
@@ -0,0 +1,40 @@
+{
+  "pft": [
+    {
+      "a_hd": 116.0,
+      "ca_ratio": 390.43,
+      "h_max": 25.33,
+      "name": "test1",
+      "par_ext": 0.5,
+      "resp_f": 0.1,
+      "resp_r": 0.913,
+      "resp_s": 0.044,
+      "rho_s": 200.0,
+      "sla": 14.0,
+      "tau_f": 4.0,
+      "tau_r": 1.04,
+      "yld": 0.17,
+      "zeta": 0.17,
+      "m": 2,
+      "n": 5
+    },
+    {
+      "a_hd": 116.0,
+      "ca_ratio": 390.43,
+      "h_max": 15.33,
+      "name": "test2",
+      "par_ext": 0.5,
+      "resp_f": 0.1,
+      "resp_r": 0.913,
+      "resp_s": 0.044,
+      "rho_s": 200.0,
+      "sla": 14.0,
+      "tau_f": 4.0,
+      "tau_r": 1.04,
+      "yld": 0.17,
+      "zeta": 0.17,
+      "m": 2,
+      "n": 5
+    }
+  ]
+}
\ No newline at end of file
diff --git a/pyrealm_build_data/community/pfts_partial.toml b/pyrealm_build_data/community/pfts_partial.toml
new file mode 100644
index 00000000..9d64a436
--- /dev/null
+++ b/pyrealm_build_data/community/pfts_partial.toml
@@ -0,0 +1,35 @@
+[[pft]]
+a_hd = 116.0
+ca_ratio = 390.43
+h_max = 25.33
+name = 'test1'
+par_ext = 0.5
+resp_f = 0.1
+resp_r = 0.913
+resp_s = 0.044
+rho_s = 200.0
+sla = 14.0
+tau_f = 4.0
+tau_r = 1.04
+yld = 0.17
+zeta = 0.17
+m = 2
+n = 5
+
+[[pft]]
+a_hd = 116.0
+ca_ratio = 390.43
+h_max = 15.33
+name = 'test2'
+par_ext = 0.5
+resp_f = 0.1
+resp_r = 0.913
+resp_s = 0.044
+rho_s = 200.0
+sla = 14.0
+tau_f = 4.0
+tau_r = 1.04
+yld = 0.17
+zeta = 0.17
+m = 2
+n = 5
diff --git a/tests/unit/demography/test_flora.py b/tests/unit/demography/test_flora.py
new file mode 100644
index 00000000..72164314
--- /dev/null
+++ b/tests/unit/demography/test_flora.py
@@ -0,0 +1,299 @@
+"""Test flora methods."""
+
+import sys
+from contextlib import nullcontext as does_not_raise
+from importlib import resources
+from json import JSONDecodeError
+
+import pytest
+from marshmallow.exceptions import ValidationError
+from pandas.errors import ParserError
+
+if sys.version_info[:2] >= (3, 11):
+    from tomllib import TOMLDecodeError
+else:
+    from tomli import TOMLDecodeError
+
+
+STRICT_PFT_ARGS = dict(
+    a_hd=116.0,
+    ca_ratio=390.43,
+    h_max=15.33,
+    lai=1.8,
+    name="broadleaf",
+    par_ext=0.5,
+    resp_f=0.1,
+    resp_r=0.913,
+    resp_s=0.044,
+    rho_s=200.0,
+    sla=14.0,
+    tau_f=4.0,
+    tau_r=1.04,
+    yld=0.17,
+    zeta=0.17,
+    m=2,
+    n=5,
+)
+"""A dictionary of the full set of arguments needed for PlantFunctionalTypeStrict."""
+
+
+#
+# Test PlantFunctionalTypeStrict dataclass
+#
+
+
+@pytest.mark.parametrize(
+    argnames="args,outcome",
+    argvalues=[
+        pytest.param(STRICT_PFT_ARGS, does_not_raise(), id="full"),
+        pytest.param({"name": "broadleaf"}, pytest.raises(TypeError), id="partial"),
+        pytest.param({}, pytest.raises(TypeError), id="empty"),
+    ],
+)
+def test_PlantFunctionalTypeStrict__init__(args, outcome):
+    """Test the plant functional type initialisation."""
+
+    from pyrealm.demography.flora import (
+        PlantFunctionalTypeStrict,
+        calculate_q_m,
+        calculate_z_max_proportion,
+    )
+
+    with outcome:
+        pft = PlantFunctionalTypeStrict(**args)
+
+        # Check name attribute and post_init attributes if instantiation succeeds.
+        if isinstance(outcome, does_not_raise):
+            assert pft.name == "broadleaf"
+            # Expected values from defaults
+            assert pft.q_m == calculate_q_m(m=2, n=5)
+            assert pft.z_max_prop == calculate_z_max_proportion(m=2, n=5)
+
+
+#
+# Test PlantFunctionalType dataclass
+#
+
+
+@pytest.mark.parametrize(
+    argnames="args,outcome",
+    argvalues=[
+        pytest.param({"name": "broadleaf"}, does_not_raise(), id="correct"),
+        pytest.param({}, pytest.raises(TypeError), id="no_name"),
+    ],
+)
+def test_PlantFunctionalType__init__(args, outcome):
+    """Test the plant functional type initialisation."""
+
+    from pyrealm.demography.flora import (
+        PlantFunctionalType,
+        calculate_q_m,
+        calculate_z_max_proportion,
+    )
+
+    with outcome:
+        pft = PlantFunctionalType(**args)
+
+        # Check name attribute and post_init attributes if instantiation succeeds.
+        if isinstance(outcome, does_not_raise):
+            assert pft.name == "broadleaf"
+            # Expected values from defaults
+            assert pft.q_m == calculate_q_m(m=2, n=5)
+            assert pft.z_max_prop == calculate_z_max_proportion(m=2, n=5)
+
+
+#
+# Test Flora initialisation
+#
+
+
+@pytest.fixture()
+def flora_inputs(request):
+    """Fixture providing flora inputs for testing.
+
+    This is using indirect parameterisation in the test largely to isolate the import of
+    PlantFunctionalType within a method and to allow a single parameterised test to
+    have a diverse set of inputs.
+    """
+
+    from pyrealm.demography.flora import PlantFunctionalType, PlantFunctionalTypeStrict
+
+    broadleaf = PlantFunctionalType(name="broadleaf")
+    conifer = PlantFunctionalType(name="conifer")
+    broadleaf_strict = PlantFunctionalTypeStrict(**STRICT_PFT_ARGS)
+    conifer_strict_args = STRICT_PFT_ARGS.copy()
+    conifer_strict_args["name"] = "conifer"
+    conifer_strict = PlantFunctionalType(**conifer_strict_args)
+
+    match request.param:
+        case "not_sequence":
+            return "Notasequence"
+        case "sequence_not_all_pfts":
+            return [1, 2, 3]
+        case "single_pft":
+            return [broadleaf]
+        case "single_pft_strict":
+            return [broadleaf_strict]
+        case "multiple_pfts":
+            return [broadleaf, conifer]
+        case "multiple_pfts_strict":
+            return [broadleaf_strict, conifer_strict]
+        case "multiple_pfts_mixed":
+            return [broadleaf_strict, conifer]
+        case "duplicated_names":
+            return [broadleaf, broadleaf]
+        case "duplicated_names_mixed":
+            return [broadleaf_strict, broadleaf]
+
+
+@pytest.mark.parametrize(
+    argnames="flora_inputs,outcome",
+    argvalues=[
+        pytest.param("not_sequence", pytest.raises(ValueError)),
+        pytest.param("sequence_not_all_pfts", pytest.raises(ValueError)),
+        pytest.param("single_pft", does_not_raise()),
+        pytest.param("single_pft_strict", does_not_raise()),
+        pytest.param("multiple_pfts", does_not_raise()),
+        pytest.param("multiple_pfts_strict", does_not_raise()),
+        pytest.param("multiple_pfts_mixed", does_not_raise()),
+        pytest.param("duplicated_names", pytest.raises(ValueError)),
+        pytest.param("duplicated_names_mixed", pytest.raises(ValueError)),
+    ],
+    indirect=["flora_inputs"],
+)
+def test_Flora__init__(flora_inputs, outcome):
+    """Test the plant functional type initialisation."""
+
+    from pyrealm.demography.flora import Flora
+
+    with outcome:
+        flora = Flora(pfts=flora_inputs)
+
+        # Simple check that PFT instances are correctly keyed by name.
+        if isinstance(outcome, does_not_raise):
+            for k, v in flora.items():
+                assert k == v.name
+
+
+#
+# Test Flora factory methods from JSON, TOML, CSV
+#
+
+
+@pytest.mark.parametrize(
+    argnames="filename,outcome",
+    argvalues=[
+        pytest.param("pfts.json", does_not_raise(), id="correct"),
+        pytest.param("pfts_partial.json", pytest.raises(ValidationError), id="partial"),
+        pytest.param("pfts.toml", pytest.raises(JSONDecodeError), id="format_wrong"),
+        pytest.param("no.pfts", pytest.raises(FileNotFoundError), id="file_missing"),
+        pytest.param("pfts_invalid.json", pytest.raises(ValidationError), id="invalid"),
+    ],
+)
+def test_flora_from_json(filename, outcome):
+    """Test JSON loading."""
+    from pyrealm.demography.flora import Flora
+
+    datapath = resources.files("pyrealm_build_data.community") / filename
+
+    with outcome:
+        flora = Flora.from_json(datapath)
+
+        if isinstance(outcome, does_not_raise):
+            # Coarse check of what got loaded
+            assert len(flora) == 2
+            for nm in ["test1", "test2"]:
+                assert nm in flora
+
+
+@pytest.mark.parametrize(
+    argnames="filename,outcome",
+    argvalues=[
+        pytest.param("pfts.toml", does_not_raise(), id="correct"),
+        pytest.param("pfts_partial.toml", pytest.raises(ValidationError), id="partial"),
+        pytest.param("pfts.json", pytest.raises(TOMLDecodeError), id="format_wrong"),
+        pytest.param("no.pfts", pytest.raises(FileNotFoundError), id="file_missing"),
+        pytest.param("pfts_invalid.toml", pytest.raises(ValidationError), id="invalid"),
+    ],
+)
+def test_flora_from_toml(filename, outcome):
+    """Test TOML loading."""
+    from pyrealm.demography.flora import Flora
+
+    datapath = resources.files("pyrealm_build_data.community") / filename
+
+    with outcome:
+        flora = Flora.from_toml(datapath)
+
+        if isinstance(outcome, does_not_raise):
+            # Coarse check of what got loaded
+            assert len(flora) == 2
+            for nm in ["test1", "test2"]:
+                assert nm in flora
+
+
+@pytest.mark.parametrize(
+    argnames="filename,outcome",
+    argvalues=[
+        pytest.param("pfts.csv", does_not_raise(), id="correct"),
+        pytest.param("pfts.json", pytest.raises(ParserError), id="format_wrong"),
+        pytest.param("no.pfts", pytest.raises(FileNotFoundError), id="file_missing"),
+        pytest.param("pfts_partial.csv", pytest.raises(ValidationError), id="partial"),
+        pytest.param("pfts_invalid.csv", pytest.raises(ValidationError), id="invalid"),
+    ],
+)
+def test_flora_from_csv(filename, outcome):
+    """Test CSV loading."""
+    from pyrealm.demography.flora import Flora
+
+    datapath = resources.files("pyrealm_build_data.community") / filename
+
+    with outcome:
+        flora = Flora.from_csv(datapath)
+
+        if isinstance(outcome, does_not_raise):
+            # Coarse check of what got loaded
+            assert len(flora) == 2
+            for nm in ["test1", "test2"]:
+                assert nm in flora
+
+
+#
+# Test PlantFunctionalType __post_init__ functions
+#
+
+
+@pytest.mark.parametrize(
+    argnames="m,n,q_m",
+    argvalues=[(2, 5, 2.9038988210485766), (3, 4, 2.3953681843215673)],
+)
+def test_calculate_q_m(m, n, q_m):
+    """Test calculation of q_m."""
+
+    from pyrealm.demography.flora import calculate_q_m
+
+    calculated_q_m = calculate_q_m(m, n)
+    assert calculated_q_m == pytest.approx(q_m)
+
+
+def test_calculate_q_m_values_raises_exception_for_invalid_input():
+    """Test unhappy path for calculating q_m.
+
+    Test that an exception is raised when invalid arguments are provided to the
+    function.
+    """
+
+    pass
+
+
+@pytest.mark.parametrize(
+    argnames="m,n,z_max_ratio",
+    argvalues=[(2, 5, 0.8502830004171938), (3, 4, 0.7226568811456053)],
+)
+def test_calculate_z_max_ratio(m, n, z_max_ratio):
+    """Test calculation of z_max proportion."""
+
+    from pyrealm.demography.flora import calculate_z_max_proportion
+
+    calculated_zmr = calculate_z_max_proportion(m, n)
+    assert calculated_zmr == pytest.approx(z_max_ratio)