diff --git a/.github/PULL_REQUEST_TEMPLATE.md b/.github/PULL_REQUEST_TEMPLATE.md new file mode 100644 index 0000000000..3b35cc8d09 --- /dev/null +++ b/.github/PULL_REQUEST_TEMPLATE.md @@ -0,0 +1,44 @@ + + +### Description: + + + +### Collaborators: + + + + +### Expectation of Answer Changes: + + + + + +### Checklist: + + +- [ ] My change requires a change to the documentation. +- [ ] I have updated the in-code documentation AND wiki accordingly. +- [ ] I have read the [**CONTRIBUTING**](https://github.com/rgknox/fates/blob/rgknox-new-PR-template/CONTRIBUTING.md) document. +- [ ] FATES PASS/FAIL regression tests were run +- [ ] If answers were expected to change, evaluation was performed and provided + + +### Test Results: + + + + +FATES-CLM (or) HLM test hash-tag: + +FATES-CLM (or) HLM baseline hash-tag: + +FATES baseline hash-tag: + +Test Output: + + + + + \ No newline at end of file diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md new file mode 100644 index 0000000000..95028b56eb --- /dev/null +++ b/CONTRIBUTING.md @@ -0,0 +1,52 @@ +# How to contribute + +Thank you for considering contributing to the development of FATES. There are a few guidelines that we need contributions to follow. + +## Read and Understand the FATES license: + +https://github.com/NGEET/fates/blob/master/LICENSE.txt + +## If working with unreleased code, please read the use policy that governs the developer repository: + +https://github.com/NGEET/fates/wiki/Use-and-distribution-policy-for-fates-developer-repository + +## Getting Started + +Those who wish to contribute code to FATES must have those changes integrated through the developer repository NGEET/fates. Changes that make it to public releases must go through this repository first, as well. Here are some basic first steps. + +* All developers should create a fork of the NGEET/fates repository into their personal space on github +* Follow the developer work-flow described here: https://github.com/NGEET/fates/wiki/FATES-Development-Workflow +* Each set of changes should have its own feature branch that encapsulates your desired changes, following the conventions outlined here: https://github.com/NGEET/fates/wiki/Feature-Branch-Naming-Convention +* The work-flow will lead you eventually to submit a Pull-Request to NGEET/fates:master, please follow the template in the Pull Request and communicate as best you can if you are unsure how to fill out the text +* It is best to create an issue to describe the work you are undertaking prior to starting. This helps the community sync with your efforts, prevents duplication of efforts, and science is not done in a vaccuum! +* Expect peers to interact, help, discuss and eventually approve your submission (pull-request) + + +## Things to Remember + +* Make commits in logical units (i.e. group changes) +* Changes that are submitted should be limited to 1 single feature (i.e. don't submit changes to the radiation code and the nutrient cycle simultaneous, pick one thing) +* Check for unnecessary whitespace with `git diff --check` before committing +* We have no standard protocol for commit messages, but try to make them meaningful, concise and succinct. +* You will most likely have to test (see workflow above), see: https://github.com/NGEET/fates/wiki/Testing-Protocols + + +## Coding Practices and Style + +Please refer to the FATES style guide: https://github.com/NGEET/fates/wiki/Coding-Practices-and-Style-Guide + + +## Trivial Changes + +If changes are trivial, it's possible testing will not be required. Conversations via the Pull Request will address if tests are not needed + +## Documentation + +Yes please! If you are creating new code, fixing existing code, anything. Please add comments in the code itself. Please also follow the style guide for comments. Also, please create and/or modify existing wiki documentation. You may be asked to add documtation prior to having a pull-request approved. + + +## Additional Resources + +* [General GitHub documentation](https://help.github.com/) +* [GitHub pull request documentation](https://help.github.com/articles/creating-a-pull-request/) +* [FATES Wiki](https://github.com/NGEET/fates/wiki) diff --git a/LICENSE.txt b/LICENSE.txt index 12f41fb609..be5cae20af 100644 --- a/LICENSE.txt +++ b/LICENSE.txt @@ -1,6 +1,6 @@ Functionally Assembled Terrestrial Ecosystem Simulator (“FATES”) -Copyright (c) 2016-2017, The Regents of the University of California, through Lawrence +Copyright (c) 2016-2018, The Regents of the University of California, through Lawrence Berkeley National Laboratory, University Corporation for Atmospheric Research, Los Alamos National Security, LLC (LANS), as operator of Los Alamos National Laboratory (LANL), and President and Fellows of Harvard College. All rights reserved. diff --git a/README.md b/README.md index f95e973e84..8bc4f8b1a2 100644 --- a/README.md +++ b/README.md @@ -1,7 +1,7 @@ # NGEE-T fates repository ------------------------------ -This is the developer repository of the Next Generation Ecosystem Experiment Tropics’ (NGEE-T) model: the Functionally Assembled Terrestrial Ecosystem Simulator (FATES). +This is the developer repository of the Next Generation Ecosystem Experiment Tropics’ (NGEE-T) model: the Functionally Assembled Terrestrial Ecosystem Simulator (FATES). (https://github.com/NGEET/fates) For more information on the FATES model, see our wiki: https://github.com/NGEET/fates/wiki @@ -20,3 +20,5 @@ http://www.cesm.ucar.edu/ The NGEE-T project maintains a mirror of CLM. That software system will automatically pull in the FATES software, and is where most users should go to clone the code: https://github.com/NGEET/fates-clm + + diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 old mode 100755 new mode 100644 index 9c769255a4..1698c5369b --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -8,8 +8,10 @@ module EDCanopyStructureMod use FatesConstantsMod , only : r8 => fates_r8 use FatesGlobals , only : fates_log use EDPftvarcon , only : EDPftvarcon_inst - use EDGrowthFunctionsMod , only : c_area + use FatesAllometryMod , only : carea_allom use EDCohortDynamicsMod , only : copy_cohort, terminate_cohorts, fuse_cohorts + use EDCohortDynamicsMod , only : tree_lai + use EDCohortDynamicsMod , only : tree_sai use EDtypesMod , only : ed_site_type, ed_patch_type, ed_cohort_type, ncwd use EDTypesMod , only : nclmax use EDTypesMod , only : nlevleaf @@ -317,7 +319,9 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) rankordered_area_sofar = 0.0_r8 currentCohort => currentPatch%tallest do while (associated(currentCohort)) - currentCohort%c_area = c_area(currentCohort) + + call carea_allom(currentCohort%dbh,currentCohort%n,currentSite%spread,currentCohort%pft,currentCohort%c_area) + if(arealayer > currentPatch%area.and.currentCohort%canopy_layer == i_lyr)then if (ED_val_comp_excln .ge. 0.0_r8 ) then ! normal (stochastic) case. weight cohort demotion by inverse size to a constant power @@ -443,9 +447,11 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) currentCohort%n = 0.0_r8 currentCohort%c_area = 0._r8 else - currentCohort%c_area = c_area(currentCohort) + call carea_allom(currentCohort%dbh,currentCohort%n,currentSite%spread,currentCohort%pft, & + currentCohort%c_area) endif - copyc%c_area = c_area(copyc) + + call carea_allom(copyc%dbh,copyc%n,currentSite%spread,copyc%pft,copyc%c_area) !----------- Insert copy into linked list ------------------------! @@ -517,7 +523,8 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) currentCohort%c_area = 0._r8 else - currentCohort%c_area = c_area(currentCohort) + call carea_allom(currentCohort%dbh,currentCohort%n,currentSite%spread, & + currentCohort%pft,currentCohort%c_area) endif endif ! matches: if (cc_loss < currentCohort%c_area)then @@ -619,8 +626,7 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) do while (associated(currentCohort)) if(currentCohort%canopy_layer == i_lyr+1)then !look at the cohorts in the canopy layer below... currentCohort%canopy_layer = i_lyr - currentCohort%c_area = c_area(currentCohort) - + call carea_allom(currentCohort%dbh,currentCohort%n,currentSite%spread,currentCohort%pft,currentCohort%c_area) ! keep track of number and biomass of promoted cohort currentSite%promotion_rate(currentCohort%size_class) = & currentSite%promotion_rate(currentCohort%size_class) + currentCohort%n @@ -639,7 +645,7 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) ! This is the opposite of the demotion weighting... currentCohort => currentPatch%tallest do while (associated(currentCohort)) - currentCohort%c_area = c_area(currentCohort) + call carea_allom(currentCohort%dbh,currentCohort%n,currentSite%spread,currentCohort%pft,currentCohort%c_area) if(currentCohort%canopy_layer == i_lyr+1)then !look at the cohorts in the canopy layer below... if (ED_val_comp_excln .ge. 0.0_r8 ) then ! normal (stochastic) case, as above. @@ -720,8 +726,8 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) currentCohort%dbh = currentCohort%dbh - 0.000000000001_r8 copyc%dbh = copyc%dbh + 0.000000000001_r8 - currentCohort%c_area = c_area(currentCohort) - copyc%c_area = c_area(copyc) + call carea_allom(currentCohort%dbh,currentCohort%n,currentSite%spread,currentCohort%pft,currentCohort%c_area) + call carea_allom(copyc%dbh,copyc%n,currentSite%spread,copyc%pft,copyc%c_area) !----------- Insert copy into linked list ------------------------! copyc%shorter => currentCohort @@ -738,8 +744,8 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) ! update area AFTER we sum up the losses. the cohort may shrink at this point, ! if the upper canopy spread is smaller. this shold be dealt with by the 'excess area' loop. - currentCohort%c_area = c_area(currentCohort) - + call carea_allom(currentCohort%dbh,currentCohort%n,currentSite%spread,currentCohort%pft,currentCohort%c_area) + ! keep track of number and biomass of promoted cohort currentSite%promotion_rate(currentCohort%size_class) = & currentSite%promotion_rate(currentCohort%size_class) + currentCohort%n @@ -826,7 +832,7 @@ subroutine canopy_spread( currentSite ) !calculate canopy area in each patch... currentCohort => currentPatch%tallest do while (associated(currentCohort)) - currentCohort%c_area = c_area(currentCohort) + call carea_allom(currentCohort%dbh,currentCohort%n,currentSite%spread,currentCohort%pft,currentCohort%c_area) if(EDPftvarcon_inst%woody(currentCohort%pft) .eq. 1 .and. currentCohort%canopy_layer .eq. 1 ) then sitelevel_canopyarea = sitelevel_canopyarea + currentCohort%c_area endif @@ -862,7 +868,6 @@ subroutine canopy_summarization( nsites, sites, bc_in ) use EDPatchDynamicsMod , only : set_patchno use EDPatchDynamicsMod , only : set_root_fraction use FatesSizeAgeTypeIndicesMod, only : sizetype_class_index - use EDGrowthFunctionsMod , only : tree_lai, c_area use EDtypesMod , only : area use EDPftvarcon , only : EDPftvarcon_inst @@ -921,9 +926,9 @@ subroutine canopy_summarization( nsites, sites, bc_in ) currentCohort%b = currentCohort%balive+currentCohort%bdead+currentCohort%bstore + call carea_allom(currentCohort%dbh,currentCohort%n,sites(s)%spread,currentCohort%pft,currentCohort%c_area) currentCohort%treelai = tree_lai(currentCohort) - - currentCohort%c_area = c_area(currentCohort) + canopy_leaf_area = canopy_leaf_area + currentCohort%treelai *currentCohort%c_area if(currentCohort%canopy_layer==1)then @@ -976,7 +981,6 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) ! ! !USES: - use EDGrowthFunctionsMod , only : tree_lai, tree_sai, c_area use EDtypesMod , only : area, dinc_ed, hitemax, n_hite_bins ! @@ -1026,8 +1030,8 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) currentPatch%canopy_layer_lai(:) = 0._r8 NC = 0 currentCohort => currentPatch%shortest - do while(associated(currentCohort)) - currentCohort%c_area = c_area(currentCohort) + do while(associated(currentCohort)) + call carea_allom(currentCohort%dbh,currentCohort%n,currentSite%spread,currentCohort%pft,currentCohort%c_area) currentPatch%canopy_area = currentPatch%canopy_area + currentCohort%c_area NC = NC+1 currentCohort => currentCohort%taller @@ -1594,8 +1598,8 @@ subroutine CanopyLayerArea(currentPatch,layer_index,layer_area) layer_area = 0.0_r8 currentCohort => currentPatch%tallest - do while (associated(currentCohort)) - currentCohort%c_area = c_area(currentCohort) ! Reassess cohort area. + do while (associated(currentCohort)) + call carea_allom(currentCohort%dbh,currentCohort%n,currentPatch%siteptr%spread,currentCohort%pft,currentCohort%c_area) if (currentCohort%canopy_layer .eq. layer_index) then layer_area = layer_area + currentCohort%c_area end if @@ -1623,7 +1627,7 @@ function NumPotentialCanopyLayers(currentPatch,include_substory) result(z) type(ed_cohort_type),pointer :: currentCohort integer :: z - + real(r8) :: c_area real(r8) :: arealayer z = 1 @@ -1638,7 +1642,8 @@ function NumPotentialCanopyLayers(currentPatch,include_substory) result(z) currentCohort => currentPatch%tallest do while (associated(currentCohort)) if(currentCohort%canopy_layer == z) then - arealayer = arealayer + c_area(currentCohort) + call carea_allom(currentCohort%dbh,currentCohort%n,currentPatch%siteptr%spread,currentCohort%pft,c_area) + arealayer = arealayer + c_area end if currentCohort => currentCohort%shorter enddo diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 old mode 100755 new mode 100644 index 105df95f78..cc20f45f88 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -10,10 +10,9 @@ module EDCohortDynamicsMod use FatesInterfaceMod , only : bc_in_type use FatesConstantsMod , only : r8 => fates_r8 use FatesConstantsMod , only : fates_unset_int - use FatesConstantsMod , only : itrue + use FatesConstantsMod , only : itrue,ifalse use FatesInterfaceMod , only : hlm_days_per_year use EDPftvarcon , only : EDPftvarcon_inst - use EDGrowthFunctionsMod , only : c_area, tree_lai use EDTypesMod , only : ed_site_type, ed_patch_type, ed_cohort_type use EDTypesMod , only : nclmax use EDTypesMod , only : ncwd @@ -21,6 +20,8 @@ module EDCohortDynamicsMod use EDTypesMod , only : AREA use EDTypesMod , only : min_npm2, min_nppatch use EDTypesMod , only : min_n_safemath + use EDTypesMod , only : nlevleaf + use EDTypesMod , only : dinc_ed use FatesInterfaceMod , only : hlm_use_planthydro use FatesPlantHydraulicsMod, only : FuseCohortHydraulics use FatesPlantHydraulicsMod, only : CopyCohortHydraulics @@ -29,7 +30,11 @@ module EDCohortDynamicsMod use FatesPlantHydraulicsMod, only : InitHydrCohort use FatesPlantHydraulicsMod, only : DeallocateHydrCohort use FatesSizeAgeTypeIndicesMod, only : sizetype_class_index - + use FatesAllometryMod , only : bsap_allom + use FatesAllometryMod , only : bleaf + use FatesAllometryMod , only : bfineroot + use FatesAllometryMod , only : h_allom + use FatesAllometryMod , only : carea_allom ! CIME globals use shr_log_mod , only : errMsg => shr_log_errMsg @@ -47,6 +52,8 @@ module EDCohortDynamicsMod public :: copy_cohort public :: count_cohorts public :: allocate_live_biomass + public :: tree_lai + public :: tree_sai logical, parameter :: DEBUG = .false. ! local debug flag @@ -143,7 +150,8 @@ subroutine create_cohort(patchptr, pft, nn, hite, dbh, & call allocate_live_biomass(new_cohort,0) ! Assign canopy extent and depth - new_cohort%c_area = c_area(new_cohort) + call carea_allom(new_cohort%dbh,new_cohort%n,new_cohort%siteptr%spread,new_cohort%pft,new_cohort%c_area) + new_cohort%treelai = tree_lai(new_cohort) new_cohort%lai = new_cohort%treelai * new_cohort%c_area/patchptr%area new_cohort%treesai = 0.0_r8 !FIX(RF,032414) @@ -210,14 +218,28 @@ subroutine allocate_live_biomass(cc_p,mode) real(r8) :: new_bl real(r8) :: new_br real(r8) :: new_bsw + real(r8) :: tar_bl ! target leaf biomass when leaves are flushed (includes trimming) + real(r8) :: tar_br ! target fineroot biomass (includes trimming) + real(r8) :: tar_bsw ! target sapwood biomass + real(r8) :: bfr_per_leaf ! ratio of fine roots to leaf mass when plants are on allometry + real(r8) :: bsw_per_leaf ! ratio of sapwood to leaf mass when plants are on allometry + real(r8) :: temp_h + integer :: ft ! functional type integer :: leaves_off_switch !---------------------------------------------------------------------- currentCohort => cc_p ft = currentcohort%pft - leaf_frac = 1.0_r8/(1.0_r8 + EDpftvarcon_inst%allom_latosa_int(ft) * currentcohort%hite + EDPftvarcon_inst%allom_l2fr(ft)) + + call bleaf(currentcohort%dbh,currentcohort%hite,ft,currentcohort%canopy_trim,tar_bl) + call bfineroot(currentcohort%dbh,currentcohort%hite,ft,currentcohort%canopy_trim,tar_br) + call bsap_allom(currentcohort%dbh,ft,currentcohort%canopy_trim,tar_bsw) + + leaf_frac = tar_bl/(tar_bl+tar_br+tar_bsw) + bfr_per_leaf = tar_br/tar_bl + bsw_per_leaf = tar_bsw/tar_bl !currentcohort%bl = currentcohort%balive*leaf_frac !for deciduous trees, there are no leaves @@ -242,14 +264,13 @@ subroutine allocate_live_biomass(cc_p,mode) endif ! Use different proportions if the leaves are on vs off - if(leaves_off_switch==0)then + if(leaves_off_switch.eq.ifalse)then ! leaves are on new_bl = currentcohort%balive*leaf_frac - new_br = EDpftvarcon_inst%allom_l2fr(ft) * (currentcohort%balive + currentcohort%laimemory) * leaf_frac + new_br = bfr_per_leaf * (currentcohort%balive + currentcohort%laimemory) * leaf_frac - new_bsw = EDpftvarcon_inst%allom_latosa_int(ft) * currentcohort%hite *(currentcohort%balive + & - currentcohort%laimemory)*leaf_frac + new_bsw = bsw_per_leaf * (currentcohort%balive + currentcohort%laimemory) * leaf_frac !diagnose the root and stem biomass from the functional balance hypothesis. This is used when the leaves are !fully on. @@ -271,7 +292,7 @@ subroutine allocate_live_biomass(cc_p,mode) currentcohort%br = new_br currentcohort%bsw = new_bsw - else ! Leaves are off (leaves_off_switch==1) + else ! Leaves are off (leaves_off_switch==.itrue.) !the purpose of this section is to figure out the root and stem biomass when the leaves are off !at this point, we know the former leaf mass (laimemory) and the current alive mass @@ -280,14 +301,13 @@ subroutine allocate_live_biomass(cc_p,mode) !not have enough live biomass to support the hypothesized root mass !thus, we use 'ratio_balive' to adjust br and bsw. Apologies that this is so complicated! RF - ideal_balive = currentcohort%laimemory * EDPftvarcon_inst%allom_l2fr(ft) + & - currentcohort%laimemory* EDpftvarcon_inst%allom_latosa_int(ft) * currentcohort%hite + ideal_balive = currentcohort%laimemory * bfr_per_leaf + currentcohort%laimemory * bsw_per_leaf + ratio_balive = currentcohort%balive / ideal_balive - new_br = EDpftvarcon_inst%allom_l2fr(ft) * (ideal_balive + currentcohort%laimemory) * & - leaf_frac * ratio_balive - new_bsw = EDpftvarcon_inst%allom_latosa_int(ft) * currentcohort%hite * & - (ideal_balive + currentcohort%laimemory) * leaf_frac * ratio_balive + new_br = bfr_per_leaf * (ideal_balive + currentcohort%laimemory) * leaf_frac * ratio_balive + + new_bsw = bsw_per_leaf * (ideal_balive + currentcohort%laimemory) * leaf_frac * ratio_balive ! Diagnostics if(mode==1)then @@ -670,7 +690,6 @@ subroutine fuse_cohorts(patchptr, bc_in) ! Join similar cohorts to reduce total number ! ! !USES: - use EDTypesMod , only : nlevleaf use EDParamsMod , only : ED_val_cohort_fusion_tol use shr_infnan_mod, only : nan => shr_infnan_nan, assignment(=) ! @@ -1334,7 +1353,83 @@ function count_cohorts( currentPatch ) result ( backcount ) end function count_cohorts + ! ===================================================================================== + + real(r8) function tree_lai( cohort_in ) + ! ============================================================================ + ! LAI of individual trees is a function of the total leaf area and the total canopy area. + ! ============================================================================ + + type(ed_cohort_type), intent(inout) :: cohort_in + + real(r8) :: leafc_per_unitarea ! KgC of leaf per m2 area of ground. + real(r8) :: slat ! the sla of the top leaf layer. m2/kgC + + if( cohort_in%bl < 0._r8 .or. cohort_in%pft == 0 ) then + write(fates_log(),*) 'problem in treelai',cohort_in%bl,cohort_in%pft + endif + + if( cohort_in%status_coh == 2 ) then ! are the leaves on? + slat = 1000.0_r8 * EDPftvarcon_inst%slatop(cohort_in%pft) ! m2/g to m2/kg + leafc_per_unitarea = cohort_in%bl/(cohort_in%c_area/cohort_in%n) !KgC/m2 + if(leafc_per_unitarea > 0.0_r8)then + tree_lai = leafc_per_unitarea * slat !kg/m2 * m2/kg = unitless LAI + else + tree_lai = 0.0_r8 + endif + else + tree_lai = 0.0_r8 + endif !status + cohort_in%treelai = tree_lai + + ! here, if the LAI exceeeds the maximum size of the possible array, then we have no way of accomodating it + ! at the moments nlevleaf default is 40, which is very large, so exceeding this would clearly illustrate a + ! huge error + if(cohort_in%treelai > nlevleaf*dinc_ed)then + write(fates_log(),*) 'too much lai' , cohort_in%treelai , cohort_in%pft , nlevleaf * dinc_ed + endif + + return + + end function tree_lai + + ! ============================================================================ + + real(r8) function tree_sai( cohort_in ) + + ! ============================================================================ + ! SAI of individual trees is a function of the total dead biomass per unit canopy area. + ! ============================================================================ + + type(ed_cohort_type), intent(inout) :: cohort_in + + real(r8) :: bdead_per_unitarea ! KgC of leaf per m2 area of ground. + real(r8) :: sai_scaler + + sai_scaler = EDPftvarcon_inst%allom_sai_scaler(cohort_in%pft) + + if( cohort_in%bdead < 0._r8 .or. cohort_in%pft == 0 ) then + write(fates_log(),*) 'problem in treesai',cohort_in%bdead,cohort_in%pft + endif + + bdead_per_unitarea = cohort_in%bdead/(cohort_in%c_area/cohort_in%n) !KgC/m2 + tree_sai = bdead_per_unitarea * sai_scaler !kg/m2 * m2/kg = unitless LAI + + cohort_in%treesai = tree_sai + + ! here, if the LAI exceeeds the maximum size of the possible array, then we have no way of accomodating it + ! at the moments nlevleaf default is 40, which is very large, so exceeding this would clearly illustrate a + ! huge error + if(cohort_in%treesai > nlevleaf*dinc_ed)then + write(fates_log(),*) 'too much sai' , cohort_in%treesai , cohort_in%pft , nlevleaf * dinc_ed + endif + + return + + end function tree_sai + + ! ============================================================================ !-------------------------------------------------------------------------------------! diff --git a/biogeochem/EDGrowthFunctionsMod.F90 b/biogeochem/EDGrowthFunctionsMod.F90 deleted file mode 100755 index 075c216a33..0000000000 --- a/biogeochem/EDGrowthFunctionsMod.F90 +++ /dev/null @@ -1,470 +0,0 @@ -module EDGrowthFunctionsMod - - ! ============================================================================ - ! Functions that control the trajectory of plant growth. - ! Ideally these would all use parameters that are fed in from the parameter file. - ! At present, there is only a single allocation trajectory. - ! ============================================================================ - - use FatesConstantsMod, only : r8 => fates_r8 - use FatesGlobals , only : fates_log - use EDPftvarcon , only : EDPftvarcon_inst - use EDTypesMod , only : ed_cohort_type, nlevleaf, dinc_ed - use FatesConstantsMod , only : itrue,ifalse - - implicit none - private - - public :: bleaf - public :: hite - public :: ddbhdbd - public :: ddbhdbl - public :: dhdbd - public :: dbh - public :: bdead - public :: tree_lai - public :: tree_sai - public :: c_area - public :: mortality_rates - - logical :: DEBUG_growth = .false. - - ! ============================================================================ - ! 10/30/09: Created by Rosie Fisher - ! ============================================================================ - -contains - - real(r8) function Dbh( cohort_in ) - - ! ============================================================================ - ! Creates diameter in cm as a function of height in m - ! Height(m) diameter(cm) relationships. O'Brien et al - for 56 patch at BCI - ! ============================================================================ - - type(ed_cohort_type), intent(in) :: cohort_in - - !FIX(SPM,040214) - move to param file - real(r8) :: m ! parameter of allometric equation - real(r8) :: c ! parameter of allometric equation - - m = EDPftvarcon_inst%allom_d2h1(cohort_in%pft) - c = EDPftvarcon_inst%allom_d2h2(cohort_in%pft) - - dbh = (10.0_r8**((log10(cohort_in%hite) - c)/m)) - - return - - end function dbh - -! ============================================================================ - - real(r8) function Hite( cohort_in ) - - ! ============================================================================ - ! Creates height in m as a function of diameter in cm. - ! Height(m) diameter(cm) relationships. O'Brien et al - for 56 pft at BCI - ! ============================================================================ - - type(ed_cohort_type), intent(inout) :: cohort_in - - real(r8) :: m - real(r8) :: c - real(r8) :: h - - m = EDPftvarcon_inst%allom_d2h1(cohort_in%pft) - c = EDPftvarcon_inst%allom_d2h2(cohort_in%pft) - - if(cohort_in%dbh <= 0._r8)then - write(fates_log(),*) 'ED: dbh less than zero problem!' - cohort_in%dbh = 0.1_r8 - endif - - ! if the hite is larger than the maximum allowable height (set by dbhmax) then - ! set the height to the maximum value. - ! this could do with at least re-factoring and probably re-thinking. RF - if(cohort_in%dbh <= EDPftvarcon_inst%allom_dbh_maxheight(cohort_in%pft)) then - h = (10.0_r8**(log10(cohort_in%dbh) * m + c)) - else - h = (10.0_r8**(log10(EDPftvarcon_inst%allom_dbh_maxheight(cohort_in%pft))*m + c)) - endif - Hite = h - - return - - end function Hite - -! ============================================================================ - - real(r8) function Bleaf( cohort_in ) - - ! ============================================================================ - ! Creates leaf biomass (kGC) as a function of tree diameter. - ! ============================================================================ - - type(ed_cohort_type), intent(in) :: cohort_in - - real(r8) :: dbh2bl_a - real(r8) :: dbh2bl_b - real(r8) :: dbh2bl_c - - dbh2bl_a = EDPftvarcon_inst%allom_d2bl1(cohort_in%pft) - dbh2bl_b = EDPftvarcon_inst%allom_d2bl2(cohort_in%pft) - dbh2bl_c = EDPftvarcon_inst%allom_d2bl3(cohort_in%pft) - - if(cohort_in%dbh < 0._r8.or.cohort_in%pft == 0.or.cohort_in%dbh > 1000.0_r8)then - write(fates_log(),*) 'problems in bleaf',cohort_in%dbh,cohort_in%pft - endif - - if(cohort_in%dbh <= EDPftvarcon_inst%allom_dbh_maxheight(cohort_in%pft))then - bleaf = dbh2bl_a * (cohort_in%dbh**dbh2bl_b) * EDPftvarcon_inst%wood_density(cohort_in%pft)**dbh2bl_c - else - bleaf = dbh2bl_a * (EDPftvarcon_inst%allom_dbh_maxheight(cohort_in%pft)**dbh2bl_b) * & - EDPftvarcon_inst%wood_density(cohort_in%pft)**dbh2bl_c - endif - - !write(fates_log(),*) 'bleaf',bleaf, slascaler,cohort_in%pft - - !Adjust for canopies that have become so deep that their bottom layer is not producing any carbon... - !nb this will change the allometry and the effects of this remain untested. RF. April 2014 - - bleaf = bleaf * cohort_in%canopy_trim - - return - end function Bleaf - -! ============================================================================ - - real(r8) function tree_lai( cohort_in ) - - ! ============================================================================ - ! LAI of individual trees is a function of the total leaf area and the total canopy area. - ! ============================================================================ - - type(ed_cohort_type), intent(inout) :: cohort_in - - real(r8) :: leafc_per_unitarea ! KgC of leaf per m2 area of ground. - real(r8) :: slat ! the sla of the top leaf layer. m2/kgC - - if( cohort_in%bl < 0._r8 .or. cohort_in%pft == 0 ) then - write(fates_log(),*) 'problem in treelai',cohort_in%bl,cohort_in%pft - endif - - if( cohort_in%status_coh == 2 ) then ! are the leaves on? - slat = 1000.0_r8 * EDPftvarcon_inst%slatop(cohort_in%pft) ! m2/g to m2/kg - cohort_in%c_area = c_area(cohort_in) ! call the tree area - leafc_per_unitarea = cohort_in%bl/(cohort_in%c_area/cohort_in%n) !KgC/m2 - if(leafc_per_unitarea > 0.0_r8)then - tree_lai = leafc_per_unitarea * slat !kg/m2 * m2/kg = unitless LAI - else - tree_lai = 0.0_r8 - endif - else - tree_lai = 0.0_r8 - endif !status - cohort_in%treelai = tree_lai - - ! here, if the LAI exceeeds the maximum size of the possible array, then we have no way of accomodating it - ! at the moments nlevleaf default is 40, which is very large, so exceeding this would clearly illustrate a - ! huge error - if(cohort_in%treelai > nlevleaf*dinc_ed)then - write(fates_log(),*) 'too much lai' , cohort_in%treelai , cohort_in%pft , nlevleaf * dinc_ed - endif - - return - - end function tree_lai - - ! ============================================================================ - - real(r8) function tree_sai( cohort_in ) - - ! ============================================================================ - ! SAI of individual trees is a function of the total dead biomass per unit canopy area. - ! ============================================================================ - - type(ed_cohort_type), intent(inout) :: cohort_in - - real(r8) :: bdead_per_unitarea ! KgC of leaf per m2 area of ground. - real(r8) :: sai_scaler - - sai_scaler = EDPftvarcon_inst%allom_sai_scaler(cohort_in%pft) - - if( cohort_in%bdead < 0._r8 .or. cohort_in%pft == 0 ) then - write(fates_log(),*) 'problem in treesai',cohort_in%bdead,cohort_in%pft - endif - - cohort_in%c_area = c_area(cohort_in) ! call the tree area - bdead_per_unitarea = cohort_in%bdead/(cohort_in%c_area/cohort_in%n) !KgC/m2 - tree_sai = bdead_per_unitarea * sai_scaler !kg/m2 * m2/kg = unitless LAI - - cohort_in%treesai = tree_sai - - ! here, if the LAI exceeeds the maximum size of the possible array, then we have no way of accomodating it - ! at the moments nlevleaf default is 40, which is very large, so exceeding this would clearly illustrate a - ! huge error - if(cohort_in%treesai > nlevleaf*dinc_ed)then - write(fates_log(),*) 'too much sai' , cohort_in%treesai , cohort_in%pft , nlevleaf * dinc_ed - endif - - return - - end function tree_sai - - -! ============================================================================ - - real(r8) function c_area( cohort_in ) - - ! ============================================================================ - ! Calculate area of ground covered by entire cohort. (m2) - ! Function of DBH (cm) canopy spread (m/cm) and number of individuals. - ! ============================================================================ - - use EDTypesMod , only : nclmax - - type(ed_cohort_type), intent(in) :: cohort_in - - real(r8) :: dbh ! Tree diameter at breat height. cm. - real(r8) :: crown_area_to_dbh_exponent - real(r8) :: spreadterm - - ! default is to use the same exponent as the dbh to bleaf exponent so that per-plant canopy depth remains invariant during growth, - ! but allowed to vary via the allom_blca_expnt_diff term (which has default value of zero) - crown_area_to_dbh_exponent = EDPftvarcon_inst%allom_d2bl2(cohort_in%pft) + & - EDPftvarcon_inst%allom_blca_expnt_diff(cohort_in%pft) - - if (DEBUG_growth) then - write(fates_log(),*) 'z_area 1',cohort_in%dbh,cohort_in%pft - write(fates_log(),*) 'z_area 2',EDPftvarcon_inst%allom_dbh_maxheight - write(fates_log(),*) 'z_area 3',EDPftvarcon_inst%woody - write(fates_log(),*) 'z_area 4',cohort_in%n - write(fates_log(),*) 'z_area 5',cohort_in%siteptr%spread - write(fates_log(),*) 'z_area 6',cohort_in%canopy_layer - end if - - dbh = min(cohort_in%dbh,EDPftvarcon_inst%allom_dbh_maxheight(cohort_in%pft)) - - ! ---------------------------------------------------------------------------------- - ! The function c_area is called during the process of canopy position demotion - ! and promotion. As such, some cohorts are temporarily elevated to canopy positions - ! that are outside the number of alloted canopy spaces. Ie, a two story canopy - ! may have a third-story plant, if only for a moment. However, these plants - ! still need to generate a crown area to complete the promotion, demotion process. - ! So we allow layer index exceedence here and force it down to max. - ! (rgk/cdk 05/2017) - ! ---------------------------------------------------------------------------------- - - ! apply site-level spread elasticity to the cohort crown allometry term - spreadterm = cohort_in%siteptr%spread * EDPftvarcon_inst%allom_d2ca_coefficient_max(cohort_in%pft) + & - (1._r8 - cohort_in%siteptr%spread) * EDPftvarcon_inst%allom_d2ca_coefficient_min(cohort_in%pft) - ! - c_area = cohort_in%n * spreadterm * dbh ** crown_area_to_dbh_exponent - - end function c_area - -! ============================================================================ - - real(r8) function Bdead( cohort_in ) - - ! ============================================================================ - ! Calculate stem biomass from height(m) dbh(cm) and wood density(g/cm3) - ! default params using allometry of J.G. Saldarriaga et al 1988 - Rio Negro - ! Journal of Ecology vol 76 p938-958 - ! - ! NOTE (RGK 07-2017) Various other biomass allometries calculate above ground - ! biomass, and it appear Saldariagga may be an outlier that calculates total - ! biomass (these parameters will have to be a placeholder for both) - ! - ! ============================================================================ - - type(ed_cohort_type), intent(in) :: cohort_in - - real(r8) :: dbh2bd_a - real(r8) :: dbh2bd_b - real(r8) :: dbh2bd_c - real(r8) :: dbh2bd_d - - dbh2bd_a = EDPftvarcon_inst%allom_agb1(cohort_in%pft) - dbh2bd_b = EDPftvarcon_inst%allom_agb2(cohort_in%pft) - dbh2bd_c = EDPftvarcon_inst%allom_agb3(cohort_in%pft) - dbh2bd_d = EDPftvarcon_inst%allom_agb4(cohort_in%pft) - - bdead = dbh2bd_a*(cohort_in%hite**dbh2bd_b)*(cohort_in%dbh**dbh2bd_c)* & - (EDPftvarcon_inst%wood_density(cohort_in%pft)** dbh2bd_d) - - end function Bdead - -! ============================================================================ - - real(r8) function dHdBd( cohort_in ) - - ! ============================================================================ - ! convert changes in structural biomass to changes in height - ! consistent with Bstem and h-dbh allometries - ! ============================================================================ - - type(ed_cohort_type), intent(in) :: cohort_in - - real(r8) :: dbddh ! rate of change of dead biomass (KgC) per unit change of height (m) - real(r8) :: dbh2bd_a - real(r8) :: dbh2bd_b - real(r8) :: dbh2bd_c - real(r8) :: dbh2bd_d - - dbh2bd_a = EDPftvarcon_inst%allom_agb1(cohort_in%pft) - dbh2bd_b = EDPftvarcon_inst%allom_agb2(cohort_in%pft) - dbh2bd_c = EDPftvarcon_inst%allom_agb3(cohort_in%pft) - dbh2bd_d = EDPftvarcon_inst%allom_agb4(cohort_in%pft) - - dbddh = dbh2bd_a*dbh2bd_b*(cohort_in%hite**(dbh2bd_b-1.0_r8))*(cohort_in%dbh**dbh2bd_c)* & - (EDPftvarcon_inst%wood_density(cohort_in%pft)**dbh2bd_d) - dHdBd = 1.0_r8/dbddh !m/KgC - - return - - end function dHdBd - -! ============================================================================ - real(r8) function dDbhdBd( cohort_in ) - - ! ============================================================================ - ! convert changes in structural biomass to changes in diameter - ! consistent with Bstem and h-dbh allometries - ! ============================================================================ - - type(ed_cohort_type), intent(in) :: cohort_in - - real(r8) :: dBD_dDBH !Rate of change of dead biomass (KgC) with change in DBH (cm) - real(r8) :: dH_dDBH !Rate of change of height (m) with change in DBH (cm) - real(r8) :: m - real(r8) :: c - real(r8) :: h - real(r8) :: dbh2bd_a - real(r8) :: dbh2bd_b - real(r8) :: dbh2bd_c - real(r8) :: dbh2bd_d - - m = EDPftvarcon_inst%allom_d2h1(cohort_in%pft) - c = EDPftvarcon_inst%allom_d2h2(cohort_in%pft) - - dbh2bd_a = EDPftvarcon_inst%allom_agb1(cohort_in%pft) - dbh2bd_b = EDPftvarcon_inst%allom_agb2(cohort_in%pft) - dbh2bd_c = EDPftvarcon_inst%allom_agb3(cohort_in%pft) - dbh2bd_d = EDPftvarcon_inst%allom_agb4(cohort_in%pft) - - dBD_dDBH = dbh2bd_c*dbh2bd_a*(cohort_in%hite**dbh2bd_b)*(cohort_in%dbh**(dbh2bd_c-1.0_r8))* & - (EDPftvarcon_inst%wood_density(cohort_in%pft)**dbh2bd_d) - - if(cohort_in%dbh < EDPftvarcon_inst%allom_dbh_maxheight(cohort_in%pft))then - dH_dDBH = (10.0_r8**c)*m*(cohort_in%dbh**(m-1.0_r8)) - - dBD_dDBH = dBD_dDBH + dbh2bd_b*dbh2bd_a*(cohort_in%hite**(dbh2bd_b - 1.0_r8))* & - (cohort_in%dbh**dbh2bd_c)*(EDPftvarcon_inst%wood_density(cohort_in%pft)**dbh2bd_d)*dH_dDBH - endif - - dDbhdBd = 1.0_r8/dBD_dDBH - - return - - end function dDbhdBd - -! ============================================================================ - - real(r8) function dDbhdBl( cohort_in ) - - ! ============================================================================ - ! convert changes in leaf biomass (KgC) to changes in DBH (cm) - ! ============================================================================ - - type(ed_cohort_type), intent(in) :: cohort_in - - real(r8) :: dblddbh ! Rate of change of leaf biomass with change in DBH - real(r8) :: dbh2bl_a - real(r8) :: dbh2bl_b - real(r8) :: dbh2bl_c - - dbh2bl_a = EDPftvarcon_inst%allom_d2bl1(cohort_in%pft) - dbh2bl_b = EDPftvarcon_inst%allom_d2bl2(cohort_in%pft) - dbh2bl_c = EDPftvarcon_inst%allom_d2bl3(cohort_in%pft) - - - dblddbh = dbh2bl_b*dbh2bl_a*(cohort_in%dbh**dbh2bl_b)*(EDPftvarcon_inst%wood_density(cohort_in%pft)**dbh2bl_c) - dblddbh = dblddbh*cohort_in%canopy_trim - - if( cohort_in%dbh 0._r8 ) then - if(Bleaf(cohort_in) > 0._r8 .and. cohort_in%bstore <= Bleaf(cohort_in))then - frac = cohort_in%bstore/(Bleaf(cohort_in)) - cmort = max(0.0_r8,ED_val_stress_mort*(1.0_r8 - frac)) - else - cmort = 0.0_r8 - endif - - else - write(fates_log(),*) 'dbh problem in mortality_rates', & - cohort_in%dbh,cohort_in%pft,cohort_in%n,cohort_in%canopy_layer - endif - - !mortality_rates = bmort + hmort + cmort - - else ! i.e. hlm_use_ed_prescribed_phys is true - if ( cohort_in%canopy_layer .eq. 1) then - bmort = EDPftvarcon_inst%prescribed_mortality_canopy(cohort_in%pft) - else - bmort = EDPftvarcon_inst%prescribed_mortality_understory(cohort_in%pft) - endif - cmort = 0._r8 - hmort = 0._r8 - endif - - end subroutine mortality_rates - -! ============================================================================ - -end module EDGrowthFunctionsMod diff --git a/biogeochem/EDLoggingMortalityMod.F90 b/biogeochem/EDLoggingMortalityMod.F90 index 4f87fcaaed..d3fe5d6b1e 100644 --- a/biogeochem/EDLoggingMortalityMod.F90 +++ b/biogeochem/EDLoggingMortalityMod.F90 @@ -225,7 +225,8 @@ subroutine logging_litter_fluxes(currentSite, currentPatch, newPatch, patch_site use EDtypesMod, only : ed_site_type use EDtypesMod, only : ed_patch_type use EDtypesMod, only : ed_cohort_type - use EDGrowthFunctionsMod, only : c_area + use FatesAllometryMod , only : carea_allom + ! !ARGUMENTS: type(ed_site_type) , intent(inout), target :: currentSite @@ -451,7 +452,7 @@ subroutine logging_litter_fluxes(currentSite, currentPatch, newPatch, patch_site currentCohort => newPatch%shortest do while(associated(currentCohort)) - currentCohort%c_area = c_area(currentCohort) + call carea_allom(currentCohort%dbh,currentCohort%n,currentSite%spread,currentCohort%pft,currentCohort%c_area) currentCohort => currentCohort%taller enddo diff --git a/biogeochem/EDMortalityFunctionsMod.F90 b/biogeochem/EDMortalityFunctionsMod.F90 new file mode 100644 index 0000000000..c4e34b0c52 --- /dev/null +++ b/biogeochem/EDMortalityFunctionsMod.F90 @@ -0,0 +1,95 @@ +module EDMortalityFunctionsMod + + ! ============================================================================ + ! Functions that control mortality. + ! ============================================================================ + + use FatesConstantsMod, only : r8 => fates_r8 + use FatesGlobals , only : fates_log + use EDPftvarcon , only : EDPftvarcon_inst + use EDTypesMod , only : ed_cohort_type + use FatesConstantsMod, only : itrue,ifalse + use FatesAllometryMod, only : bleaf + use EDParamsMod , only : ED_val_stress_mort + use FatesInterfaceMod, only : hlm_use_ed_prescribed_phys + + implicit none + private + + + public :: mortality_rates + + logical :: DEBUG_growth = .false. + + ! ============================================================================ + ! 10/30/09: Created by Rosie Fisher + ! ============================================================================ + +contains + + + + subroutine mortality_rates( cohort_in,cmort,hmort,bmort ) + + ! ============================================================================ + ! Calculate mortality rates as a function of carbon storage + ! ============================================================================ + + + + type (ed_cohort_type), intent(in) :: cohort_in + real(r8),intent(out) :: bmort ! background mortality : Fraction per year + real(r8),intent(out) :: cmort ! carbon starvation mortality + real(r8),intent(out) :: hmort ! hydraulic failure mortality + + real(r8) :: frac ! relativised stored carbohydrate + real(r8) :: b_leaf ! leaf biomass kgC + real(r8) :: hf_sm_threshold ! hydraulic failure soil moisture threshold + + + if (hlm_use_ed_prescribed_phys .eq. ifalse) then + + ! 'Background' mortality (can vary as a function of density as in ED1.0 and ED2.0, but doesn't here for tractability) + bmort = EDPftvarcon_inst%bmort(cohort_in%pft) + + ! Proxy for hydraulic failure induced mortality. + hf_sm_threshold = EDPftvarcon_inst%hf_sm_threshold(cohort_in%pft) + + if(cohort_in%patchptr%btran_ft(cohort_in%pft) <= hf_sm_threshold)then + hmort = ED_val_stress_mort + else + hmort = 0.0_r8 + endif + + ! Carbon Starvation induced mortality. + if ( cohort_in%dbh > 0._r8 ) then + call bleaf(cohort_in%dbh,cohort_in%hite,cohort_in%pft,cohort_in%canopy_trim,b_leaf) + if( b_leaf > 0._r8 .and. cohort_in%bstore <= b_leaf )then + frac = cohort_in%bstore/ b_leaf + cmort = max(0.0_r8,ED_val_stress_mort*(1.0_r8 - frac)) + else + cmort = 0.0_r8 + endif + + else + write(fates_log(),*) 'dbh problem in mortality_rates', & + cohort_in%dbh,cohort_in%pft,cohort_in%n,cohort_in%canopy_layer + endif + + !mortality_rates = bmort + hmort + cmort + + else ! i.e. hlm_use_ed_prescribed_phys is true + if ( cohort_in%canopy_layer .eq. 1) then + bmort = EDPftvarcon_inst%prescribed_mortality_canopy(cohort_in%pft) + else + bmort = EDPftvarcon_inst%prescribed_mortality_understory(cohort_in%pft) + endif + cmort = 0._r8 + hmort = 0._r8 + endif + + end subroutine mortality_rates + +! ============================================================================ + +end module EDMortalityFunctionsMod diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index bde59dafcd..051baf689a 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -31,6 +31,7 @@ module EDPatchDynamicsMod use EDLoggingMortalityMod, only : logging_litter_fluxes use EDLoggingMortalityMod, only : logging_time use EDParamsMod , only : fates_mortality_disturbance_fraction + use FatesAllometryMod , only : carea_allom use FatesConstantsMod , only : g_per_kg use FatesConstantsMod , only : ha_per_m2 use FatesConstantsMod , only : days_per_sec @@ -78,7 +79,7 @@ subroutine disturbance_rates( site_in) ! Modify to add logging disturbance ! !USES: - use EDGrowthFunctionsMod , only : c_area, mortality_rates + use EDMortalityFunctionsMod , only : mortality_rates ! loging flux use EDLoggingMortalityMod , only : LoggingMortality_frac @@ -98,7 +99,7 @@ subroutine disturbance_rates( site_in) real(r8) :: lmort_collateral real(r8) :: lmort_infra - integer :: threshold_sizeclass + integer :: threshold_sizeclass !---------------------------------------------------------------------------------------------- ! Calculate Mortality Rates (these were previously calculated during growth derivatives) @@ -115,7 +116,8 @@ subroutine disturbance_rates( site_in) call mortality_rates(currentCohort,cmort,hmort,bmort) currentCohort%dmort = cmort+hmort+bmort - currentCohort%c_area = c_area(currentCohort) + + call carea_allom(currentCohort%dbh,currentCohort%n,site_in%spread,currentCohort%pft,currentCohort%c_area) ! Initialize diagnostic mortality rates currentCohort%cmort = cmort @@ -779,7 +781,6 @@ subroutine fire_litter_fluxes(currentSite, cp_target, new_patch_target, patch_si ! ! !USES: use SFParamsMod, only : SF_VAL_CWD_FRAC - use EDGrowthFunctionsMod, only : c_area use EDtypesMod , only : dl_sf ! ! !ARGUMENTS: @@ -937,7 +938,7 @@ subroutine fire_litter_fluxes(currentSite, cp_target, new_patch_target, patch_si currentCohort => new_patch%shortest do while(associated(currentCohort)) - currentCohort%c_area = c_area(currentCohort) + call carea_allom(currentCohort%dbh,currentCohort%n,currentSite%spread,currentCohort%pft,currentCohort%c_area) if(EDPftvarcon_inst%woody(currentCohort%pft) == 1)then burned_leaves = (currentCohort%bl+currentCohort%bsw) * currentCohort%cfa else diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 old mode 100755 new mode 100644 index db2f258b93..154de40050 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -17,6 +17,8 @@ module EDPhysiologyMod use FatesInterfaceMod, only : bc_in_type use EDCohortDynamicsMod , only : allocate_live_biomass, zero_cohort use EDCohortDynamicsMod , only : create_cohort, sort_cohorts + use EDCohortDynamicsMod , only : tree_lai + use EDCohortDynamicsMod , only : tree_sai use EDTypesMod , only : numWaterMem use EDTypesMod , only : dl_sf, dinc_ed @@ -33,6 +35,17 @@ module EDPhysiologyMod use EDParamsMod , only : fates_mortality_disturbance_fraction use FatesConstantsMod , only : itrue,ifalse + use FatesAllometryMod , only : h_allom + use FatesAllometryMod , only : h2d_allom + use FatesAllometryMod , only : bag_allom + use FatesAllometryMod , only : bsap_allom + use FatesAllometryMod , only : bleaf + use FatesAllometryMod , only : bfineroot + use FatesAllometryMod , only : bdead_allom + use FatesAllometryMod , only : bcr_allom + use FatesAllometryMod , only : carea_allom + + implicit none private @@ -54,6 +67,8 @@ module EDPhysiologyMod logical, parameter :: DEBUG = .false. ! local debug flag character(len=*), parameter, private :: sourcefile = & __FILE__ + + ! ============================================================================ contains @@ -156,7 +171,6 @@ subroutine trim_canopy( currentSite ) ! ! !USES: ! - use EDGrowthFunctionsMod, only : tree_lai ! ! !ARGUMENTS type (ed_site_type),intent(inout), target :: currentSite @@ -166,7 +180,11 @@ subroutine trim_canopy( currentSite ) type (ed_patch_type) , pointer :: currentPatch integer :: z ! leaf layer + integer :: ipft ! pft index integer :: trimmed ! was this layer trimmed in this year? If not expand the canopy. + real(r8) :: tar_bl ! target leaf biomass (leaves flushed, trimmed) + real(r8) :: tar_bfr ! target fine-root biomass (leaves flushed, trimmed) + real(r8) :: bfr_per_bleaf ! ratio of fine root per leaf biomass !---------------------------------------------------------------------- @@ -176,6 +194,8 @@ subroutine trim_canopy( currentSite ) currentCohort => currentPatch%tallest do while (associated(currentCohort)) trimmed = 0 + ipft = currentCohort%pft + call carea_allom(currentCohort%dbh,currentCohort%n,currentSite%spread,currentCohort%pft,currentCohort%c_area) currentCohort%treelai = tree_lai(currentCohort) currentCohort%nv = ceiling((currentCohort%treelai+currentCohort%treesai)/dinc_ed) if (currentCohort%nv > nlevleaf)then @@ -183,39 +203,44 @@ subroutine trim_canopy( currentSite ) currentCohort%c_area,currentCohort%n,currentCohort%bl endif + call bleaf(currentcohort%dbh,currentcohort%hite,ipft,currentcohort%canopy_trim,tar_bl) + call bfineroot(currentcohort%dbh,currentcohort%hite,ipft,currentcohort%canopy_trim,tar_bfr) + + bfr_per_bleaf = tar_bfr/tar_bl + !Leaf cost vs netuptake for each leaf layer. do z = 1,nlevleaf if (currentCohort%year_net_uptake(z) /= 999._r8)then !there was activity this year in this leaf layer. !Leaf Cost kgC/m2/year-1 !decidous costs. - if (EDPftvarcon_inst%season_decid(currentCohort%pft) == 1.or. & - EDPftvarcon_inst%stress_decid(currentCohort%pft) == 1)then - currentCohort%leaf_cost = 1._r8/(EDPftvarcon_inst%slatop(currentCohort%pft)*1000.0_r8) + if (EDPftvarcon_inst%season_decid(ipft) == 1.or. & + EDPftvarcon_inst%stress_decid(ipft) == 1)then + + + currentCohort%leaf_cost = 1._r8/(EDPftvarcon_inst%slatop(ipft)*1000.0_r8) currentCohort%leaf_cost = currentCohort%leaf_cost + & - 1.0_r8/(EDPftvarcon_inst%slatop(currentCohort%pft)*1000.0_r8) * & - EDPftvarcon_inst%allom_l2fr(currentCohort%pft) / EDPftvarcon_inst%root_long(currentCohort%pft) - currentCohort%leaf_cost = currentCohort%leaf_cost * (EDPftvarcon_inst%grperc(currentCohort%pft) + 1._r8) + 1.0_r8/(EDPftvarcon_inst%slatop(ipft)*1000.0_r8) * bfr_per_bleaf / EDPftvarcon_inst%root_long(ipft) + + currentCohort%leaf_cost = currentCohort%leaf_cost * (EDPftvarcon_inst%grperc(ipft) + 1._r8) else !evergreen costs - currentCohort%leaf_cost = 1.0_r8/(EDPftvarcon_inst%slatop(currentCohort%pft)* & - EDPftvarcon_inst%leaf_long(currentCohort%pft)*1000.0_r8) !convert from sla in m2g-1 to m2kg-1 + currentCohort%leaf_cost = 1.0_r8/(EDPftvarcon_inst%slatop(ipft)* & + EDPftvarcon_inst%leaf_long(ipft)*1000.0_r8) !convert from sla in m2g-1 to m2kg-1 currentCohort%leaf_cost = currentCohort%leaf_cost + & - 1.0_r8/(EDPftvarcon_inst%slatop(currentCohort%pft)*1000.0_r8) * & - EDPftvarcon_inst%allom_l2fr(currentCohort%pft) / EDPftvarcon_inst%root_long(currentCohort%pft) - currentCohort%leaf_cost = currentCohort%leaf_cost * (EDPftvarcon_inst%grperc(currentCohort%pft) + 1._r8) + 1.0_r8/(EDPftvarcon_inst%slatop(ipft)*1000.0_r8) * bfr_per_bleaf / EDPftvarcon_inst%root_long(ipft) + currentCohort%leaf_cost = currentCohort%leaf_cost * (EDPftvarcon_inst%grperc(ipft) + 1._r8) endif if (currentCohort%year_net_uptake(z) < currentCohort%leaf_cost)then - if (currentCohort%canopy_trim > EDPftvarcon_inst%trim_limit(currentCohort%pft))then + if (currentCohort%canopy_trim > EDPftvarcon_inst%trim_limit(ipft))then if ( DEBUG ) then write(fates_log(),*) 'trimming leaves',currentCohort%canopy_trim,currentCohort%leaf_cost endif ! keep trimming until none of the canopy is in negative carbon balance. - if (currentCohort%hite > EDPftvarcon_inst%hgt_min(currentCohort%pft))then - currentCohort%canopy_trim = currentCohort%canopy_trim - EDPftvarcon_inst%trim_inc(currentCohort%pft) - if (EDPftvarcon_inst%evergreen(currentCohort%pft) /= 1)then - currentCohort%laimemory = currentCohort%laimemory*(1.0_r8 - & - EDPftvarcon_inst%trim_inc(currentCohort%pft)) + if (currentCohort%hite > EDPftvarcon_inst%hgt_min(ipft))then + currentCohort%canopy_trim = currentCohort%canopy_trim - EDPftvarcon_inst%trim_inc(ipft) + if (EDPftvarcon_inst%evergreen(ipft) /= 1)then + currentCohort%laimemory = currentCohort%laimemory*(1.0_r8 - EDPftvarcon_inst%trim_inc(ipft)) endif trimmed = 1 endif @@ -231,7 +256,7 @@ subroutine trim_canopy( currentSite ) currentCohort%year_net_uptake(:) = 999.0_r8 if (trimmed == 0.and.currentCohort%canopy_trim < 1.0_r8)then - currentCohort%canopy_trim = currentCohort%canopy_trim + EDPftvarcon_inst%trim_inc(currentCohort%pft) + currentCohort%canopy_trim = currentCohort%canopy_trim + EDPftvarcon_inst%trim_inc(ipft) endif if ( DEBUG ) then @@ -762,7 +787,8 @@ subroutine Growth_Derivatives( currentSite, currentCohort, bc_in) ! Main subroutine controlling growth and allocation derivatives ! ! !USES: - use EDGrowthFunctionsMod , only : Bleaf, dDbhdBd, dhdbd, hite, mortality_rates,dDbhdBl + + use EDMortalityFunctionsMod , only : mortality_rates use FatesInterfaceMod, only : hlm_use_ed_prescribed_phys use EDLoggingMortalityMod, only : LoggingMortality_frac @@ -779,7 +805,7 @@ subroutine Growth_Derivatives( currentSite, currentCohort, bc_in) real(r8) :: dhdbd_fn !rate of change of height per unit dbh real(r8) :: va !fraction of growth going to alive biomass real(r8) :: vs !fraction of growth going to structural biomass - real(r8) :: u,h !intermediates + real(r8) :: u real(r8) :: frac !fraction the stored carbon is of target store amount real(r8) :: f_store !fraction of NPP allocated to storage in this timestep (functionf of stored pool) real(r8) :: gr_fract !fraction of carbon balance that is allocated to growth (not reproduction) @@ -792,14 +818,36 @@ subroutine Growth_Derivatives( currentSite, currentCohort, bc_in) real(r8) :: lmort_collateral ! Mortality fraction associated with logging collateral damage real(r8) :: lmort_infra ! Mortality fraction associated with logging infrastructure real(r8) :: dndt_logging ! Mortality rate (per day) associated with the a logging event - - real(r8) :: balive_loss + real(r8) :: balive_loss ! Carbon that will be removed from the alive pool due to things + ! maintenance turnover + real(r8) :: height ! plant height + + ! Per plant allocation variables + + real(r8) :: b_leaf ! leaf biomass (kgC) + real(r8) :: db_leaf_dd ! change in leaf biomass wrt diameter (kgC/cm) + real(r8) :: b_fineroot ! fine root biomass (kgC) + real(r8) :: db_fineroot_dd ! change in fine root biomass wrt diameter (kgC/cm) + real(r8) :: b_sap ! sapwood biomass (kgC) + real(r8) :: db_sap_dd ! change in sapwood biomass wrt diameter (kgC/cm) + real(r8) :: b_ag ! above ground biomass (kgC/cm) + real(r8) :: db_ag_dd ! change in above ground biomass wrt diameter (kgC/cm) + real(r8) :: b_cr ! coarse root biomass (kgC) + real(r8) :: db_cr_dd ! change in coarse root biomass (kgC/cm) + real(r8) :: b_dead ! dead (structural) biomass (kgC) + real(r8) :: db_dead_dd ! change in dead biomass wrt diameter (kgC/cm) + real(r8) :: dbalivedbd ! change in alive biomass wrt dead biomass (kgC/kgC) + real(r8) :: jh ! plant height (unused) + real(r8) :: dh_dd ! change in plant height WRT diameter (m/cm) + integer :: ipft ! local copy of the pft index !---------------------------------------------------------------------- + ipft = currentCohort%pft + ! Mortality for trees in the understorey. !if trees are in the canopy, then their death is 'disturbance'. This probably needs a different terminology call mortality_rates(currentCohort,cmort,hmort,bmort) - call LoggingMortality_frac(currentCohort%pft, currentCohort%dbh, & + call LoggingMortality_frac(ipft, currentCohort%dbh, & currentCohort%lmort_logging, & currentCohort%lmort_collateral, & currentCohort%lmort_infra ) @@ -818,23 +866,34 @@ subroutine Growth_Derivatives( currentSite, currentCohort, bc_in) endif ! Height - currentCohort%hite = Hite(currentCohort) - h = currentCohort%hite + + call h_allom(currentCohort%dbh,currentCohort%pft,currentCohort%hite) call allocate_live_biomass(currentCohort,0) - ! calculate target size of living biomass compartment for a given dbh. - target_balive = Bleaf(currentCohort) * (1.0_r8 + EDPftvarcon_inst%allom_l2fr(currentCohort%pft) + & - EDpftvarcon_inst%allom_latosa_int(currentCohort%pft)*h) + ! ----------------------------------------------------------------------------------- + ! calculate target size of living biomass compartment for a given dbh. + ! ----------------------------------------------------------------------------------- + + ! Calculate leaf biomass, this wrapper finds the maximum per allometry, and then + ! applies trimming + call bleaf(currentCohort%dbh,currentCohort%hite,ipft,currentCohort%canopy_trim,b_leaf) + + ! Calculate the fine root biomass, this wrapper finds the maximum per allometry, + ! and in the current default case, will trim fine root biomass at the same proportion + ! that it trims leaves + call bfineroot(currentCohort%dbh,currentCohort%hite,ipft,currentCohort%canopy_trim,b_fineroot) + + ! Calculate sapwood biomass + call bsap_allom(currentCohort%dbh,ipft,currentCohort%canopy_trim,b_sap) + + target_balive = b_leaf + b_fineroot + b_sap + !target balive without leaves. if (currentCohort%status_coh == 1)then - target_balive = Bleaf(currentCohort) * (EDPftvarcon_inst%allom_l2fr(currentCohort%pft) + & - EDpftvarcon_inst%allom_latosa_int(currentCohort%pft) * h) + target_balive = b_fineroot + b_sap endif - ! NPP - if ( DEBUG ) write(fates_log(),*) 'EDphys 716 ',currentCohort%npp_acc - ! convert from kgC/indiv/day into kgC/indiv/year ! TODO: CONVERT DAYS_PER_YEAR TO DBLE (HOLDING FOR B4B COMPARISONS, RGK-01-2017) currentCohort%npp_acc_hold = currentCohort%npp_acc * hlm_days_per_year @@ -843,12 +902,12 @@ subroutine Growth_Derivatives( currentSite, currentCohort, bc_in) if (hlm_use_ed_prescribed_phys .eq. itrue) then if (currentCohort%canopy_layer .eq. 1) then - currentCohort%npp_acc_hold = EDPftvarcon_inst%prescribed_npp_canopy(currentCohort%pft) & - * currentCohort%c_area / currentCohort%n + currentCohort%npp_acc_hold = EDPftvarcon_inst%prescribed_npp_canopy(ipft) & + * currentCohort%c_area / currentCohort%n currentCohort%npp_acc = currentCohort%npp_acc_hold / hlm_days_per_year ! add these for balance checking purposes else - currentCohort%npp_acc_hold = EDPftvarcon_inst%prescribed_npp_understory(currentCohort%pft) & - * currentCohort%c_area / currentCohort%n + currentCohort%npp_acc_hold = EDPftvarcon_inst%prescribed_npp_understory(ipft) & + * currentCohort%c_area / currentCohort%n currentCohort%npp_acc = currentCohort%npp_acc_hold / hlm_days_per_year ! add these for balance checking purposes endif endif @@ -856,9 +915,9 @@ subroutine Growth_Derivatives( currentSite, currentCohort, bc_in) currentSite%flux_in = currentSite%flux_in + currentCohort%npp_acc * currentCohort%n ! Maintenance demands - if (EDPftvarcon_inst%evergreen(currentCohort%pft) == 1)then !grass and EBT - currentCohort%leaf_md = currentCohort%bl / EDPftvarcon_inst%leaf_long(currentCohort%pft) - currentCohort%root_md = currentCohort%br / EDPftvarcon_inst%root_long(currentCohort%pft) + if (EDPftvarcon_inst%evergreen(ipft) == 1)then !grass and EBT + currentCohort%leaf_md = currentCohort%bl / EDPftvarcon_inst%leaf_long(ipft) + currentCohort%root_md = currentCohort%br / EDPftvarcon_inst%root_long(ipft) currentCohort%md = currentCohort%root_md + currentCohort%leaf_md endif @@ -867,24 +926,24 @@ subroutine Growth_Derivatives( currentSite, currentCohort, bc_in) !with which I am not especially comfortable, particularly as the concept of sapwood turnover is unclear for trees that !are still in an expansion phase. - if (EDPftvarcon_inst%season_decid(currentCohort%pft) == 1)then - currentCohort%root_md = currentCohort%br /EDPftvarcon_inst%root_long(currentCohort%pft) + if (EDPftvarcon_inst%season_decid(ipft) == 1)then + currentCohort%root_md = currentCohort%br /EDPftvarcon_inst%root_long(ipft) currentCohort%leaf_md = 0._r8 currentCohort%md = currentCohort%root_md + currentCohort%leaf_md endif - if (EDPftvarcon_inst%stress_decid(currentCohort%pft) == 1)then - currentCohort%root_md = currentCohort%br /EDPftvarcon_inst%root_long(currentCohort%pft) + if (EDPftvarcon_inst%stress_decid(ipft) == 1)then + currentCohort%root_md = currentCohort%br /EDPftvarcon_inst%root_long(ipft) currentCohort%leaf_md = 0._r8 currentCohort%md = currentCohort%root_md + currentCohort%leaf_md endif - if (EDPftvarcon_inst%stress_decid(currentCohort%pft) /= 1 & - .and.EDPftvarcon_inst%season_decid(currentCohort%pft) /= 1.and. & - EDPftvarcon_inst%evergreen(currentCohort%pft) /= 1)then - write(fates_log(),*) 'problem with phenology definitions',currentCohort%pft, & - EDPftvarcon_inst%stress_decid(currentCohort%pft), & - EDPftvarcon_inst%season_decid(currentCohort%pft),EDPftvarcon_inst%evergreen(currentCohort%pft) + if (EDPftvarcon_inst%stress_decid(ipft) /= 1 & + .and.EDPftvarcon_inst%season_decid(ipft) /= 1.and. & + EDPftvarcon_inst%evergreen(ipft) /= 1)then + write(fates_log(),*) 'problem with phenology definitions',ipft, & + EDPftvarcon_inst%stress_decid(ipft), & + EDPftvarcon_inst%season_decid(ipft),EDPftvarcon_inst%evergreen(ipft) endif ! FIX(RF,032414) -turned off for now as it makes balive go negative.... @@ -896,26 +955,26 @@ subroutine Growth_Derivatives( currentSite, currentCohort, bc_in) ! this is the fraction of maintenance demand we -have- to do... if ( DEBUG ) write(fates_log(),*) 'EDphys 760 ',currentCohort%npp_acc_hold, currentCohort%md, & - EDPftvarcon_inst%leaf_stor_priority(currentCohort%pft) + EDPftvarcon_inst%leaf_stor_priority(ipft) currentCohort%carbon_balance = currentCohort%npp_acc_hold - & - currentCohort%md * EDPftvarcon_inst%leaf_stor_priority(currentCohort%pft) + currentCohort%md * EDPftvarcon_inst%leaf_stor_priority(ipft) ! Allowing only carbon from NPP pool to account for npp flux into the maintenance turnover pools ! ie this does not include any use of storage carbon or balive to make up for missing carbon balance in the transfer currentCohort%npp_leaf = max(0.0_r8,min(currentCohort%npp_acc_hold*currentCohort%leaf_md/currentCohort%md, & - currentCohort%leaf_md*EDPftvarcon_inst%leaf_stor_priority(currentCohort%pft))) + currentCohort%leaf_md*EDPftvarcon_inst%leaf_stor_priority(ipft))) currentCohort%npp_froot = max(0.0_r8,min(currentCohort%npp_acc_hold*currentCohort%root_md/currentCohort%md, & - currentCohort%root_md*EDPftvarcon_inst%leaf_stor_priority(currentCohort%pft))) + currentCohort%root_md*EDPftvarcon_inst%leaf_stor_priority(ipft))) - if (Bleaf(currentCohort) > 0._r8)then + if (b_leaf > 0._r8)then if ( DEBUG ) write(fates_log(),*) 'EDphys A ',currentCohort%carbon_balance if (currentCohort%carbon_balance > 0._r8)then !spend C on growing and storing !what fraction of the target storage do we have? - frac = max(0.0_r8,currentCohort%bstore/(Bleaf(currentCohort) * EDPftvarcon_inst%cushion(currentCohort%pft))) + frac = max(0.0_r8,currentCohort%bstore/( b_leaf * EDPftvarcon_inst%cushion(ipft))) ! FIX(SPM,080514,fstore never used ) f_store = max(exp(-1.*frac**4._r8) - exp( -1.0_r8 ),0.0_r8) !what fraction of allocation do we divert to storage? @@ -940,21 +999,21 @@ subroutine Growth_Derivatives( currentSite, currentCohort, bc_in) else - write(fates_log(),*) 'No target leaf area in GrowthDerivs? Bleaf(cohort) <= 0?' + write(fates_log(),*) 'No target leaf area in GrowthDerivs? b_leaf <= 0?' call endrun(msg=errMsg(sourcefile, __LINE__)) endif !Do we have enough carbon left over to make up the rest of the turnover demand? balive_loss = 0._r8 - if (currentCohort%carbon_balance > currentCohort%md*(1.0_r8- EDPftvarcon_inst%leaf_stor_priority(currentCohort%pft)))then ! Yes... + if (currentCohort%carbon_balance > currentCohort%md*(1.0_r8- EDPftvarcon_inst%leaf_stor_priority(ipft)))then ! Yes... currentCohort%carbon_balance = currentCohort%carbon_balance - currentCohort%md * (1.0_r8 - & - EDPftvarcon_inst%leaf_stor_priority(currentCohort%pft)) + EDPftvarcon_inst%leaf_stor_priority(ipft)) currentCohort%npp_leaf = currentCohort%npp_leaf + & - currentCohort%leaf_md * (1.0_r8-EDPftvarcon_inst%leaf_stor_priority(currentCohort%pft)) + currentCohort%leaf_md * (1.0_r8-EDPftvarcon_inst%leaf_stor_priority(ipft)) currentCohort%npp_froot = currentCohort%npp_froot + & - currentCohort%root_md * (1.0_r8-EDPftvarcon_inst%leaf_stor_priority(currentCohort%pft)) + currentCohort%root_md * (1.0_r8-EDPftvarcon_inst%leaf_stor_priority(ipft)) else ! we can't maintain constant leaf area and root area. Balive is reduced @@ -963,8 +1022,7 @@ subroutine Growth_Derivatives( currentSite, currentCohort, bc_in) currentCohort%npp_froot = currentCohort%npp_froot + & max(0.0_r8,currentCohort%carbon_balance*(currentCohort%root_md/currentCohort%md)) - balive_loss = currentCohort%md *(1.0_r8- EDPftvarcon_inst%leaf_stor_priority(currentCohort%pft)) - & - currentCohort%carbon_balance + balive_loss = currentCohort%md *(1.0_r8- EDPftvarcon_inst%leaf_stor_priority(ipft))- currentCohort%carbon_balance currentCohort%carbon_balance = 0._r8 endif @@ -973,38 +1031,81 @@ subroutine Growth_Derivatives( currentSite, currentCohort, bc_in) !********************************************/ !Use remaining carbon to refill balive or to get larger. + ! Tally up the relative change in dead biomass WRT diameter + call bag_allom(currentCohort%dbh,currentCohort%hite,ipft,b_ag,db_ag_dd) + call bcr_allom(currentCohort%dbh,currentCohort%hite,ipft,b_cr,db_cr_dd) + call bdead_allom( b_ag, b_cr, b_sap, ipft, b_dead, db_ag_dd, db_cr_dd, db_sap_dd, db_dead_dd ) + !only if carbon balance is +ve - if ((currentCohort%balive >= target_balive).AND.(currentCohort%carbon_balance > 0._r8))then + + if ((currentCohort%balive >= target_balive).and.(currentCohort%carbon_balance > 0._r8))then ! fraction of carbon going into active vs structural carbon - if (currentCohort%dbh <= EDPftvarcon_inst%allom_dbh_maxheight(currentCohort%pft))then ! cap on leaf biomass - dbldbd = dDbhdBd(currentCohort)/dDbhdBl(currentCohort) - dbrdbd = EDPftvarcon_inst%allom_l2fr(currentCohort%pft) * dbldbd - dhdbd_fn = dhdbd(currentCohort) - dbswdbd = EDpftvarcon_inst%allom_latosa_int(currentCohort%pft) * (h*dbldbd + currentCohort%bl*dhdbd_fn) - u = 1.0_r8 / (dbldbd + dbrdbd + dbswdbd) + + ! fraction of carbon not going towards reproduction + if (currentCohort%dbh <= EDPftvarcon_inst%dbh_repro_threshold(ipft)) then ! cap on leaf biomass + gr_fract = 1.0_r8 - EDPftvarcon_inst%seed_alloc(ipft) + else + gr_fract = 1.0_r8 - (EDPftvarcon_inst%seed_alloc(ipft) + EDPftvarcon_inst%clone_alloc(ipft)) + end if + + ! Tally up the relative change in alive biomass WRT diameter + ! These calculations will take into account any height capping + ! (if the user wanted it) and its implications to these pools + call bleaf(currentCohort%dbh,currentCohort%hite,ipft, & + currentCohort%canopy_trim,b_leaf,db_leaf_dd) + call bfineroot(currentCohort%dbh,currentCohort%hite,ipft, & + currentCohort%canopy_trim,b_fineroot,db_fineroot_dd) + call bsap_allom(currentCohort%dbh,ipft, & + currentCohort%canopy_trim,b_sap,db_sap_dd) + + ! Total change in alive biomass relative to dead biomass [kgC/kgC] + dbalivedbd = (db_leaf_dd + db_fineroot_dd + db_sap_dd)/db_dead_dd + + if(dbalivedbd>tiny(dbalivedbd))then + + ! In this case, the plant allometry module is telling us that + ! the plant is still expected to gain live (leaf,froot,sap) + ! biomass as it grows in size, and therfore it should be + ! allocated in proportion with structural + + u = 1.0_r8 / dbalivedbd va = 1.0_r8 / (1.0_r8 + u) vs = u / (1.0_r8 + u) - gr_fract = 1.0_r8 - EDPftvarcon_inst%seed_alloc(currentCohort%pft) + else - dbldbd = 0._r8; dbrdbd = 0._r8 ;dbswdbd = 0._r8 + + ! If there is no change in alive biomass per change in dead, + ! most likely we are dealing with plants that have surpassed + ! an allometric threshold that limits alive biomass. + va = 0.0_r8 vs = 1.0_r8 - gr_fract = 1.0_r8 - (EDPftvarcon_inst%seed_alloc(currentCohort%pft) + EDPftvarcon_inst%clone_alloc(currentCohort%pft)) - endif + + end if - !FIX(RF,032414) - to fix high bl's. needed to prevent numerical errors without the ODEINT. + ! FIX(RF,032414) - to fix high bl's. needed to + ! prevent numerical errors without the ODEINT. if (currentCohort%balive > target_balive*1.1_r8)then va = 0.0_r8; vs = 1._r8 - if (DEBUG) write(fates_log(),*) 'using high bl cap',target_balive,currentCohort%balive + if (DEBUG) write(fates_log(),*) 'using high bl cap',target_balive,currentCohort%balive endif else - dbldbd = 0._r8; dbrdbd = 0._r8; dbswdbd = 0._r8 - va = 1.0_r8; vs = 0._r8 + + ! -------------------------------------------------------------------------------- + ! In this case, either there was not enough carbon generated, or the plant is off + ! allometry (ie the alive pools are smaller than the maximums dictated by allometry + ! and timming). So push all carbon into the alive pool (va = 1.0), and none + ! into structural (vs = 0.0) or seed (ie non-growth, gr_fract = 1.0). + ! -------------------------------------------------------------------------------- + + va = 1.0_r8 + vs = 0.0_r8 gr_fract = 1.0_r8 - endif - ! calculate derivatives of living and dead carbon pools + endif + + ! calculate derivatives of living and dead carbon pools currentCohort%dbalivedt = gr_fract * va * currentCohort%carbon_balance - balive_loss currentCohort%dbdeaddt = gr_fract * vs * currentCohort%carbon_balance currentCohort%dbstoredt = currentCohort%storage_flux @@ -1038,14 +1139,17 @@ subroutine Growth_Derivatives( currentSite, currentCohort, bc_in) currentCohort%npp_bseed = currentCohort%seed_prod - ! calculate change in diameter and height - currentCohort%ddbhdt = currentCohort%dbdeaddt * dDbhdBd(currentCohort) - currentCohort%dhdt = currentCohort%dbdeaddt * dHdBd(currentCohort) + ! calculate change in diameter + currentCohort%ddbhdt = currentCohort%dbdeaddt / db_dead_dd + + ! calculate change in hite + call h_allom(currentCohort%dbh,ipft,height,dh_dd) + currentCohort%dhdt = currentCohort%ddbhdt * dh_dd ! If the cohort has grown, it is not new currentCohort%isnew=.false. - end subroutine Growth_Derivatives + end subroutine Growth_Derivatives ! ============================================================================ subroutine recruitment( currentSite, currentPatch, bc_in ) @@ -1054,7 +1158,6 @@ subroutine recruitment( currentSite, currentPatch, bc_in ) ! spawn new cohorts of juveniles of each PFT ! ! !USES: - use EDGrowthFunctionsMod, only : bdead,dbh, Bleaf use FatesInterfaceMod, only : hlm_use_ed_prescribed_phys ! ! !ARGUMENTS @@ -1066,6 +1169,12 @@ subroutine recruitment( currentSite, currentPatch, bc_in ) integer :: ft type (ed_cohort_type) , pointer :: temp_cohort integer :: cohortstatus + real(r8) :: b_leaf + real(r8) :: b_fineroot + real(r8) :: b_sapwood + real(r8) :: b_aboveground + real(r8) :: b_coarseroot + !---------------------------------------------------------------------- allocate(temp_cohort) ! create temporary cohort @@ -1076,12 +1185,19 @@ subroutine recruitment( currentSite, currentPatch, bc_in ) temp_cohort%canopy_trim = 0.8_r8 !starting with the canopy not fully expanded temp_cohort%pft = ft temp_cohort%hite = EDPftvarcon_inst%hgt_min(ft) - temp_cohort%dbh = Dbh(temp_cohort) - temp_cohort%bdead = Bdead(temp_cohort) - temp_cohort%balive = Bleaf(temp_cohort)*(1.0_r8 + EDPftvarcon_inst%allom_l2fr(ft) & - + EDpftvarcon_inst%allom_latosa_int(ft)*temp_cohort%hite) - temp_cohort%bstore = EDPftvarcon_inst%cushion(ft)*(temp_cohort%balive/ (1.0_r8 + EDPftvarcon_inst%allom_l2fr(ft) & - + EDpftvarcon_inst%allom_latosa_int(ft)*temp_cohort%hite)) + call h2d_allom(temp_cohort%hite,ft,temp_cohort%dbh) + + ! Initialize balive (leaf+fineroot+sapwood) + call bleaf(temp_cohort%dbh,temp_cohort%hite,ft,temp_cohort%canopy_trim,b_leaf) + call bfineroot(temp_cohort%dbh,temp_cohort%hite,ft,temp_cohort%canopy_trim,b_fineroot) + call bsap_allom(temp_cohort%dbh,ft,temp_cohort%canopy_trim,b_sapwood) + + call bag_allom(temp_cohort%dbh,temp_cohort%hite,ft,b_aboveground) + call bcr_allom(temp_cohort%dbh,temp_cohort%hite,ft,b_coarseroot) + call bdead_allom(b_aboveground,b_coarseroot,b_sapwood,ft,temp_cohort%bdead) + + temp_cohort%balive = b_leaf + b_sapwood + b_fineroot + temp_cohort%bstore = EDPftvarcon_inst%cushion(ft) * b_leaf if (hlm_use_ed_prescribed_phys .eq. ifalse .or. EDPftvarcon_inst%prescribed_recruitment(ft) .lt. 0. ) then temp_cohort%n = currentPatch%area * currentPatch%seed_germination(ft)*hlm_freq_day & @@ -1098,12 +1214,10 @@ subroutine recruitment( currentSite, currentPatch, bc_in ) temp_cohort%laimemory = 0.0_r8 if (EDPftvarcon_inst%season_decid(temp_cohort%pft) == 1.and.currentSite%status == 1)then - temp_cohort%laimemory = (1.0_r8/(1.0_r8 + EDPftvarcon_inst%allom_l2fr(ft) + & - EDpftvarcon_inst%allom_latosa_int(ft)*temp_cohort%hite))*temp_cohort%balive + temp_cohort%laimemory = b_leaf endif if (EDPftvarcon_inst%stress_decid(temp_cohort%pft) == 1.and.currentSite%dstatus == 1)then - temp_cohort%laimemory = (1.0_r8/(1.0_r8 + EDPftvarcon_inst%allom_l2fr(ft) + & - EDpftvarcon_inst%allom_latosa_int(ft)*temp_cohort%hite))*temp_cohort%balive + temp_cohort%laimemory = b_leaf endif cohortstatus = currentSite%status diff --git a/biogeochem/FatesAllometryMod.F90 b/biogeochem/FatesAllometryMod.F90 new file mode 100644 index 0000000000..38e0631124 --- /dev/null +++ b/biogeochem/FatesAllometryMod.F90 @@ -0,0 +1,1465 @@ +!=============================================================================== +! +! FatesAllometryMod.F90 +! +! A library of functions that calculate plant allometry and their +! derivatives. Most relationships are related to diameter [cm]. All +! derivatives with respect to change in diameter have same units. +! In some cases multiple areguments will be provided, yet +! those arguments in those cases are trivially determined from diameter. +! +! Each function presented in this library is written to also return the +! derivative with respect to diameter, if the argument is provided. +! +! The name convention of the functions follows the form d... Which +! indicates "diameter to ...". Allometries for the following variables are +! calculated: +! h: height [m] +! bag: biomass above ground [kgC] (aka AGB) +! blmax: biomass in leaves when leaves are "on allometry" +! this also is the "pre-trimmed" value, which is the maximum +! or potential leaf mass a tree may have [kgC] +! bcr: biomass in coarse roots [kgC] (belowground sap+dead wood, no fines) +! bfrmax: biomass in fine roots when "on allometry" [kgC] +! bsap: biomass in sapwood (above and below) [kgC] +! bdead: biomass (above and below) in the structural pool [kgC] +! +! +! The following function switches are rused +! +! allom_hmode, integer, height allometry function type +! allom_lmode, integer, maximum leaf allometry function type +! allom_rmode, integer, maximum root allometry function type +! allom_amode, integer, AGB allometry function type +! allom_cmode, integer, coarse root allometry function type +! allom_smode, integer, sapwood allometry function type +! +! The following parameters (traits) are used +! +! wood_density, mean stem wood specific gravity (heart,sap,bark) +! allom_latosa_int, leaf area to sap area ratio, intercept [m2/cm2] +! allom_latosa_slp, leaf area to sap area ratio, slope on diameter [m2/cm2/cm] +! c2b, real, carbon to biomass multiplier (~2.0) +! allom_l2fr, fine root biomass per leaf biomass ratio [kgC/kgC] +! allom_agb_frac, the fraction of stem above ground [-] +! allom_d2h1, parameter 1 for d2h allometry (intercept) +! allom_d2h2, parameter 2 for d2h allometry (slope) +! allom_d2h3, parameter 3 for d2h allometry (optional) +! allom_d2bl1, parameter 1 for d2bl allometry (intercept) +! allom_d2bl2, parameter 2 for d2bl allometry (slope) +! allom_d2bl3, parameter 3 for d2bl allometry (optional) +! allom_agb1 +! allom_agb2 +! allom_agb3 +! allom_dbh_maxheight, dbh at maximum height [cm] +! h_min, the height associated with newly recruited plant [m] +! +! Note - i4 types are expressed explicitly to accomodate unit testing calls +! to this module +! +! +! OPEN QUESTIONS: +! SHOULD SAPWOOD ALLOMETRY BE EFFECTED BY TRIMMING, OR BE OFF OF BLMAX? +! WHAT POOLS DO WE ASSUME ARE CONTAINED in AGB ALLOMETRY? +! WE ARE NOT EXTENDING SAPWOOD BELOW GROUND, WE PROBABLY SHOULD +! +! +! Carbon Pool Configurations are as follows, and assume a constant proportionality +! between above and below-ground pools. Sapwood (bsap) is both above and below +! ground. Above ground biomass contains above ground dead wood, above ground +! sapwood and leaves. Coarse roots (bcr) contain only the below ground +! deadwood, not the fine roots or the sapwood. Coarse roots are typically +! a proportion of above ground deadwood. +! Leaf biomass, height and above ground biomass typically have non-linear +! allometry models. The default for sapwood is the pipe model. +! +! We ignore leaf biomass contributions in allometry to AGB: +! +! bag = (bdead+bsap)*agb_frac +! bdead = bag - (bsap*agb_frac) + bcr +! bcr = bdead * (1-agb_frac) = (bag - (bsap*agb_frac) + bcr)*(1-agb_frac) +! +! +! Initial Implementation: Ryan Knox July 2017 +! +!=============================================================================== + +module FatesAllometryMod + + ! If this is a unit-test, these globals will be provided by a wrapper + + use EDPFTvarcon , only : EDPftvarcon_inst + use FatesConstantsMod, only : r8 => fates_r8 + use FatesConstantsMod, only : i4 => fates_int + use shr_log_mod , only : errMsg => shr_log_errMsg + use FatesGlobals , only : fates_log + use FatesGlobals , only : endrun => fates_endrun + + implicit none + + private + public :: h2d_allom ! Generic height to diameter wrapper + public :: h_allom ! Generic diameter to height wrapper + public :: bag_allom ! Generic AGB wrapper + public :: blmax_allom ! Generic maximum leaf biomass wrapper + public :: bleaf ! Generic actual leaf biomass wrapper + public :: bsap_allom ! Generic sapwood wrapper + public :: bcr_allom ! Generic coarse root wrapper + public :: bfineroot ! Generic actual fine root biomass wrapper + public :: bdead_allom ! Generic bdead wrapper + public :: carea_allom ! Generic crown area wrapper + + character(len=*), parameter :: sourcefile = __FILE__ + + ! If testing b4b with older versions, do not remove sapwood + ! Our old methods with saldarriaga did not remove sapwood from the + ! bdead pool. But newer allometries are providing total agb + ! which includes sapwood. Although a small quantity, it needs to be removed + ! from the agb pool. + ! Additionally, our calculation of sapwood biomass may be missing some unite conversions + ! + logical,parameter :: test_b4b = .true. + +contains + + ! ============================================================================ + ! Parameter Checks + ! ============================================================================ + + ! Checks to make sure parameters are not within expected ranges for each + ! functions + + ! Check to make sure Martinez-Cano height cap is not on, or explicitly allowed + + + ! ============================================================================ + ! Generic height to diameter interface + ! ============================================================================ + + subroutine h2d_allom(h,ipft,d,dddh) + + real(r8),intent(in) :: h ! height of plant [m] + integer(i4),intent(in) :: ipft ! PFT index + real(r8),intent(out) :: d ! plant diameter [cm] + real(r8),intent(out),optional :: dddh ! change in diameter per height [cm/m] + + associate( p1 => EDPftvarcon_inst%allom_d2h1(ipft), & + p2 => EDPftvarcon_inst%allom_d2h2(ipft), & + p3 => EDPftvarcon_inst%allom_d2h3(ipft), & + allom_hmode => EDPftvarcon_inst%allom_hmode(ipft)) + + select case(int(allom_hmode)) + case (1) ! Obrien et al. 199X BCI + call h2d_obrien(h,p1,p2,d,dddh) + case (2) ! poorter 2006 + call h2d_poorter2006(h,p1,p2,p3,d,dddh) + case (3) ! 2 parameter power function + call h2d_2pwr(h,p1,p2,d,dddh) + case (4) ! chave 2014 + call h2d_chave2014(h,p1,p2,p3,d,dddh) + case (5) ! Martinez-Cano + call h2d_martcano(h,p1,p2,p3,d,dddh) + case DEFAULT + write(fates_log(),*) 'An undefined h2d allometry was specified: ',allom_hmode + write(fates_log(),*) 'Aborting' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end select + + end associate + return + end subroutine h2d_allom + + ! ============================================================================ + ! Generic height interface + ! ============================================================================ + + subroutine h_allom(d,ipft,h,dhdd) + + ! Arguments + real(r8),intent(in) :: d ! plant diameter [cm] + integer(i4),intent(in) :: ipft ! PFT index + real(r8),intent(out) :: h ! plant height [m] + real(r8),intent(out),optional :: dhdd ! change in height per diameter [m/cm] + + associate( dbh_maxh => EDPftvarcon_inst%allom_dbh_maxheight(ipft), & + p1 => EDPftvarcon_inst%allom_d2h1(ipft), & + p2 => EDPftvarcon_inst%allom_d2h2(ipft), & + p3 => EDPftvarcon_inst%allom_d2h3(ipft), & + allom_hmode => EDPftvarcon_inst%allom_hmode(ipft)) + + select case(int(allom_hmode)) + case (1) ! "obrien" + call d2h_obrien(d,p1,p2,dbh_maxh,h,dhdd) + case (2) ! "poorter06" + call d2h_poorter2006(d,p1,p2,p3,dbh_maxh,h,dhdd) + case (3) ! "2parameter power function h=a*d^b " + call d2h_2pwr(d,p1,p2,dbh_maxh,h,dhdd) + case (4) ! "chave14") + call d2h_chave2014(d,p1,p2,p3,dbh_maxh,h,dhdd) + case (5) ! Martinez-Cano + call d2h_martcano(d,p1,p2,p3,dbh_maxh,h,dhdd) + case DEFAULT + write(fates_log(),*) 'An undefined height allometry was specified: ',allom_hmode + write(fates_log(),*) 'Aborting' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end select + + end associate + return + end subroutine h_allom + + ! ============================================================================ + ! Generic AGB interface + ! ============================================================================ + + subroutine bag_allom(d,h,ipft,bag,dbagdd) + + + real(r8),intent(in) :: d ! plant diameter [cm] + real(r8),intent(in) :: h ! plant height [m] + integer(i4),intent(in) :: ipft ! PFT index + real(r8),intent(out) :: bag ! plant height [m] + real(r8),intent(out),optional :: dbagdd ! change in agb per diameter [kgC/cm] + + real(r8) :: hj ! height (dummy arg) + real(r8) :: dhdd ! change in height wrt d + + associate( p1 => EDPftvarcon_inst%allom_agb1(ipft), & + p2 => EDPftvarcon_inst%allom_agb2(ipft), & + p3 => EDPftvarcon_inst%allom_agb3(ipft), & + p4 => EDPftvarcon_inst%allom_agb4(ipft), & + wood_density => EDPftvarcon_inst%wood_density(ipft), & + c2b => EDPftvarcon_inst%c2b(ipft), & + agb_frac => EDPftvarcon_inst%allom_agb_frac(ipft), & + allom_amode => EDPftvarcon_inst%allom_amode(ipft)) + + select case(int(allom_amode)) + case (1) !"salda") + call h_allom(d,ipft,hj,dhdd) + call dh2bag_salda(d,h,dhdd,p1,p2,p3,p4,wood_density,c2b,agb_frac,bag,dbagdd) + case (2) !"2par_pwr") + ! Switch for woodland dbh->drc + call d2bag_2pwr(d,p1,p2,c2b,bag,dbagdd) + case (3) !"chave14") + call h_allom(d,ipft,hj,dhdd) + call dh2bag_chave2014(d,h,dhdd,p1,p2,wood_density,c2b,bag,dbagdd) + case DEFAULT + write(fates_log(),*) 'An undefined AGB allometry was specified: ',allom_amode + write(fates_log(),*) 'Aborting' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end select + + end associate + return + end subroutine bag_allom + + ! ============================================================================ + ! Generic diameter to maximum leaf biomass interface + ! ============================================================================ + + subroutine blmax_allom(d,h,ipft,blmax,dblmaxdd) + + real(r8),intent(in) :: d ! plant diameter [cm] + real(r8),intent(in) :: h ! plant height [m] + integer(i4),intent(in) :: ipft ! PFT index + real(r8),intent(out) :: blmax ! plant leaf biomass [kg] + real(r8),intent(out),optional :: dblmaxdd ! change leaf bio per diameter [kgC/cm] + + associate( dbh_maxh => EDPftvarcon_inst%allom_dbh_maxheight(ipft), & + rho => EDPftvarcon_inst%wood_density(ipft), & + c2b => EDPftvarcon_inst%c2b(ipft), & + allom_lmode => EDPftvarcon_inst%allom_lmode(ipft), & + p1 => EDPftvarcon_inst%allom_d2bl1(ipft), & + p2 => EDPftvarcon_inst%allom_d2bl2(ipft), & + p3 => EDPftvarcon_inst%allom_d2bl3(ipft)) + + select case(int(allom_lmode)) + case(1) !"salda") + call d2blmax_salda(d,p1,p2,p3,rho,dbh_maxh,c2b,blmax,dblmaxdd) + case(2) !"2par_pwr") + call d2blmax_2pwr(d,p1,p2,c2b,blmax,dblmaxdd) + case(3) ! dh2blmax_2pwr + call dh2blmax_2pwr(d,p1,p2,dbh_maxh,c2b,blmax,dblmaxdd) + case DEFAULT + write(fates_log(),*) 'An undefined leaf allometry was specified: ', & + allom_lmode + write(fates_log(),*) 'Aborting' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end select + end associate + return + end subroutine blmax_allom + + ! ============================================================================ + ! Generic crown area allometry wrapper + ! ============================================================================ + + subroutine carea_allom(d,nplant,site_spread,ipft,c_area) + + real(r8),intent(in) :: d ! plant diameter [cm] + real(r8),intent(in) :: site_spread ! site level spread factor (crowdedness) + real(r8),intent(in) :: nplant ! number of plants [1/ha] + integer(i4),intent(in) :: ipft ! PFT index + real(r8),intent(out) :: c_area ! crown area per plant (m2) + + real(r8) :: d_eff ! Effective diameter (cm) + + associate( dbh_maxh => EDPftvarcon_inst%allom_dbh_maxheight(ipft), & + allom_lmode => EDPftvarcon_inst%allom_lmode(ipft), & + d2bl_p2 => EDPftvarcon_inst%allom_d2bl2(ipft), & + d2bl_ediff => EDPftvarcon_inst%allom_blca_expnt_diff(ipft), & + d2ca_min => EDPftvarcon_inst%allom_d2ca_coefficient_min(ipft), & + d2ca_max => EDPftvarcon_inst%allom_d2ca_coefficient_max(ipft)) + + select case(int(allom_lmode)) + case(1,3) ! "salda" and "height capped generic two power" + d_eff = min(d,dbh_maxh) + call carea_2pwr(d_eff,site_spread,d2bl_p2,d2bl_ediff,d2ca_min,d2ca_max,c_area) + case(2) ! "2par_pwr") + call carea_2pwr(d,site_spread,d2bl_p2,d2bl_ediff,d2ca_min,d2ca_max,c_area) + case DEFAULT + write(fates_log(),*) 'An undefined leaf allometry was specified: ', & + allom_lmode + write(fates_log(),*) 'Aborting' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end select + + c_area = c_area * nplant + + + end associate + return + end subroutine carea_allom + + + + ! ===================================================================================== + + subroutine bleaf(d,h,ipft,canopy_trim,bl,dbldd) + + ! ------------------------------------------------------------------------- + ! This subroutine calculates the actual target bleaf + ! based on trimming. Because trimming + ! is not allometry and rather an emergent property, + ! this routine is not name-spaced with allom_ + ! ------------------------------------------------------------------------- + + real(r8),intent(in) :: d ! plant diameter [cm] + real(r8),intent(in) :: h ! plant height [m] + integer(i4),intent(in) :: ipft ! PFT index + real(r8),intent(in) :: canopy_trim ! trimming function + real(r8),intent(out) :: bl ! plant leaf biomass [kg] + real(r8),intent(out),optional :: dbldd ! change leaf bio per diameter [kgC/cm] + + real(r8) :: blmax + real(r8) :: dblmaxdd + + call blmax_allom(d,h,ipft,blmax,dblmaxdd) + + ! ------------------------------------------------------------------------- + ! Adjust for canopies that have become so deep that their bottom layer is + ! not producing any carbon... + ! nb this will change the allometry and the effects of this remain untested + ! RF. April 2014 + ! ------------------------------------------------------------------------- + + bl = blmax * canopy_trim + + if(present(dbldd))then + dbldd = dblmaxdd * canopy_trim + end if + + return + end subroutine bleaf + + + ! ============================================================================ + ! Generic sapwood biomass interface + ! ============================================================================ + + subroutine bsap_allom(d,ipft,canopy_trim,bsap,dbsapdd) + + real(r8),intent(in) :: d ! plant diameter [cm] + integer(i4),intent(in) :: ipft ! PFT index + real(r8),intent(in) :: canopy_trim + real(r8),intent(out) :: bsap ! plant leaf biomass [kgC] + real(r8),intent(out),optional :: dbsapdd ! change leaf bio per d [kgC/cm] + + real(r8) :: h ! Plant height [m] + real(r8) :: dhdd + real(r8) :: blmax + real(r8) :: dblmaxdd + real(r8) :: bag + real(r8) :: dbagdd + + select case(int(EDPftvarcon_inst%allom_smode(ipft))) + ! --------------------------------------------------------------------- + ! Currently both sapwood area proportionality methods use the same + ! machinery. The only differences are related to the parameter + ! checking at the beginning. For constant proportionality, the slope + ! of the la:sa to diameter line is zero. + ! --------------------------------------------------------------------- + case(1,2) !"constant","dlinear") + + call h_allom(d,ipft,h,dhdd) + if(test_b4b)then + call bleaf(d,h,ipft,canopy_trim,blmax,dblmaxdd) + else + call blmax_allom(d,h,ipft,blmax,dblmaxdd) + end if + call bag_allom(d,h,ipft,bag,dbagdd) + call bsap_dlinear(d,h,dhdd,blmax,dblmaxdd,bag,dbagdd,ipft,bsap,dbsapdd) + case DEFAULT + write(fates_log(),*) 'An undefined sapwood allometry was specified: ', & + EDPftvarcon_inst%allom_smode(ipft) + write(fates_log(),*) 'Aborting' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end select + return + end subroutine bsap_allom + + ! ============================================================================ + ! Generic coarse root biomass interface + ! ============================================================================ + + subroutine bcr_allom(d,h,ipft,bcr,dbcrdd) + + + real(r8),intent(in) :: d ! plant diameter [cm] + real(r8),intent(in) :: h ! plant height [m] + integer(i4),intent(in) :: ipft ! PFT index + real(r8),intent(out) :: bcr ! coarse root biomass [kgC] + real(r8),intent(out),optional :: dbcrdd ! change croot bio per diam [kgC/cm] + + real(r8) :: bag ! above ground biomass [kgC] + real(r8) :: dbagdd ! change in agb per diameter [kgC/cm] + + select case(int(EDPftvarcon_inst%allom_cmode(ipft))) + case(1) !"constant") + call bag_allom(d,h,ipft,bag,dbagdd) + call bcr_const(d,bag,dbagdd,ipft,bcr,dbcrdd) + case DEFAULT + write(fates_log(),*) 'An undefined coarse root allometry was specified: ', & + EDPftvarcon_inst%allom_cmode(ipft) + write(fates_log(),*) 'Aborting' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end select + return + end subroutine bcr_allom + + ! ============================================================================ + ! Fine root biomass allometry wrapper + ! ============================================================================ + + subroutine bfineroot(d,h,ipft,canopy_trim,bfr,dbfrdd) + + ! ------------------------------------------------------------------------- + ! This subroutine calculates the actual target fineroot biomass + ! based on functions that may or may not have prognostic properties. + ! ------------------------------------------------------------------------- + + real(r8),intent(in) :: d ! plant diameter [cm] + real(r8),intent(in) :: h ! plant height [m] + integer(i4),intent(in) :: ipft ! PFT index + real(r8),intent(in) :: canopy_trim ! trimming function + real(r8),intent(out) :: bfr ! fine root biomass [kgC] + real(r8),intent(out),optional :: dbfrdd ! change leaf bio per diameter [kgC/cm] + + real(r8) :: blmax ! maximum leaf biomss per allometry + real(r8) :: dblmaxdd + real(r8) :: bfrmax + real(r8) :: dbfrmaxdd + real(r8) :: slascaler + + select case(int(EDPftvarcon_inst%allom_fmode(ipft))) + case(1) ! "constant proportionality with bleaf" + + call blmax_allom(d,h,ipft,blmax,dblmaxdd) + call bfrmax_const(d,blmax,dblmaxdd,ipft,bfrmax,dbfrmaxdd) + bfr = bfrmax * canopy_trim + if(present(dbfrdd))then + dbfrdd = dbfrmaxdd * canopy_trim + end if + case DEFAULT + write(fates_log(),*) 'An undefined fine root allometry was specified: ', & + EDPftvarcon_inst%allom_fmode(ipft) + write(fates_log(),*) 'Aborting' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end select + + return + end subroutine bfineroot + + + ! ============================================================================ + ! Dead biomass interface + ! ============================================================================ + + subroutine bdead_allom(bag,bcr,bsap,ipft,bdead,dbagdd,dbcrdd,dbsapdd,dbdeaddd) + + + real(r8),intent(in) :: bag ! agb [kgC] + real(r8),intent(in) :: bcr ! coarse root biomass [kgC] + real(r8),intent(in) :: bsap ! sapwood biomass [kgC] + integer(i4),intent(in) :: ipft ! PFT index + real(r8),intent(out) :: bdead ! dead biomass (heartw/struct) [kgC] + + real(r8),intent(in),optional :: dbagdd ! change in agb per d [kgC/cm] + real(r8),intent(in),optional :: dbcrdd ! change in croot per d [kgC/cm] + real(r8),intent(in),optional :: dbsapdd ! change in bsap per d [kgC/cm] + real(r8),intent(out),optional :: dbdeaddd ! change in bdead per d [kgC/cm] + + ! bdead is diagnosed as the mass balance from all other pools + ! and therefore, no options are necessary + + associate( agb_fraction => EDPftvarcon_inst%allom_agb_frac(ipft)) + + select case(int(EDPftvarcon_inst%allom_amode(ipft))) + case(1) ! Saldariagga mass allometry originally calculated bdead directly. + ! we assume proportionality between bdead and bag + + bdead = bag/agb_fraction + if(present(dbagdd) .and. present(dbdeaddd))then + dbdeaddd = dbagdd/agb_fraction + end if + + case(2,3) + + bdead = bag+bcr-bsap + if(present(dbagdd) .and. present(dbcrdd) .and. & + present(dbdeaddd) .and. present(dbsapdd) )then + dbdeaddd = dbagdd+dbcrdd-dbsapdd + end if + + case DEFAULT + + write(fates_log(),*) 'An undefined AGB allometry was specified: ',& + EDPftvarcon_inst%allom_amode(ipft) + write(fates_log(),*) 'Aborting' + call endrun(msg=errMsg(sourcefile, __LINE__)) + + end select + end associate + return + end subroutine bdead_allom + + ! ============================================================================ + ! Specific bfrmax relationships + ! ============================================================================ + + subroutine bfrmax_const(d,blmax,dblmaxdd,ipft,bfrmax,dbfrmaxdd) + + + real(r8),intent(in) :: d ! plant diameter [cm] + real(r8),intent(in) :: blmax ! max leaf biomass [kgC] + real(r8),intent(in) :: dblmaxdd ! change in blmax per diam [kgC/cm] + integer(i4),intent(in) :: ipft ! PFT index + real(r8),intent(out) :: bfrmax ! max fine-root root biomass [kgC] + real(r8),intent(out),optional :: dbfrmaxdd ! change frmax bio per diam [kgC/cm] + + associate( l2fr => EDPftvarcon_inst%allom_l2fr(ipft) ) + + bfrmax = blmax*l2fr + + ! dbfr/dd = dbfrmax/dblmax * dblmax/dd + if(present(dbfrmaxdd))then + dbfrmaxdd = dblmaxdd*l2fr + end if + + end associate + return + end subroutine bfrmax_const + + ! ============================================================================ + ! Specific bcr relationships + ! ============================================================================ + + subroutine bcr_const(d,bag,dbagdd,ipft,bcr,dbcrdd) + + real(r8),intent(in) :: d ! plant diameter [cm] + real(r8),intent(in) :: bag ! above ground biomass [kg] + real(r8),intent(in) :: dbagdd ! change in agb per diameter [kg/cm] + integer(i4),intent(in) :: ipft ! PFT index + real(r8),intent(out) :: bcr ! coarse root biomass [kg] + real(r8),intent(out),optional :: dbcrdd ! change croot bio per diam [kg/cm] + + associate( agb_fraction => EDPftvarcon_inst%allom_agb_frac(ipft) ) + + ! btot = bag + bcr + ! bag = btot*agb_fraction + ! bag/agb_fraction = bag + bcr + ! bcr = bag*(1/agb_fraction-1) + bcr = bag*(1.0_r8/agb_fraction-1.0_r8) + + ! Derivative + ! dbcr/dd = dbcr/dbag * dbag/dd + if(present(dbcrdd))then + dbcrdd = (1.0_r8/agb_fraction-1.0_r8)*dbagdd + end if + + end associate + return + end subroutine bcr_const + + ! ============================================================================ + ! Specific d2bsap relationships + ! ============================================================================ + + subroutine bsap_dlinear(d,h,dhdd,blmax,dblmaxdd,bag,dbagdd,ipft,bsap,dbsapdd) + + use FatesConstantsMod, only : g_per_kg + use FatesConstantsMod, only : cm2_per_m2 + use FatesConstantsMod, only : kg_per_Megag + + ! ------------------------------------------------------------------------- + ! Calculate sapwood biomass based on leaf area to sapwood area + ! proportionality. In this function, the leaftosapwood area is a function + ! of plant size, see Calvo-Alvarado and Bradley Christoferson + ! In this case: parameter latosa (from constant proportionality) + ! is the intercept of the diameter function. + ! + ! For very small plants, the fraction can get very large, so cap the amount + ! of sapwood at X! of agb-bleaf + ! ------------------------------------------------------------------------- + + real(r8),intent(in) :: d ! plant diameter [cm] + real(r8),intent(in) :: h ! plant height [m] + real(r8),intent(in) :: dhdd ! change in height per diameter [m/cm] + real(r8),intent(in) :: blmax ! plant leaf biomass [kgC] + real(r8),intent(in) :: dblmaxdd ! change in blmax per diam [kgC/cm] + real(r8),intent(in) :: bag ! aboveground biomass [kgC] + real(r8),intent(in) :: dbagdd ! change in agb per diam [kgC/cm] + integer(i4),intent(in) :: ipft ! PFT index + real(r8),intent(out) :: bsap ! plant leaf biomass [kgC] + real(r8),intent(out),optional :: dbsapdd ! change leaf bio per diameter [kgC/cm] + + + real(r8) :: latosa ! applied leaf area to sap area + ! may or may not contain diameter correction + real(r8) :: hbl2bsap ! sapwood biomass per lineal height and kg of leaf + + ! Constrain sapwood to be no larger than 75% of total agb + real(r8),parameter :: max_agbfrac = 0.75_r8 + + associate ( latosa_int => EDPftvarcon_inst%allom_latosa_int(ipft), & + latosa_slp => EDPftvarcon_inst%allom_latosa_slp(ipft), & + sla => EDPftvarcon_inst%slatop(ipft), & + wood_density => EDPftvarcon_inst%wood_density(ipft), & + c2b => EDPftvarcon_inst%c2b(ipft)) + + ! ------------------------------------------------------------------------ + ! Calculate sapwood biomass per linear height and kgC of leaf [m-1] + ! Units: + ! (1/latosa)* slatop* gtokg * cm2tom2 / c2b * mg2kg * dens + ! [cm2/m2]*[m2/gC]*[1000gC/1kgC]*[1m2/10000cm2] /[kg/kgC]*[kg/Mg]*[Mg/m3] + ! ->[cm2/gC] + ! ->[cm2/kgC] + ! ->[m2/kgC] + ! ->[m2/kg] + ! ->[m2/Mg] + ! ->[/m] + ! ------------------------------------------------------------------------ + + if(test_b4b) then + + bsap = blmax * latosa_int * h + + if(present(dbsapdd))then + dbsapdd = latosa_int*(h*dblmaxdd + blmax*dhdd) + end if + + else + + latosa = latosa_int + d*latosa_slp + hbl2bsap = sla*g_per_kg*wood_density*kg_per_Megag/(latosa*c2b*cm2_per_m2 ) + + ! Force sapwood to be less than a maximum fraction of total alive biomass + ! (this comes into play typically in very small plants) + bsap = min(max_agbfrac*bag,hbl2bsap * h * blmax) + + ! Derivative + ! dbldmaxdd is deriv of blmax wrt dbh (use directives to check oop) + ! dhdd is deriv of height wrt dbh (use directives to check oop) + if(present(dbsapdd))then + if (bsap=dbh_maxh)then + dblmaxdd = 0.0_r8 + else + dblmaxdd = p1*p2*d**(p2-1.0_r8) / c2b + end if + end if + + return + end subroutine dh2blmax_2pwr + + ! ========================================================================= + ! Diameter to height (D2H) functions + ! ========================================================================= + + subroutine d2h_chave2014(d,p1,p2,p3,dbh_maxh,h,dhdd) + + ! "d2h_chave2014" + ! "d to height via Chave et al. 2014" + + ! This function calculates tree height based on tree diameter and the + ! environmental stress factor "E", as per Chave et al. 2015 GCB + ! As opposed to previous allometric models in ED, in this formulation + ! we do not impose a hard cap on tree height. But, maximum_height + ! is an important parameter, but instead of imposing a hard limit, in + ! the new methodology, it will be used to trigger a change in carbon + ! balance accounting. Such that a tree that hits its maximum height will + ! begin to route available NPP into seed and defense respiration. + ! + ! The stress function is based on the geographic location of the site. If + ! a user decides to use Chave2015 allometry, the E factor will be read in + ! from a global gridded dataset and assigned for each ED patch (note it + ! will be the same for each ED patch, but this distinction will help in + ! porting ED into different models (patches are pure ED). It + ! assumes that the site is within the pan-tropics, and is a linear function + ! of climatic water deficit, temperature seasonality and precipitation + ! seasonality. See equation 6b of Chave et al. + ! The relevant equation for height in this function is 6a of the same + ! manuscript, and is intended to pair with diameter to relate with + ! structural biomass as per equation 7 (in which H is implicit). + ! + ! Chave et al. Improved allometric models to estimate the abovegroud + ! biomass of tropical trees. Global Change Biology. V20, p3177-3190. 2015. + ! + ! ========================================================================= + + real(r8),intent(in) :: d ! plant diameter [cm] + real(r8),intent(in) :: p1 ! parameter a + real(r8),intent(in) :: p2 ! parameter b + real(r8),intent(in) :: p3 ! parameter c + real(r8),intent(in) :: dbh_maxh ! dbh at maximum height [cm] + + real(r8),intent(out) :: h ! plant height [m] + real(r8),intent(out),optional :: dhdd ! change in height per diameter [m/cm] + real(r8) :: p1e + + p1e = p1 ! -eclim (assumed that p1 already has eclim removed) + if(d>=dbh_maxh) then + h = exp( p1e + p2*log(dbh_maxh) + p3*log(dbh_maxh)**2.0 ) + else + h = exp( p1e + p2*log(d) + p3*log(d)**2.0 ) + end if + + if(present(dhdd))then + if(d>=dbh_maxh ) then + dhdd = 0.0_r8 + else + dhdd = exp(p1e) * ( 2.0_r8*p3*d**(p2-1.0_r8+p3*log(d))*log(d) + & + p2*d**(p2-1.0_r8+p3*log(d)) ) + end if + end if + return + end subroutine d2h_chave2014 + + ! =========================================================================== + + subroutine d2h_poorter2006(d,p1,p2,p3,dbh_maxh,h,dhdd) + + ! "d2h_poorter2006" + ! "d to height via Poorter et al. 2006, these routines use natively + ! asymtotic functions" + ! + ! Poorter et al calculated height diameter allometries over a variety of + ! species in Bolivia, including those that could be classified in guilds + ! that were Partial-shade-tolerants, long-lived pioneers, shade-tolerants + ! and short-lived pioneers. There results between height and diameter + ! found that small stature trees had less of a tendency to asymotote in + ! height and showed more linear relationships, and the largest stature + ! trees tended to show non-linear relationships that asymtote. + ! + ! h = h_max*(1-exp(-a*d**b)) + ! + ! Poorter L, L Bongers and F Bongers. Architecture of 54 moist-forest tree + ! species: traits, trade-offs, and functional groups. Ecology 87(5), 2006. + ! + ! ========================================================================= + + real(r8),intent(in) :: d ! plant diameter [cm] + real(r8),intent(in) :: p1 ! parameter a = h_max + real(r8),intent(in) :: p2 ! parameter b + real(r8),intent(in) :: p3 ! parameter c + real(r8),intent(in) :: dbh_maxh ! dbh at maximum height + real(r8),intent(out) :: h ! plant height [m] + real(r8),intent(out),optional :: dhdd ! change in height per diameter [m/cm] + + h = p1*(1.0_r8 - exp(p2*min(d,dbh_maxh)**p3)) + + !h = h_max - h_max (exp(a*d**b)) + !f(x) = -h_max*exp(g(x)) + !g(x) = a*d**b + !d/dx f(g(x) = f'(g(x))*g'(x) = -a1*exp(a2*d**a3) * a3*a2*d**(a3-1) + + if(present(dhdd))then + if( d>=dbh_maxh ) then + dhdd = 0.0_r8 + else + dhdd = -p1*exp(p2*d**p3) * p3*p2*d**(p3-1.0_r8) + end if + end if + + return + end subroutine d2h_poorter2006 + + ! =========================================================================== + + subroutine d2h_2pwr(d,p1,p2,dbh_maxh,h,dhdd) + + ! ========================================================================= + ! "d2h_2pwr" + ! "d to height via 2 parameter power function" + ! where height h is related to diameter by a linear relationship on the log + ! transform where log(a) is the intercept and b is the slope parameter. + ! + ! log(h) = log(a) + b*log(d) + ! h = exp(log(a)) * exp(log(d))**b + ! h = a*d**b + ! + ! This functional form is used often in temperate allometries + ! Therefore, no base reference is cited. Although, the reader is pointed + ! towards Dietze et al. 2008, King 1991, Ducey 2012 and many others for + ! reasonable parameters. Note that this subroutine is intended only for + ! trees well below their maximum height, ie initialization. + ! + ! ========================================================================= + ! From King et al. 1990 at BCI for saplings + ! log(d) = a + b*log(h) + ! d = exp(a) * h**b + ! h = (1/exp(a)) * d**(1/b) + ! h = p1*d**p2 where p1 = 1/exp(a) = 1.07293 p2 = 1/b = 1.4925 + ! d = (h/p1)**(1/p2) + ! For T. tuberculata (canopy tree) a = -0.0704, b = 0.67 + ! ========================================================================= + + ! args + ! ========================================================================= + ! d: diameter at breast height + ! p1: the intercept parameter + ! (however exponential of the fitted log trans) + ! p2: the slope parameter + ! return: + ! h: total tree height [m] + ! ========================================================================= + + real(r8),intent(in) :: d ! plant diameter [cm] + real(r8),intent(in) :: p1 ! parameter a + real(r8),intent(in) :: p2 ! parameter b + real(r8),intent(in) :: dbh_maxh ! dbh where max height occurs [cm] + real(r8),intent(out) :: h ! plant height [m] + real(r8),intent(out),optional :: dhdd ! change in height per diameter [m/cm] + + h = p1*min(d,dbh_maxh)**p2 + + if(present(dhdd))then + if( d>=dbh_maxh ) then + dhdd = 0.0_r8 + else + dhdd = (p2*p1)*d**(p2-1.0_r8) + end if + end if + + return + end subroutine d2h_2pwr + + ! ============================================================================ + + subroutine d2h_obrien(d,p1,p2,dbh_maxh,h,dhdd) + + real(r8),intent(in) :: d ! plant diameter [cm] + real(r8),intent(in) :: p1 ! parameter a + real(r8),intent(in) :: p2 ! parameter b + real(r8),intent(in) :: dbh_maxh ! dbh where max height occurs [cm] + real(r8),intent(out) :: h ! plant height [m] + real(r8),intent(out),optional :: dhdd ! change in height per diameter [m/cm] + + !p1 = 0.64 + !p2 = 0.37 + h = 10.0_r8**(log10(min(d,dbh_maxh))*p1+p2) + + if(present(dhdd))then + if(d>=dbh_maxh ) then + dhdd = 0.0_r8 + else + dhdd = p1*10.0_r8**p2*d**(p1-1.0_r8) + end if + end if + + + return + end subroutine d2h_obrien + + ! =========================================================================== + + subroutine d2h_martcano(d,p1,p2,p3,dbh_maxh,h,dhdd) + + ! ========================================================================= + ! "d2h_martcano" + ! "d to height via 3 parameter Michaelis-Menten following work at BCI + ! by Martinez-Cano et al. 2016 + ! + ! h = (a*d**b)/(c+d**b) + ! + ! h' = [(a*d**b)'(c+d**b) - (c+d**b)'(a*d**b)]/(c+d**b)**2 + ! dhdd = h' = [(ba*d**(b-1))(c+d**b) - (b*d**(b-1))(a*d**b)]/(c+d**b)**2 + ! + ! args + ! ========================================================================= + ! d: diameter at breast height + ! h: total tree height [m] + ! ========================================================================= + + real(r8),intent(in) :: d ! plant diameter [cm] + real(r8),intent(in) :: p1 ! parameter a + real(r8),intent(in) :: p2 ! parameter b + real(r8),intent(in) :: p3 ! parameter c + real(r8),intent(in) :: dbh_maxh ! diameter at maximum height + real(r8),intent(out) :: h ! plant height [m] + real(r8),intent(out),optional :: dhdd ! change in height per diameter [m/cm] + + h = (p1*min(d,dbh_maxh)**p2)/(p3+min(d,dbh_maxh)**p2) + + if(present(dhdd))then + if(d>=dbh_maxh ) then + dhdd = 0.0 + else + dhdd = ((p2*p1*d**(p2-1._r8))*(p3+d**p2) - & + (p2*d**(p2-1._r8))*(p1*d**p2) )/ & + (p3+d**p2)**2._r8 + end if + end if + + return + end subroutine d2h_martcano + + + ! ========================================================================= + ! Diameter to (2) above-ground biomass + ! ========================================================================= + + subroutine dh2bag_chave2014(d,h,dhdd,p1,p2,wood_density,c2b,bag,dbagdd) + + ! ========================================================================= + ! This function calculates tree structural biomass from tree diameter, + ! height and wood density. + ! + ! Chave et al. Improved allometric models to estimate the abovegroud + ! biomass of tropical trees. Global Change Biology. V20, p3177-3190. 2015. + ! + ! Input arguments: + ! d: Diameter at breast height [cm] + ! rho: wood specific gravity (dry wood mass per green volume) + ! height: total tree height [m] + ! a1: structural biomass allometry parameter 1 (intercept) + ! a2: structural biomass allometry parameter 2 (slope) + ! Output: + ! bag: Total above ground biomass [kgC] + ! + ! ========================================================================= + + + real(r8),intent(in) :: d ! plant diameter [cm] + real(r8),intent(in) :: h ! plant height [m] + real(r8),intent(in) :: dhdd ! change in height wrt diameter + real(r8),intent(in) :: p1 ! allometry parameter 1 + real(r8),intent(in) :: p2 ! allometry parameter 2 + real(r8),intent(in) :: wood_density + real(r8),intent(in) :: c2b + real(r8),intent(out) :: bag ! plant height [m] + real(r8),intent(out),optional :: dbagdd ! change in agb per diameter [kgC/cm] + + real(r8) :: dbagdd1,dbagdd2,dbagdd3 + + bag = (p1 * (wood_density*d**2.0_r8*h)**p2)/c2b + + + if(present(dbagdd))then + ! Need the the derivative of height to diameter to + ! solve the derivative of agb with height + dbagdd1 = (p1*wood_density**p2)/c2b + dbagdd2 = p2*d**(2.0_r8*p2)*h**(p2-1.0_r8)*dhdd + dbagdd3 = h**p2*2.0_r8*p2*d**(2.0_r8*p2-1.0_r8) + dbagdd = dbagdd1*(dbagdd2 + dbagdd3) + end if + + return + end subroutine dh2bag_chave2014 + + subroutine d2bag_2pwr(d,p1,p2,c2b,bag,dbagdd) + + ! ========================================================================= + ! This function calculates tree above ground biomass according to 2 + ! parameter power functions. (slope and intercepts of a log-log + ! diameter-agb fit: + ! + ! These relationships are typical for temperate/boreal plants in North + ! America. Parameters are available from Chojnacky 2014 and Jenkins 2003 + ! + ! Note that we are using an effective diameter here, as Chojnacky 2014 + ! and Jenkins use diameter at the base of the plant for "woodland" species + ! The diameters should be converted prior to this routine if drc. + ! + ! Input arguments: + ! diam: effective diameter (d or drc) in cm + ! FOR SPECIES THAT EXPECT DCM, THIS NEEDS TO BE PRE-CALCULATED!!!! + ! Output: + ! agb: Total above ground biomass [kgC] + ! + ! ========================================================================= + ! Aceraceae, Betulaceae, Fagaceae and Salicaceae comprised nearly + ! three-fourths of the hardwood species (Table 3) + ! + ! Fabaceae and Juglandaceae had specific gravities .0.60 and were + ! combined, as were Hippocastanaceae and Tilaceae with specific gravities + ! near 0.30. The remaining 9 families, which included mostly species with + ! specific gravity 0.45–0.55, were initially grouped to construct a general + ! hardwood taxon for those families having few published biomass equa- + ! tions however, 3 warranted separation, leaving 6 families for the general + ! taxon. + ! bag = exp(b0 + b1*ln(diameter))/c2b + ! ========================================================================= + + + real(r8),intent(in) :: d ! plant diameter [cm] + real(r8),intent(in) :: p1 ! allometry parameter 1 + real(r8),intent(in) :: p2 ! allometry parameter 2 + real(r8),intent(in) :: c2b ! carbon to biomass multiplier ~2 + real(r8),intent(out) :: bag ! plant height [m] + real(r8),intent(out),optional :: dbagdd ! change in agb per diameter [kgC/cm] + + bag = (p1 * d**p2)/c2b + if(present(dbagdd))then + dbagdd = (p2*p1*d**(p2-1.0_r8))/c2b + end if + + return + end subroutine d2bag_2pwr + + + subroutine dh2bag_salda(d,h,dhdd,p1,p2,p3,p4, & + wood_density,c2b,allom_agb_frac,bag,dbagdd) + + ! -------------------------------------------------------------------- + ! Calculate stem biomass from height(m) dbh(cm) and wood density(g/cm3) + ! default params using allometry of J.G. Saldarriaga et al 1988 - Rio Negro + ! Journal of Ecology vol 76 p938-958 + ! Saldarriaga 1988 provided calculations on total dead biomass + ! So here, we calculate total dead, and then remove the below-ground + ! fraction + ! -------------------------------------------------------------------- + + real(r8),intent(in) :: d ! plant diameter [cm] + real(r8),intent(in) :: h ! plant height [m] + real(r8),intent(in) :: dhdd ! change in height wrt diameter + real(r8),intent(in) :: p1 ! = 0.06896_r8 + real(r8),intent(in) :: p2 ! = 0.572_r8 + real(r8),intent(in) :: p3 ! = 1.94_r8 + real(r8),intent(in) :: p4 ! = 0.931_r8 + real(r8),intent(in) :: c2b ! carbon 2 biomass ratio + real(r8),intent(in) :: wood_density + real(r8),intent(in) :: allom_agb_frac + real(r8),intent(out) :: bag ! plant biomass [kgC/indiv] + real(r8),intent(out),optional :: dbagdd ! change in agb per diameter [kgC/cm] + + real(r8) :: term1,term2,term3 + + + bag = allom_agb_frac*p1*(h**p2)*(d**p3)*(wood_density**p4) + + ! Add sapwood calculation to this? + + ! bag = a1 * h**a2 * d**a3 * r**a4 + ! dbag/dd = a1*r**a4 * d/dd (h**a2 * d**a3) + ! dbag/dd = a1*r**a4 * [ d**a3 *d/dd(h**a2) + h**a2*d/dd(d**a3) ] + ! dbag/dd = a1*r**a4 * [ d**a3 * a2*h**(a2-1)dh/dd + h**a2*a3*d**(a3-1)] + + ! From code + ! dbag/dd = a3 * a1 *(h**a2)*(d**(a3-1))* (r**a4) + a2*a1*(h**(a2-1))*(d**a3)*(r**a4)*dhdd + ! dbag/dd = a1*r**a4 * [ d**a3 * a2* h**(a2-1)*dhdd + a3 * (h**a2)*(d**(a3-1)) ] + + + if(present(dbagdd)) then + + term1 = allom_agb_frac*p1*(wood_density**p4) + term2 = (h**p2)*p3*d**(p3-1.0_r8) + term3 = p2*h**(p2-1.0_r8)*(d**p3)*dhdd + dbagdd = term1*(term2+term3) + + end if + + return + end subroutine dh2bag_salda + + ! ============================================================================ + ! height to diameter conversions + ! Note that these conversions will only back-calculate the actual diameter + ! for plants that have not started to experience height capping or an + ! asymptote. In these cases they may be called effective diameter. + ! ============================================================================ + + subroutine h2d_chave2014(h,p1,p2,p3,de,ddedh) + + + real(r8),intent(in) :: h ! plant height [m] + real(r8),intent(in) :: p1 + real(r8),intent(in) :: p2 + real(r8),intent(in) :: p3 + + real(r8),intent(out) :: de ! effective plant diameter [cm] + real(r8),intent(out),optional :: ddedh ! effective change in d per height [cm/m] + + real(r8) :: p1e, eroot, dbh1,dhpdd + + p1e = p1 !-eclim (assumed that p1 already has eclim removed) + eroot = (-p2 + sqrt(p2**2.0_r8 + 4.0_r8*log(h*exp(-p1e))*p3)) & + /(2.0_r8*p3) + + de = exp(eroot) + + if(present(ddedh))then + ! Invert the derivative at d without asymtote + dhpdd = exp(p1e)*( p3*2.0_r8*de**(p2-1.0_r8)*log(de)* & + exp(p3*log(de)**2) + p2*de**(p2-1.0_r8)* & + exp(p3*log(de)**2.0_r8) ) + + ddedh = 1.0_r8/dhpdd + end if + + ! term1 = exp(-p2/(2*p3)) + ! term2 = exp(p2**2/(4*p3**2)) + ! term3 = exp(-p1e/p3) + ! term4 = h**(1/p3-1.0_r8)/(p3) + ! d = term1*term2*term3*term4 + return + end subroutine h2d_chave2014 + + ! ============================================================================ + + subroutine h2d_poorter2006(h,p1,p2,p3,d,dddh) + + ! ------------------------------------------------------------------------- + ! Note that the height to diameter conversion in poorter is only necessary + ! when initializing saplings. In other methods, height to diameter is + ! useful when defining the dbh at which point to asymptote, as maximum + ! height is the user set parameter. This function should not need to set a + ! dbh_max parameter for instance, but it may end up doing so anyway, even + ! if it is not used, do to poor filtering. The poorter et al. d2h and h2d + ! functions are already asymptotic, and the parameter governing maximum + ! height is the p1 parameter. Note as dbh gets very large, the + ! exponential goes to zero, and the maximum height approaches p1. + ! However, the asymptote has a much different shape than the logistic, so + ! therefore in the Poorter et al functions, we do not set p1 == h_max. + ! That being said, if an h_max that is greater than p1 is passed to this + ! function, it will return a complex number. During parameter + ! initialization, a check will be placed that forces: + ! h_max = p1*0.98 + ! ------------------------------------------------------------------------- + + real(r8),intent(in) :: h ! plant height [m] + real(r8),intent(in) :: p1 + real(r8),intent(in) :: p2 + real(r8),intent(in) :: p3 + + real(r8),intent(out) :: d ! plant diameter [cm] + real(r8),intent(out),optional :: dddh ! change in d per height [cm/m] + + ! ------------------------------------------------------------------------- + ! h = a1*(1 - exp(a2*d**a3)) + ! h = a1 - a1*exp(a2*d**a3) + ! a1-h = a1*exp(a2*d**a3) + ! (a1-h)/a1 = exp(a2*d**a3) + ! log(1-h/a1) = a2*d**a3 + ! [log(1-h/a1)/a2]**(1/a3) = d + ! + ! derivative dd/dh + ! dd/dh = [log((a1-h)/a1)/a2]**(1/a3)' + ! = (1/a3)*[log((a1-h)/a1)/a2]**(1/a3-1)* [(log(a1-h)-log(a1))/a2]' + ! = (1/a3)*[log((a1-h)/a1)/a2]**(1/a3-1) * (1/(a2*(h-a1)) + ! dd/dh = -((log(1-h/a1)/a2)**(1/a3-1))/(a2*a3*(a1-h)) + ! ------------------------------------------------------------------------- + + d = (log(1.0_r8-h/p1)/p2)**(1.0_r8/p3) + + if(present(dddh))then + dddh = -((log(1-h/p1)/p2)**(1.0_r8/p3-1.0_r8))/ & + (p2*p3*(p1-h)) + end if + + return + end subroutine h2d_poorter2006 + + ! ============================================================================ + + subroutine h2d_2pwr(h,p1,p2,d,dddh) + + + real(r8),intent(in) :: h ! plant height [m] + real(r8),intent(in) :: p1 ! parameter 1 + real(r8),intent(in) :: p2 ! parameter 2 + + real(r8),intent(out) :: d ! plant diameter [cm] + real(r8),intent(out),optional :: dddh ! change in d per height [cm/m] + + !h = a1*d**a2 + d = (h/p1)**(1.0_r8/p2) + + ! d = (1/a1)**(1/a2)*h**(1/a2) + if(present(dddh)) then + dddh = (1.0_r8/p2)*(1.0_r8/p1)**(1.0_r8/p2) & + *h**(1.0_r8/p2-1.0_r8) + end if + + return + end subroutine h2d_2pwr + + ! ============================================================================ + + subroutine h2d_obrien(h,p1,p2,d,dddh) + + real(r8),intent(in) :: h ! plant height [m] + real(r8),intent(in) :: p1 + real(r8),intent(in) :: p2 + + real(r8),intent(out) :: d ! plant diameter [cm] + real(r8),intent(out),optional :: dddh ! change in d per height [cm/m] + + d = 10.0_r8**((log10(h)-p2)/p1) + + if(present(dddh))then + dddh = 1.0_r8/(p1*10.0_r8**p2*d**(p1-1.0_r8)) + end if + + return + end subroutine h2d_obrien + + ! ============================================================================ + + subroutine h2d_martcano(h,p1,p2,p3,d,dddh) + + ! ========================================================================= + ! "d2h_martcano" + ! "d to height via 3 parameter Michaelis-Menten following work at BCI + ! by Martinez-Cano et al. 2016 + ! + ! h = (a*d**b)/(c+d**b) + ! + ! d = [(h*c)/(a-h)]**(1/b) + ! d = [(h*c)**(1/b)] / [(a-h)**(1/b)] + ! d' = {[(h*c)**(1/b)]' [(a-h)**(1/b)] - [(a-h)**(1/b)]'[(h*c)**(1/b)]} / + ! [(a-h)**(1/b)]**2 + ! dddh = d' = {[(1/b)(h*c)**(1/b-1)] [(a-h)**(1/b)] - + ! [(1/b)(a-h)**(1/b-1)] [(h*c)**(1/b)]} / + ! [(a-h)**(1/b)]**2 + ! + ! ========================================================================= + + real(r8),intent(in) :: h ! plant height [m] + real(r8),intent(in) :: p1 + real(r8),intent(in) :: p2 + real(r8),intent(in) :: p3 + + real(r8),intent(out) :: d ! plant diameter [cm] + real(r8),intent(out),optional :: dddh ! change in diameter per height [cm/m] + + d = ((h*p3)/(p1-h))**(1._r8/p2) + + if(present(dddh))then + dddh = (((1._r8/p2)*(h*p3)**(1._r8/p2-1._r8))*((p1-h)**(1._r8/p2)) - & + ((1._r8/p2)*(p1-h)**(1._r8/p2-1._r8))* ((h*p3)**(1._r8/p2)) ) / & + ((p1-h)**(1._r8/p2))**2._r8 + end if + return + end subroutine h2d_martcano + + ! ============================================================================= + ! Specific diameter to crown area allometries + ! ============================================================================= + + + subroutine carea_2pwr(d,spread,d2bl_p2,d2bl_ediff,d2ca_min,d2ca_max,c_area) + + ! ============================================================================ + ! Calculate area of ground covered by entire cohort. (m2) + ! Function of DBH (cm) canopy spread (m/cm) and number of individuals. + ! ============================================================================ + + real(r8),intent(in) :: d ! diameter [cm] + real(r8),intent(in) :: spread ! site level relative spread score [0-1] + real(r8),intent(in) :: d2bl_p2 ! parameter 2 in the diameter->bleaf allometry (exponent) + real(r8),intent(in) :: d2bl_ediff ! area difference factor in the diameter-bleaf allometry (exponent) + real(r8),intent(in) :: d2ca_min ! minimum diameter to crown area scaling factor + real(r8),intent(in) :: d2ca_max ! maximum diameter to crown area scaling factor + real(r8),intent(out) :: c_area ! crown area for one plant [m2] + + real(r8) :: crown_area_to_dbh_exponent + real(r8) :: spreadterm ! Effective 2bh to crown area scaling factor + + ! default is to use the same exponent as the dbh to bleaf exponent so that per-plant + ! canopy depth remains invariant during growth, but allowed to vary via the + ! allom_blca_expnt_diff term (which has default value of zero) + crown_area_to_dbh_exponent = d2bl_p2 + d2bl_ediff + + ! ---------------------------------------------------------------------------------- + ! The function c_area is called during the process of canopy position demotion + ! and promotion. As such, some cohorts are temporarily elevated to canopy positions + ! that are outside the number of alloted canopy spaces. Ie, a two story canopy + ! may have a third-story plant, if only for a moment. However, these plants + ! still need to generate a crown area to complete the promotion, demotion process. + ! So we allow layer index exceedence here and force it down to max. + ! (rgk/cdk 05/2017) + ! ---------------------------------------------------------------------------------- + + ! apply site-level spread elasticity to the cohort crown allometry term + + spreadterm = spread * d2ca_max + (1._r8 - spread) * d2ca_min + + c_area = spreadterm * d ** crown_area_to_dbh_exponent + + end subroutine carea_2pwr + + + ! =========================================================================== + + subroutine cspline(x1,x2,y1,y2,dydx1,dydx2,x,y,dydx) + + ! ============================================================================ + ! This subroutine performs a cubic spline interpolation between known + ! endpoints. The endpoints have known coordinats and slopes + ! + ! This routine may come in handy if we ever want to start using allometries + ! that are different for juvenile and adult plants, and then connect + ! the curves smoothly + ! + ! ============================================================================ + + ! Arguments + + real(r8),intent(in) :: x1 ! Lower endpoint independent + real(r8),intent(in) :: x2 ! Upper endpoint independent + real(r8),intent(in) :: y1 ! Lower endpoint dependent + real(r8),intent(in) :: y2 ! Upper endpoint dependent + real(r8),intent(in) :: dydx1 ! Lower endpoint slope + real(r8),intent(in) :: dydx2 ! Upper endpoint slope + real(r8),intent(in) :: x ! Independent + real(r8),intent(out) :: y ! Dependent + real(r8),intent(out) :: dydx ! Slope + + ! Temps + real(r8) :: t + real(r8) :: a + real(r8) :: b + + t = (x-x1)/(x2-x1) + a = dydx1*(x2-x1) - (y2-y1) + b = -dydx2*(x2-x1) + (y2-y1) + + y = (1.0_r8-t)*y1 + t*y2 + t*(1.0_r8-t)*(a*(1.0_r8-t) + b*t) + dydx = (y2-y1)/(x2-x1) + (1.0_r8-2.0_r8*t)*(a*(1.0_r8-t)+b*t)/(x2-x1) + t*(1.0_r8-t)*(b-a)/(x2-x1) + return + end subroutine cspline + +end module FatesAllometryMod diff --git a/biogeophys/FatesPlantRespPhotosynthMod.F90 b/biogeophys/FatesPlantRespPhotosynthMod.F90 index 0cfd778323..2154e20352 100644 --- a/biogeophys/FatesPlantRespPhotosynthMod.F90 +++ b/biogeophys/FatesPlantRespPhotosynthMod.F90 @@ -193,6 +193,10 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) ! ----------------------------------------------------------------------------------- real(r8), dimension(2) :: bbbopt + logical, parameter :: test_b4b = .true. ! Leave in place some questionable allometry + ! while modular allometry is being introduce + ! to preserve results with previous release + associate( & c3psn => EDPftvarcon_inst%c3psn , & @@ -493,15 +497,16 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) ! fine root pools. ! ------------------------------------------------------------------ - leaf_frac = 1.0_r8/(currentCohort%canopy_trim + & - EDPftvarcon_inst%allom_latosa_int(currentCohort%pft) * & - currentCohort%hite + EDPftvarcon_inst%allom_l2fr(currentCohort%pft)) - - - currentCohort%bsw = EDPftvarcon_inst%allom_latosa_int(currentCohort%pft) * & - currentCohort%hite * & - (currentCohort%balive + currentCohort%laimemory)*leaf_frac - + if(test_b4b)then + leaf_frac = 1.0_r8/(currentCohort%canopy_trim + & + EDPftvarcon_inst%allom_latosa_int(currentCohort%pft) * & + currentCohort%hite + EDPftvarcon_inst%allom_l2fr(currentCohort%pft)) + + + currentCohort%bsw = EDPftvarcon_inst%allom_latosa_int(currentCohort%pft) * & + currentCohort%hite * & + (currentCohort%balive + currentCohort%laimemory)*leaf_frac + end if ! Calculate the amount of nitrogen in the above and below ground ! stem and root pools, used for maint resp diff --git a/fire/SFMainMod.F90 b/fire/SFMainMod.F90 old mode 100755 new mode 100644 diff --git a/main/EDInitMod.F90 b/main/EDInitMod.F90 old mode 100755 new mode 100644 index a7a7a34e5f..a6eed3097c --- a/main/EDInitMod.F90 +++ b/main/EDInitMod.F90 @@ -12,7 +12,6 @@ module EDInitMod use FatesGlobals , only : fates_log use FatesInterfaceMod , only : hlm_is_restart use EDPftvarcon , only : EDPftvarcon_inst - use EDGrowthFunctionsMod , only : bdead, bleaf, dbh use EDCohortDynamicsMod , only : create_cohort, fuse_cohorts, sort_cohorts use EDPatchDynamicsMod , only : create_patch use EDTypesMod , only : ed_site_type, ed_patch_type, ed_cohort_type @@ -26,6 +25,13 @@ module EDInitMod use FatesInterfaceMod , only : numpft use ChecksBalancesMod , only : SiteCarbonStock use FatesInterfaceMod , only : nlevsclass + use FatesAllometryMod , only : h2d_allom + use FatesAllometryMod , only : bag_allom + use FatesAllometryMod , only : bcr_allom + use FatesAllometryMod , only : bleaf + use FatesAllometryMod , only : bfineroot + use FatesAllometryMod , only : bsap_allom + use FatesAllometryMod , only : bdead_allom ! CIME GLOBALS use shr_log_mod , only : errMsg => shr_log_errMsg @@ -336,8 +342,13 @@ subroutine init_cohorts( patch_in, bc_in) ! ! !LOCAL VARIABLES: type(ed_cohort_type),pointer :: temp_cohort - integer :: cstatus - integer :: pft + integer :: cstatus + integer :: pft + real(r8) :: b_ag ! biomass above ground [kgC] + real(r8) :: b_cr ! biomass in coarse roots [kgC] + real(r8) :: b_leaf ! biomass in leaves [kgC] + real(r8) :: b_fineroot ! biomass in fine roots [kgC] + real(r8) :: b_sapwood ! biomass in sapwood [kgC] !---------------------------------------------------------------------- patch_in%tallest => null() @@ -352,27 +363,47 @@ subroutine init_cohorts( patch_in, bc_in) temp_cohort%pft = pft temp_cohort%n = EDPftvarcon_inst%initd(pft) * patch_in%area temp_cohort%hite = EDPftvarcon_inst%hgt_min(pft) - !temp_cohort%n = 0.5_r8 * 0.0028_r8 * patch_in%area ! BOC for fixed size runs EDPftvarcon_inst%initd(pft) * patch_in%area - !temp_cohort%hite = 28.65_r8 ! BOC translates to DBH of 50cm. EDPftvarcon_inst%hgt_min(pft) - temp_cohort%dbh = Dbh(temp_cohort) ! FIX(RF, 090314) - comment out addition of ' + 0.0001_r8*pft ' - seperate out PFTs a little bit... + + ! Calculate the plant diameter from height + call h2d_allom(temp_cohort%hite,pft,temp_cohort%dbh) + temp_cohort%canopy_trim = 1.0_r8 - temp_cohort%bdead = Bdead(temp_cohort) - temp_cohort%balive = Bleaf(temp_cohort)*(1.0_r8 + EDPftvarcon_inst%allom_l2fr(pft) & - + EDPftvarcon_inst%allom_latosa_int(temp_cohort%pft)*temp_cohort%hite) - temp_cohort%b = temp_cohort%balive + temp_cohort%bdead + + ! Calculate total above-ground biomass from allometry + call bag_allom(temp_cohort%dbh,temp_cohort%hite,pft,b_ag) + + ! Calculate coarse root biomass from allometry + call bcr_allom(temp_cohort%dbh,temp_cohort%hite,pft,b_cr) + + ! Calculate the leaf biomass + ! (calculates a maximum first, then applies canopy trim) + call bleaf(temp_cohort%dbh,temp_cohort%hite,pft,temp_cohort%canopy_trim,b_leaf) + + ! Calculate fine root biomass + ! (calculates a maximum and then trimming value) + call bfineroot(temp_cohort%dbh,temp_cohort%hite,pft,temp_cohort%canopy_trim,b_fineroot) + + ! Calculate sapwood biomass + call bsap_allom(temp_cohort%dbh,pft,temp_cohort%canopy_trim,b_sapwood) + + temp_cohort%balive = b_leaf + b_fineroot + b_sapwood + + call bdead_allom( b_ag, b_cr, b_sapwood, pft, temp_cohort%bdead ) + + temp_cohort%b = temp_cohort%balive + temp_cohort%bdead if( EDPftvarcon_inst%evergreen(pft) == 1) then - temp_cohort%bstore = Bleaf(temp_cohort) * EDPftvarcon_inst%cushion(pft) + temp_cohort%bstore = b_leaf * EDPftvarcon_inst%cushion(pft) temp_cohort%laimemory = 0._r8 cstatus = 2 endif if( EDPftvarcon_inst%season_decid(pft) == 1 ) then !for dorment places - temp_cohort%bstore = Bleaf(temp_cohort) * EDPftvarcon_inst%cushion(pft) !stored carbon in new seedlings. + temp_cohort%bstore = b_leaf * EDPftvarcon_inst%cushion(pft) !stored carbon in new seedlings. if(patch_in%siteptr%status == 2)then temp_cohort%laimemory = 0.0_r8 else - temp_cohort%laimemory = Bleaf(temp_cohort) + temp_cohort%laimemory = b_leaf endif ! reduce biomass according to size of store, this will be recovered when elaves com on. temp_cohort%balive = temp_cohort%balive - temp_cohort%laimemory @@ -380,8 +411,8 @@ subroutine init_cohorts( patch_in, bc_in) endif if ( EDPftvarcon_inst%stress_decid(pft) == 1 ) then - temp_cohort%bstore = Bleaf(temp_cohort) * EDPftvarcon_inst%cushion(pft) - temp_cohort%laimemory = Bleaf(temp_cohort) + temp_cohort%bstore = b_leaf * EDPftvarcon_inst%cushion(pft) + temp_cohort%laimemory = b_leaf temp_cohort%balive = temp_cohort%balive - temp_cohort%laimemory cstatus = patch_in%siteptr%dstatus endif diff --git a/main/EDMainMod.F90 b/main/EDMainMod.F90 old mode 100755 new mode 100644 index dc3a7fd9f5..b7d488ff32 --- a/main/EDMainMod.F90 +++ b/main/EDMainMod.F90 @@ -271,6 +271,7 @@ subroutine ed_integrate_state_variables(currentSite, bc_in ) currentCohort%dbh = max(small_no,currentCohort%dbh + currentCohort%ddbhdt * hlm_freq_day ) currentCohort%balive = currentCohort%balive + currentCohort%dbalivedt * hlm_freq_day currentCohort%bdead = max(small_no,currentCohort%bdead + currentCohort%dbdeaddt * hlm_freq_day ) + if ( DEBUG ) then write(fates_log(),*) 'EDMainMod dbstoredt I ',currentCohort%bstore, & currentCohort%dbstoredt,hlm_freq_day @@ -298,6 +299,7 @@ subroutine ed_integrate_state_variables(currentSite, bc_in ) currentCohort%gpp_acc = 0.0_r8 currentCohort%resp_acc = 0.0_r8 + call allocate_live_biomass(currentCohort,1) ! BOC...update tree 'hydraulic geometry' diff --git a/main/EDPftvarcon.F90 b/main/EDPftvarcon.F90 index 465a28fc59..c8f09739d8 100644 --- a/main/EDPftvarcon.F90 +++ b/main/EDPftvarcon.F90 @@ -69,8 +69,6 @@ module EDPftvarcon real(r8), allocatable :: smpso(:) real(r8), allocatable :: smpsc(:) real(r8), allocatable :: grperc(:) - - real(r8), allocatable :: bmort(:) real(r8), allocatable :: hf_sm_threshold(:) real(r8), allocatable :: vcmaxha(:) @@ -486,7 +484,7 @@ subroutine Register_PFT(this, fates_params) name = 'fates_allom_d2h3' call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) - + name = 'fates_allom_d2bl1' call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) @@ -1408,7 +1406,6 @@ subroutine FatesReportPFTParams(is_master) write(fates_log(),fmt0) 'root_long = ',EDPftvarcon_inst%root_long write(fates_log(),fmt0) 'clone_alloc = ',EDPftvarcon_inst%clone_alloc write(fates_log(),fmt0) 'seed_alloc = ',EDPftvarcon_inst%seed_alloc - write(fates_log(),fmt0) 'C2B = ',EDPftvarcon_inst%c2b write(fates_log(),fmt0) 'woody = ',EDPftvarcon_inst%woody write(fates_log(),fmt0) 'stress_decid = ',EDPftvarcon_inst%stress_decid write(fates_log(),fmt0) 'season_decid = ',EDPftvarcon_inst%season_decid @@ -1431,6 +1428,7 @@ subroutine FatesReportPFTParams(is_master) write(fates_log(),fmt0) 'smpso = ',EDPftvarcon_inst%smpso write(fates_log(),fmt0) 'smpsc = ',EDPftvarcon_inst%smpsc write(fates_log(),fmt0) 'grperc = ',EDPftvarcon_inst%grperc + write(fates_log(),fmt0) 'c2b = ',EDPftvarcon_inst%c2b write(fates_log(),fmt0) 'bmort = ',EDPftvarcon_inst%bmort write(fates_log(),fmt0) 'hf_sm_threshold = ',EDPftvarcon_inst%hf_sm_threshold write(fates_log(),fmt0) 'vcmaxha = ',EDPftvarcon_inst%vcmaxha diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 old mode 100755 new mode 100644 diff --git a/main/FatesConstantsMod.F90 b/main/FatesConstantsMod.F90 index b536b3e11b..f7980b4bff 100644 --- a/main/FatesConstantsMod.F90 +++ b/main/FatesConstantsMod.F90 @@ -8,6 +8,7 @@ module FatesConstantsMod ! kinds integer, parameter :: fates_r8 = selected_real_kind(12) ! 8 byte real + integer, parameter :: fates_int = selected_int_kind(8) ! 4 byte int ! string lengths integer, parameter :: fates_avg_flag_length = 3 @@ -37,6 +38,9 @@ module FatesConstantsMod ! Conversion factor: miligrams per grams real(fates_r8), parameter :: mg_per_g = 1000.0_fates_r8 + ! Conversion factor: kilograms per Megagram + real(fates_r8), parameter :: kg_per_Megag = 1000.0_fates_r8 + ! Conversion factor: micromoles per milimole real(fates_r8), parameter :: umol_per_mmol = 1000.0_fates_r8 @@ -49,6 +53,9 @@ module FatesConstantsMod ! Conversion factor: m2 per ha real(fates_r8), parameter :: m2_per_ha = 1.0e4_fates_r8 + ! Conversion factor: cm2 per m2 + real(fates_r8), parameter :: cm2_per_m2 = 10000.0_fates_r8 + ! Conversion factor :: ha per m2 real(fates_r8), parameter :: ha_per_m2 = 1.0e-4_fates_r8 diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index d32732339b..8287713fe6 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -71,6 +71,8 @@ module FatesHistoryInterfaceMod integer, private :: ih_bdead_pa integer, private :: ih_balive_pa integer, private :: ih_bleaf_pa + integer, private :: ih_bsapwood_pa + integer, private :: ih_bfineroot_pa integer, private :: ih_btotal_pa integer, private :: ih_npp_pa integer, private :: ih_gpp_pa @@ -159,8 +161,6 @@ module FatesHistoryInterfaceMod integer, private :: ih_m4_si_scpf integer, private :: ih_m5_si_scpf integer, private :: ih_m6_si_scpf - - !LOGGING , make sure to add ih_m7_si_scpf and hio_m7_si_scpf integer, private :: ih_m7_si_scpf integer, private :: ih_ar_si_scpf @@ -173,6 +173,7 @@ module FatesHistoryInterfaceMod ! indices to (site x scls) variables integer, private :: ih_ba_si_scls + integer, private :: ih_nplant_si_scls integer, private :: ih_nplant_canopy_si_scls integer, private :: ih_nplant_understory_si_scls integer, private :: ih_mortality_canopy_si_scls @@ -185,6 +186,17 @@ module FatesHistoryInterfaceMod integer, private :: ih_crown_area_understory_si_scls integer, private :: ih_ddbh_canopy_si_scls integer, private :: ih_ddbh_understory_si_scls + integer, private :: ih_agb_si_scls + integer, private :: ih_biomass_si_scls + + ! mortality vars + integer, private :: ih_m1_si_scls + integer, private :: ih_m2_si_scls + integer, private :: ih_m3_si_scls + integer, private :: ih_m4_si_scls + integer, private :: ih_m5_si_scls + integer, private :: ih_m6_si_scls + integer, private :: ih_m7_si_scls ! lots of non-default diagnostics for understanding canopy versus understory carbon balances integer, private :: ih_rdark_canopy_si_scls @@ -1170,6 +1182,8 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_bdead_pa => this%hvars(ih_bdead_pa)%r81d, & hio_balive_pa => this%hvars(ih_balive_pa)%r81d, & hio_bleaf_pa => this%hvars(ih_bleaf_pa)%r81d, & + hio_bsapwood_pa => this%hvars(ih_bsapwood_pa)%r81d, & + hio_bfineroot_pa => this%hvars(ih_bfineroot_pa)%r81d, & hio_btotal_pa => this%hvars(ih_btotal_pa)%r81d, & hio_canopy_biomass_pa => this%hvars(ih_canopy_biomass_pa)%r81d, & hio_understory_biomass_pa => this%hvars(ih_understory_biomass_pa)%r81d, & @@ -1208,10 +1222,18 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_m4_si_scpf => this%hvars(ih_m4_si_scpf)%r82d, & hio_m5_si_scpf => this%hvars(ih_m5_si_scpf)%r82d, & hio_m6_si_scpf => this%hvars(ih_m6_si_scpf)%r82d, & - hio_m7_si_scpf => this%hvars(ih_m7_si_scpf)%r82d, & - + hio_m1_si_scls => this%hvars(ih_m1_si_scls)%r82d, & + hio_m2_si_scls => this%hvars(ih_m2_si_scls)%r82d, & + hio_m3_si_scls => this%hvars(ih_m3_si_scls)%r82d, & + hio_m4_si_scls => this%hvars(ih_m4_si_scls)%r82d, & + hio_m5_si_scls => this%hvars(ih_m5_si_scls)%r82d, & + hio_m6_si_scls => this%hvars(ih_m6_si_scls)%r82d, & + hio_m7_si_scls => this%hvars(ih_m7_si_scls)%r82d, & hio_ba_si_scls => this%hvars(ih_ba_si_scls)%r82d, & + hio_agb_si_scls => this%hvars(ih_agb_si_scls)%r82d, & + hio_biomass_si_scls => this%hvars(ih_biomass_si_scls)%r82d, & + hio_nplant_si_scls => this%hvars(ih_nplant_si_scls)%r82d, & hio_nplant_canopy_si_scls => this%hvars(ih_nplant_canopy_si_scls)%r82d, & hio_nplant_understory_si_scls => this%hvars(ih_nplant_understory_si_scls)%r82d, & hio_mortality_canopy_si_scls => this%hvars(ih_mortality_canopy_si_scls)%r82d, & @@ -1376,11 +1398,13 @@ subroutine update_history_dyn(this,nc,nsites,sites) + ccohort%c_area * AREA_INV ! Update biomass components - hio_bleaf_pa(io_pa) = hio_bleaf_pa(io_pa) + n_density * ccohort%bl * g_per_kg - hio_bstore_pa(io_pa) = hio_bstore_pa(io_pa) + n_density * ccohort%bstore * g_per_kg - hio_btotal_pa(io_pa) = hio_btotal_pa(io_pa) + n_density * ccohort%b * g_per_kg - hio_bdead_pa(io_pa) = hio_bdead_pa(io_pa) + n_density * ccohort%bdead * g_per_kg - hio_balive_pa(io_pa) = hio_balive_pa(io_pa) + n_density * ccohort%balive * g_per_kg + hio_bleaf_pa(io_pa) = hio_bleaf_pa(io_pa) + n_density * ccohort%bl * g_per_kg + hio_bstore_pa(io_pa) = hio_bstore_pa(io_pa) + n_density * ccohort%bstore * g_per_kg + hio_btotal_pa(io_pa) = hio_btotal_pa(io_pa) + n_density * ccohort%b * g_per_kg + hio_bdead_pa(io_pa) = hio_bdead_pa(io_pa) + n_density * ccohort%bdead * g_per_kg + hio_balive_pa(io_pa) = hio_balive_pa(io_pa) + n_density * ccohort%balive * g_per_kg + hio_bsapwood_pa(io_pa) = hio_bsapwood_pa(io_pa) + n_density * ccohort%bsw * g_per_kg + hio_bfineroot_pa(io_pa) = hio_bfineroot_pa(io_pa) + n_density * ccohort%br * g_per_kg ! Update PFT partitioned biomass components hio_leafbiomass_si_pft(io_si,ft) = hio_leafbiomass_si_pft(io_si,ft) + & @@ -1416,42 +1440,44 @@ subroutine update_history_dyn(this,nc,nsites,sites) n_perm2*ccohort%gpp_acc_hold ! [kgC/m2/yr] hio_npp_totl_si_scpf(io_si,scpf) = hio_npp_totl_si_scpf(io_si,scpf) + & ccohort%npp_acc_hold *n_perm2 - hio_npp_leaf_si_scpf(io_si,scpf) = hio_npp_leaf_si_scpf(io_si,scpf) + & - ccohort%npp_leaf*n_perm2 - hio_npp_fnrt_si_scpf(io_si,scpf) = hio_npp_fnrt_si_scpf(io_si,scpf) + & - ccohort%npp_froot*n_perm2 - hio_npp_bgsw_si_scpf(io_si,scpf) = hio_npp_bgsw_si_scpf(io_si,scpf) + & - ccohort%npp_bsw*n_perm2* & - (1._r8-EDPftvarcon_inst%allom_agb_frac(ccohort%pft)) - hio_npp_agsw_si_scpf(io_si,scpf) = hio_npp_agsw_si_scpf(io_si,scpf) + & - ccohort%npp_bsw*n_perm2* & - EDPftvarcon_inst%allom_agb_frac(ccohort%pft) - hio_npp_bgdw_si_scpf(io_si,scpf) = hio_npp_bgdw_si_scpf(io_si,scpf) + & - ccohort%npp_bdead*n_perm2* & - (1._r8-EDPftvarcon_inst%allom_agb_frac(ccohort%pft)) - hio_npp_agdw_si_scpf(io_si,scpf) = hio_npp_agdw_si_scpf(io_si,scpf) + & - ccohort%npp_bdead*n_perm2* & - EDPftvarcon_inst%allom_agb_frac(ccohort%pft) - hio_npp_seed_si_scpf(io_si,scpf) = hio_npp_seed_si_scpf(io_si,scpf) + & - ccohort%npp_bseed*n_perm2 - hio_npp_stor_si_scpf(io_si,scpf) = hio_npp_stor_si_scpf(io_si,scpf) + & - ccohort%npp_store*n_perm2 - - if( abs(ccohort%npp_acc_hold-(ccohort%npp_leaf+ccohort%npp_froot+ & - ccohort%npp_bsw+ccohort%npp_bdead+ & - ccohort%npp_bseed+ccohort%npp_store))>1.e-9) then - write(fates_log(),*) 'NPP Partitions are not balancing' - write(fates_log(),*) 'Fractional Error: ', & - abs(ccohort%npp_acc_hold-(ccohort%npp_leaf+ccohort%npp_froot+ & - ccohort%npp_bsw+ccohort%npp_bdead+ & - ccohort%npp_bseed+ccohort%npp_store))/ccohort%npp_acc_hold - write(fates_log(),*) 'Terms: ',ccohort%npp_acc_hold,ccohort%npp_leaf,ccohort%npp_froot, & - ccohort%npp_bsw,ccohort%npp_bdead, & - ccohort%npp_bseed,ccohort%npp_store - write(fates_log(),*) ' NPP components during FATES-HLM linking does not balance ' - stop ! we need termination control for FATES!!! - ! call endrun(msg=errMsg(__FILE__, __LINE__)) - end if +! hio_npp_leaf_si_scpf(io_si,scpf) = hio_npp_leaf_si_scpf(io_si,scpf) + & +! ccohort%npp_leaf*n_perm2 +! hio_npp_fnrt_si_scpf(io_si,scpf) = hio_npp_fnrt_si_scpf(io_si,scpf) + & +! ccohort%npp_froot*n_perm2 +! hio_npp_bgsw_si_scpf(io_si,scpf) = hio_npp_bgsw_si_scpf(io_si,scpf) + & +! ccohort%npp_bsw*n_perm2* & +! (1._r8-EDPftvarcon_inst%allom_agb_frac(ccohort%pft)) +! hio_npp_agsw_si_scpf(io_si,scpf) = hio_npp_agsw_si_scpf(io_si,scpf) + & +! ccohort%npp_bsw*n_perm2* & +! EDPftvarcon_inst%allom_agb_frac(ccohort%pft) +! hio_npp_bgdw_si_scpf(io_si,scpf) = hio_npp_bgdw_si_scpf(io_si,scpf) + & +! ccohort%npp_bdead*n_perm2* & +! (1._r8-EDPftvarcon_inst%allom_agb_frac(ccohort%pft)) +! hio_npp_agdw_si_scpf(io_si,scpf) = hio_npp_agdw_si_scpf(io_si,scpf) + & +! ccohort%npp_bdead*n_perm2* & +! EDPftvarcon_inst%allom_agb_frac(ccohort%pft) +! hio_npp_seed_si_scpf(io_si,scpf) = hio_npp_seed_si_scpf(io_si,scpf) + & +! ccohort%npp_bseed*n_perm2 +! hio_npp_stor_si_scpf(io_si,scpf) = hio_npp_stor_si_scpf(io_si,scpf) + & +! ccohort%npp_store*n_perm2 + + ! TEMPORARILY DISABLING THIS UNTIL THE ALLOCATION REFACTOR IS COMPLETE + ! RGK Feb-2017 +! if( abs(ccohort%npp_acc_hold-(ccohort%npp_leaf+ccohort%npp_froot+ & +! ccohort%npp_bsw+ccohort%npp_bdead+ & +! ccohort%npp_bseed+ccohort%npp_store))>1.e-9) then +! write(fates_log(),*) 'NPP Partitions are not balancing' +! write(fates_log(),*) 'Fractional Error: ', & +! abs(ccohort%npp_acc_hold-(ccohort%npp_leaf+ccohort%npp_froot+ & +! ccohort%npp_bsw+ccohort%npp_bdead+ & +! ccohort%npp_bseed+ccohort%npp_store))/ccohort%npp_acc_hold +! write(fates_log(),*) 'Terms: ',ccohort%npp_acc_hold,ccohort%npp_leaf,ccohort%npp_froot, & +! ccohort%npp_bsw,ccohort%npp_bdead, & +! ccohort%npp_bseed,ccohort%npp_store +! write(fates_log(),*) ' NPP components during FATES-HLM linking does not balance ' +! stop ! we need termination control for FATES!!! +! ! call endrun(msg=errMsg(__FILE__, __LINE__)) +! end if ! Woody State Variables (basal area and number density and mortality) if (EDPftvarcon_inst%woody(ft) == 1) then @@ -1460,12 +1486,16 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_m2_si_scpf(io_si,scpf) = hio_m2_si_scpf(io_si,scpf) + ccohort%hmort*ccohort%n hio_m3_si_scpf(io_si,scpf) = hio_m3_si_scpf(io_si,scpf) + ccohort%cmort*ccohort%n hio_m5_si_scpf(io_si,scpf) = hio_m5_si_scpf(io_si,scpf) + ccohort%fmort*ccohort%n - - - !Y.X. hio_m7_si_scpf(io_si,scpf) = hio_m7_si_scpf(io_si,scpf) + & (ccohort%lmort_logging+ccohort%lmort_collateral+ccohort%lmort_infra) * ccohort%n + hio_m1_si_scls(io_si,scls) = hio_m1_si_scls(io_si,scls) + ccohort%bmort*ccohort%n + hio_m2_si_scls(io_si,scls) = hio_m2_si_scls(io_si,scls) + ccohort%hmort*ccohort%n + hio_m3_si_scls(io_si,scls) = hio_m3_si_scls(io_si,scls) + ccohort%cmort*ccohort%n + hio_m5_si_scls(io_si,scls) = hio_m5_si_scls(io_si,scls) + ccohort%fmort*ccohort%n + hio_m7_si_scls(io_si,scls) = hio_m7_si_scls(io_si,scls) + & + (ccohort%lmort_logging+ccohort%lmort_collateral+ccohort%lmort_infra) * ccohort%n + ! basal area [m2/ha] hio_ba_si_scpf(io_si,scpf) = hio_ba_si_scpf(io_si,scpf) + & @@ -1482,12 +1512,20 @@ subroutine update_history_dyn(this,nc,nsites,sites) ccohort%ddbhdt*ccohort%n end if + hio_agb_si_scls(io_si,scls) = hio_agb_si_scls(io_si,scls) + & + ccohort%b * ccohort%n * EDPftvarcon_inst%allom_agb_frac(ccohort%pft) * AREA_INV + + hio_biomass_si_scls(io_si,scls) = hio_biomass_si_scls(io_si,scls) + & + ccohort%b * ccohort%n * AREA_INV + ! update size-class x patch-age related quantities iscag = get_sizeage_class_index(ccohort%dbh,cpatch%age) hio_nplant_si_scag(io_si,iscag) = hio_nplant_si_scag(io_si,iscag) + ccohort%n + hio_nplant_si_scls(io_si,scls) = hio_nplant_si_scls(io_si,scls) + ccohort%n + ! update SCPF/SCLS- and canopy/subcanopy- partitioned quantities if (ccohort%canopy_layer .eq. 1) then hio_nplant_canopy_si_scag(io_si,iscag) = hio_nplant_canopy_si_scag(io_si,iscag) + ccohort%n @@ -1554,18 +1592,19 @@ subroutine update_history_dyn(this,nc,nsites,sites) ccohort%dbstoredt * ccohort%n hio_storage_flux_canopy_si_scls(io_si,scls) = hio_storage_flux_canopy_si_scls(io_si,scls) + & ccohort%storage_flux * ccohort%n - hio_npp_leaf_canopy_si_scls(io_si,scls) = hio_npp_leaf_canopy_si_scls(io_si,scls) + & - ccohort%npp_leaf * ccohort%n - hio_npp_froot_canopy_si_scls(io_si,scls) = hio_npp_froot_canopy_si_scls(io_si,scls) + & - ccohort%npp_froot * ccohort%n - hio_npp_bsw_canopy_si_scls(io_si,scls) = hio_npp_bsw_canopy_si_scls(io_si,scls) + & - ccohort%npp_bsw * ccohort%n - hio_npp_bdead_canopy_si_scls(io_si,scls) = hio_npp_bdead_canopy_si_scls(io_si,scls) + & - ccohort%npp_bdead * ccohort%n - hio_npp_bseed_canopy_si_scls(io_si,scls) = hio_npp_bseed_canopy_si_scls(io_si,scls) + & - ccohort%npp_bseed * ccohort%n - hio_npp_store_canopy_si_scls(io_si,scls) = hio_npp_store_canopy_si_scls(io_si,scls) + & - ccohort%npp_store * ccohort%n + +! hio_npp_leaf_canopy_si_scls(io_si,scls) = hio_npp_leaf_canopy_si_scls(io_si,scls) + & +! ccohort%npp_leaf * ccohort%n +! hio_npp_froot_canopy_si_scls(io_si,scls) = hio_npp_froot_canopy_si_scls(io_si,scls) + & +! ccohort%npp_froot * ccohort%n +! hio_npp_bsw_canopy_si_scls(io_si,scls) = hio_npp_bsw_canopy_si_scls(io_si,scls) + & +! ccohort%npp_bsw * ccohort%n +! hio_npp_bdead_canopy_si_scls(io_si,scls) = hio_npp_bdead_canopy_si_scls(io_si,scls) + & +! ccohort%npp_bdead * ccohort%n +! hio_npp_bseed_canopy_si_scls(io_si,scls) = hio_npp_bseed_canopy_si_scls(io_si,scls) + & +! ccohort%npp_bseed * ccohort%n +! hio_npp_store_canopy_si_scls(io_si,scls) = hio_npp_store_canopy_si_scls(io_si,scls) + & +! ccohort%npp_store * ccohort%n hio_yesterdaycanopylevel_canopy_si_scls(io_si,scls) = & hio_yesterdaycanopylevel_canopy_si_scls(io_si,scls) + & ccohort%canopy_layer_yesterday * ccohort%n @@ -1633,18 +1672,18 @@ subroutine update_history_dyn(this,nc,nsites,sites) ccohort%dbstoredt * ccohort%n hio_storage_flux_understory_si_scls(io_si,scls) = hio_storage_flux_understory_si_scls(io_si,scls) + & ccohort%storage_flux * ccohort%n - hio_npp_leaf_understory_si_scls(io_si,scls) = hio_npp_leaf_understory_si_scls(io_si,scls) + & - ccohort%npp_leaf * ccohort%n - hio_npp_froot_understory_si_scls(io_si,scls) = hio_npp_froot_understory_si_scls(io_si,scls) + & - ccohort%npp_froot * ccohort%n - hio_npp_bsw_understory_si_scls(io_si,scls) = hio_npp_bsw_understory_si_scls(io_si,scls) + & - ccohort%npp_bsw * ccohort%n - hio_npp_bdead_understory_si_scls(io_si,scls) = hio_npp_bdead_understory_si_scls(io_si,scls) + & - ccohort%npp_bdead * ccohort%n - hio_npp_bseed_understory_si_scls(io_si,scls) = hio_npp_bseed_understory_si_scls(io_si,scls) + & - ccohort%npp_bseed * ccohort%n - hio_npp_store_understory_si_scls(io_si,scls) = hio_npp_store_understory_si_scls(io_si,scls) + & - ccohort%npp_store * ccohort%n +! hio_npp_leaf_understory_si_scls(io_si,scls) = hio_npp_leaf_understory_si_scls(io_si,scls) + & +! ccohort%npp_leaf * ccohort%n +! hio_npp_froot_understory_si_scls(io_si,scls) = hio_npp_froot_understory_si_scls(io_si,scls) + & +! ccohort%npp_froot * ccohort%n +! hio_npp_bsw_understory_si_scls(io_si,scls) = hio_npp_bsw_understory_si_scls(io_si,scls) + & +! ccohort%npp_bsw * ccohort%n +! hio_npp_bdead_understory_si_scls(io_si,scls) = hio_npp_bdead_understory_si_scls(io_si,scls) + & +! ccohort%npp_bdead * ccohort%n +! hio_npp_bseed_understory_si_scls(io_si,scls) = hio_npp_bseed_understory_si_scls(io_si,scls) + & +! ccohort%npp_bseed * ccohort%n +! hio_npp_store_understory_si_scls(io_si,scls) = hio_npp_store_understory_si_scls(io_si,scls) + & +! ccohort%npp_store * ccohort%n hio_yesterdaycanopylevel_understory_si_scls(io_si,scls) = & hio_yesterdaycanopylevel_understory_si_scls(io_si,scls) + & ccohort%canopy_layer_yesterday * ccohort%n @@ -1752,8 +1791,15 @@ subroutine update_history_dyn(this,nc,nsites,sites) do i_pft = 1, numpft do i_scls = 1,nlevsclass i_scpf = (i_pft-1)*nlevsclass + i_scls + ! + ! termination mortality. sum of canopy and understory indices hio_m6_si_scpf(io_si,i_scpf) = (sites(s)%terminated_nindivs(i_scls,i_pft,1) + & sites(s)%terminated_nindivs(i_scls,i_pft,2)) * days_per_year + hio_m6_si_scls(io_si,i_scls) = hio_m6_si_scls(io_si,i_scls) + & + (sites(s)%terminated_nindivs(i_scls,i_pft,1) + & + sites(s)%terminated_nindivs(i_scls,i_pft,2)) * days_per_year + ! + ! add termination mortality to canopy and understory mortality hio_mortality_canopy_si_scls(io_si,i_scls) = hio_mortality_canopy_si_scls(io_si,i_scls) + & sites(s)%terminated_nindivs(i_scls,i_pft,1) * days_per_year hio_mortality_understory_si_scls(io_si,i_scls) = hio_mortality_understory_si_scls(io_si,i_scls) + & @@ -1764,7 +1810,8 @@ subroutine update_history_dyn(this,nc,nsites,sites) sites(s)%terminated_nindivs(i_scls,i_pft,2) * days_per_year ! ! imort on its own - hio_m4_si_scpf(io_si,i_scpf) = hio_m4_si_scpf(io_si,i_scpf) + sites(s)%imort_rate(i_scls, i_pft) + hio_m4_si_scpf(io_si,i_scpf) = sites(s)%imort_rate(i_scls, i_pft) + hio_m4_si_scls(io_si,i_scls) = hio_m4_si_scls(io_si,i_scls) + sites(s)%imort_rate(i_scls, i_pft) ! ! add imort to other mortality terms. consider imort as understory mortality even if it happens in ! cohorts that may have been promoted as part of the patch creation, and use the pre-calculated site-level @@ -2759,6 +2806,16 @@ subroutine define_history_vars(this, initialize_variables) avgflag='A', vtype=patch_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, & ivar=ivar, initialize=initialize_variables, index = ih_bleaf_pa ) + call this%set_history_var(vname='ED_bsapwood', units='gC m-2', & + long='Sapwood biomass', use_default='active', & + avgflag='A', vtype=patch_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, & + ivar=ivar, initialize=initialize_variables, index = ih_bsapwood_pa ) + + call this%set_history_var(vname='ED_bfineroot', units='gC m-2', & + long='Fine root biomass', use_default='active', & + avgflag='A', vtype=patch_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, & + ivar=ivar, initialize=initialize_variables, index = ih_bfineroot_pa ) + call this%set_history_var(vname='ED_biomass', units='gC m-2', & long='Total biomass', use_default='active', & avgflag='A', vtype=patch_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, & @@ -3022,18 +3079,18 @@ subroutine define_history_vars(this, initialize_variables) upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_promotion_carbonflux_si ) call this%set_history_var(vname='MORTALITY_CARBONFLUX_CANOPY', units = 'gC/m2/s', & - long='flux of biomass carbon from live to dead pools from mortality of canopy plants', use_default='inactive', & + long='flux of biomass carbon from live to dead pools from mortality of canopy plants', use_default='active', & avgflag='A', vtype=site_r8, hlms='CLM:ALM', flushval=0.0_r8, & upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_canopy_mortality_carbonflux_si ) call this%set_history_var(vname='MORTALITY_CARBONFLUX_UNDERSTORY', units = 'gC/m2/s', & - long='flux of biomass carbon from live to dead pools from mortality of understory plants',use_default='inactive',& + long='flux of biomass carbon from live to dead pools from mortality of understory plants',use_default='active',& avgflag='A', vtype=site_r8, hlms='CLM:ALM', flushval=0.0_r8, & upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_understory_mortality_carbonflux_si ) - + ! size class by age dimensioned variables call this%set_history_var(vname='NPLANT_SCAG',units = 'plants/ha', & - long='number of plants per hectare in each size x age class', use_default='inactive', & + long='number of plants per hectare in each size x age class', use_default='active', & avgflag='A', vtype=site_scag_r8, hlms='CLM:ALM', flushval=0.0_r8, & upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_nplant_si_scag ) @@ -3103,45 +3160,45 @@ subroutine define_history_vars(this, initialize_variables) upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_npp_totl_si_scpf ) - call this%set_history_var(vname='NPP_LEAF_SCPF', units='kgC/m2/yr', & - long='NPP flux into leaves by pft/size', use_default='inactive', & - avgflag='A', vtype=site_size_pft_r8, hlms='CLM:ALM', flushval=0.0_r8, & - upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_npp_leaf_si_scpf ) + call this%set_history_var(vname='NPP_LEAF_SCPF', units='kgC/m2/yr', & + long='NPP flux into leaves by pft/size', use_default='inactive', & + avgflag='A', vtype=site_size_pft_r8, hlms='CLM:ALM', flushval=-999.9_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_npp_leaf_si_scpf ) - call this%set_history_var(vname='NPP_SEED_SCPF', units='kgC/m2/yr', & - long='NPP flux into seeds by pft/size', use_default='inactive', & - avgflag='A', vtype=site_size_pft_r8, hlms='CLM:ALM', flushval=0.0_r8, & - upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_npp_seed_si_scpf ) + call this%set_history_var(vname='NPP_SEED_SCPF', units='kgC/m2/yr', & + long='NPP flux into seeds by pft/size', use_default='inactive', & + avgflag='A', vtype=site_size_pft_r8, hlms='CLM:ALM', flushval=-999.9_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_npp_seed_si_scpf ) - call this%set_history_var(vname='NPP_FNRT_SCPF', units='kgC/m2/yr', & - long='NPP flux into fine roots by pft/size', use_default='inactive', & - avgflag='A', vtype=site_size_pft_r8, hlms='CLM:ALM', flushval=0.0_r8, & - upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_npp_fnrt_si_scpf ) + call this%set_history_var(vname='NPP_FNRT_SCPF', units='kgC/m2/yr', & + long='NPP flux into fine roots by pft/size', use_default='inactive', & + avgflag='A', vtype=site_size_pft_r8, hlms='CLM:ALM', flushval=-999.9_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_npp_fnrt_si_scpf ) - call this%set_history_var(vname='NPP_BGSW_SCPF', units='kgC/m2/yr', & - long='NPP flux into below-ground sapwood by pft/size', use_default='inactive', & - avgflag='A', vtype=site_size_pft_r8, hlms='CLM:ALM', flushval=0.0_r8, & - upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_npp_bgsw_si_scpf ) + call this%set_history_var(vname='NPP_BGSW_SCPF', units='kgC/m2/yr', & + long='NPP flux into below-ground sapwood by pft/size', use_default='inactive', & + avgflag='A', vtype=site_size_pft_r8, hlms='CLM:ALM', flushval=-999.9_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_npp_bgsw_si_scpf ) - call this%set_history_var(vname='NPP_BGDW_SCPF', units='kgC/m2/yr', & - long='NPP flux into below-ground deadwood by pft/size', use_default='inactive', & - avgflag='A', vtype=site_size_pft_r8, hlms='CLM:ALM', flushval=0.0_r8, & - upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_npp_bgdw_si_scpf ) + call this%set_history_var(vname='NPP_BGDW_SCPF', units='kgC/m2/yr', & + long='NPP flux into below-ground deadwood by pft/size', use_default='inactive', & + avgflag='A', vtype=site_size_pft_r8, hlms='CLM:ALM', flushval=-999.9_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_npp_bgdw_si_scpf ) - call this%set_history_var(vname='NPP_AGSW_SCPF', units='kgC/m2/yr', & - long='NPP flux into above-ground sapwood by pft/size', use_default='inactive', & - avgflag='A', vtype=site_size_pft_r8, hlms='CLM:ALM', flushval=0.0_r8, & - upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_npp_agsw_si_scpf ) + call this%set_history_var(vname='NPP_AGSW_SCPF', units='kgC/m2/yr', & + long='NPP flux into above-ground sapwood by pft/size', use_default='inactive', & + avgflag='A', vtype=site_size_pft_r8, hlms='CLM:ALM', flushval=-999.9_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_npp_agsw_si_scpf ) - call this%set_history_var(vname = 'NPP_AGDW_SCPF', units='kgC/m2/yr', & - long='NPP flux into above-ground deadwood by pft/size', use_default='inactive', & - avgflag='A', vtype=site_size_pft_r8, hlms='CLM:ALM', flushval=0.0_r8, & - upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_npp_agdw_si_scpf ) + call this%set_history_var(vname = 'NPP_AGDW_SCPF', units='kgC/m2/yr', & + long='NPP flux into above-ground deadwood by pft/size', use_default='inactive', & + avgflag='A', vtype=site_size_pft_r8, hlms='CLM:ALM', flushval=-999.9_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_npp_agdw_si_scpf ) - call this%set_history_var(vname = 'NPP_STOR_SCPF', units='kgC/m2/yr', & - long='NPP flux into storage by pft/size', use_default='inactive', & - avgflag='A', vtype=site_size_pft_r8, hlms='CLM:ALM', flushval=0.0_r8, & - upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_npp_stor_si_scpf ) + call this%set_history_var(vname = 'NPP_STOR_SCPF', units='kgC/m2/yr', & + long='NPP flux into storage by pft/size', use_default='inactive', & + avgflag='A', vtype=site_size_pft_r8, hlms='CLM:ALM', flushval=-999.9_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_npp_stor_si_scpf ) call this%set_history_var(vname='DDBH_SCPF', units = 'cm/yr/ha', & long='diameter growth increment by pft/size',use_default='inactive', & @@ -3198,14 +3255,11 @@ subroutine define_history_vars(this, initialize_variables) avgflag='A', vtype=site_size_pft_r8, hlms='CLM:ALM', flushval=0.0_r8, & upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_m6_si_scpf ) - - !Logging call this%set_history_var(vname='M7_SCPF', units = 'N/ha/event', & - long='logging mortalities by pft/size',use_default='inactive', & + long='logging mortality by pft/size',use_default='inactive', & avgflag='A', vtype=site_size_pft_r8, hlms='CLM:ALM', flushval=0.0_r8, & upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_m7_si_scpf ) - call this%set_history_var(vname='MORTALITY_CANOPY_SCPF', units = 'N/ha/yr', & long='total mortality of canopy plants by pft/size', use_default='inactive', & avgflag='A', vtype=site_size_pft_r8, hlms='CLM:ALM', flushval=0.0_r8, & @@ -3316,12 +3370,12 @@ subroutine define_history_vars(this, initialize_variables) ! size-class only variables call this%set_history_var(vname='DDBH_CANOPY_SCLS', units = 'cm/yr/ha', & - long='diameter growth increment by pft/size',use_default='inactive', & + long='diameter growth increment by pft/size',use_default='active', & avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_ddbh_canopy_si_scls ) call this%set_history_var(vname='DDBH_UNDERSTORY_SCLS', units = 'cm/yr/ha', & - long='diameter growth increment by pft/size',use_default='inactive', & + long='diameter growth increment by pft/size',use_default='active', & avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_ddbh_understory_si_scls ) @@ -3336,10 +3390,20 @@ subroutine define_history_vars(this, initialize_variables) upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_yesterdaycanopylevel_understory_si_scls ) call this%set_history_var(vname='BA_SCLS', units = 'm2/ha', & - long='basal area by size class', use_default='inactive', & + long='basal area by size class', use_default='active', & avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_ba_si_scls ) + call this%set_history_var(vname='AGB_SCLS', units = 'kgC/m2', & + long='Aboveground biomass by size class', use_default='active', & + avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_agb_si_scls ) + + call this%set_history_var(vname='BIOMASS_SCLS', units = 'kgC/m2', & + long='Total biomass by size class', use_default='inactive', & + avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_biomass_si_scls ) + call this%set_history_var(vname='DEMOTION_RATE_SCLS', units = 'indiv/ha/yr', & long='demotion rate from canopy to understory by size class', use_default='inactive', & avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & @@ -3351,22 +3415,62 @@ subroutine define_history_vars(this, initialize_variables) upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_promotion_rate_si_scls ) call this%set_history_var(vname='NPLANT_CANOPY_SCLS', units = 'indiv/ha', & - long='number of canopy plants by size class', use_default='inactive', & + long='number of canopy plants by size class', use_default='active', & avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_nplant_canopy_si_scls ) call this%set_history_var(vname='MORTALITY_CANOPY_SCLS', units = 'indiv/ha/yr', & - long='total mortality of canopy trees by size class', use_default='inactive', & + long='total mortality of canopy trees by size class', use_default='active', & avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_mortality_canopy_si_scls ) call this%set_history_var(vname='NPLANT_UNDERSTORY_SCLS', units = 'indiv/ha', & - long='number of understory plants by size class', use_default='inactive', & + long='number of understory plants by size class', use_default='active', & avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_nplant_understory_si_scls ) + call this%set_history_var(vname='NPLANT_SCLS', units = 'indiv/ha', & + long='number of plants by size class', use_default='active', & + avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_nplant_si_scls ) + + call this%set_history_var(vname='M1_SCLS', units = 'N/ha/yr', & + long='background mortality by size', use_default='active', & + avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_m1_si_scls ) + + call this%set_history_var(vname='M2_SCLS', units = 'N/ha/yr', & + long='hydraulic mortality by size',use_default='active', & + avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_m2_si_scls ) + + call this%set_history_var(vname='M3_SCLS', units = 'N/ha/yr', & + long='carbon starvation mortality by size', use_default='active', & + avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_m3_si_scls ) + + call this%set_history_var(vname='M4_SCLS', units = 'N/ha/yr', & + long='impact mortality by size',use_default='active', & + avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_m4_si_scls ) + + call this%set_history_var(vname='M5_SCLS', units = 'N/ha/yr', & + long='fire mortality by size',use_default='active', & + avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_m5_si_scls ) + + call this%set_history_var(vname='M6_SCLS', units = 'N/ha/yr', & + long='termination mortality by size',use_default='active', & + avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_m6_si_scls ) + + call this%set_history_var(vname='M7_SCLS', units = 'N/ha/event', & + long='logging mortality by size',use_default='active', & + avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_m7_si_scls ) + call this%set_history_var(vname='MORTALITY_UNDERSTORY_SCLS', units = 'indiv/ha/yr', & - long='total mortality of understory trees by size class', use_default='inactive', & + long='total mortality of understory trees by size class', use_default='active', & avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_mortality_understory_si_scls ) @@ -3430,35 +3534,35 @@ subroutine define_history_vars(this, initialize_variables) avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_storage_flux_canopy_si_scls ) - call this%set_history_var(vname='NPP_LEAF_CANOPY_SCLS', units = 'kg C / ha / yr', & - long='NPP_LEAF for canopy plants by size class', use_default='inactive', & - avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & - upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_npp_leaf_canopy_si_scls ) + call this%set_history_var(vname='NPP_LEAF_CANOPY_SCLS', units = 'kg C / ha / yr', & + long='NPP_LEAF for canopy plants by size class', use_default='inactive', & + avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=-999.9_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_npp_leaf_canopy_si_scls ) - call this%set_history_var(vname='NPP_FROOT_CANOPY_SCLS', units = 'kg C / ha / yr', & - long='NPP_FROOT for canopy plants by size class', use_default='inactive', & - avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & - upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_npp_froot_canopy_si_scls ) + call this%set_history_var(vname='NPP_FROOT_CANOPY_SCLS', units = 'kg C / ha / yr', & + long='NPP_FROOT for canopy plants by size class', use_default='inactive', & + avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=-999.9_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_npp_froot_canopy_si_scls ) - call this%set_history_var(vname='NPP_BSW_CANOPY_SCLS', units = 'kg C / ha / yr', & - long='NPP_BSW for canopy plants by size class', use_default='inactive', & - avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & - upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_npp_bsw_canopy_si_scls ) + call this%set_history_var(vname='NPP_BSW_CANOPY_SCLS', units = 'kg C / ha / yr', & + long='NPP_BSW for canopy plants by size class', use_default='inactive', & + avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=-999.9_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_npp_bsw_canopy_si_scls ) - call this%set_history_var(vname='NPP_BDEAD_CANOPY_SCLS', units = 'kg C / ha / yr', & - long='NPP_BDEAD for canopy plants by size class', use_default='inactive', & - avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & - upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_npp_bdead_canopy_si_scls ) + call this%set_history_var(vname='NPP_BDEAD_CANOPY_SCLS', units = 'kg C / ha / yr', & + long='NPP_BDEAD for canopy plants by size class', use_default='inactive', & + avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=-999.9_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_npp_bdead_canopy_si_scls ) - call this%set_history_var(vname='NPP_BSEED_CANOPY_SCLS', units = 'kg C / ha / yr', & - long='NPP_BSEED for canopy plants by size class', use_default='inactive', & - avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & - upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_npp_bseed_canopy_si_scls ) + call this%set_history_var(vname='NPP_BSEED_CANOPY_SCLS', units = 'kg C / ha / yr', & + long='NPP_BSEED for canopy plants by size class', use_default='inactive', & + avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=-999.9_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_npp_bseed_canopy_si_scls ) - call this%set_history_var(vname='NPP_STORE_CANOPY_SCLS', units = 'kg C / ha / yr', & - long='NPP_STORE for canopy plants by size class', use_default='inactive', & - avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & - upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_npp_store_canopy_si_scls ) + call this%set_history_var(vname='NPP_STORE_CANOPY_SCLS', units = 'kg C / ha / yr', & + long='NPP_STORE for canopy plants by size class', use_default='inactive', & + avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=-999.9_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_npp_store_canopy_si_scls ) call this%set_history_var(vname='RDARK_CANOPY_SCLS', units = 'kg C / ha / yr', & long='RDARK for canopy plants by size class', use_default='inactive', & @@ -3530,35 +3634,35 @@ subroutine define_history_vars(this, initialize_variables) avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_storage_flux_understory_si_scls ) - call this%set_history_var(vname='NPP_LEAF_UNDERSTORY_SCLS', units = 'kg C / ha / yr', & - long='NPP_LEAF for understory plants by size class', use_default='inactive', & - avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & - upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_npp_leaf_understory_si_scls ) - - call this%set_history_var(vname='NPP_FROOT_UNDERSTORY_SCLS', units = 'kg C / ha / yr', & - long='NPP_FROOT for understory plants by size class', use_default='inactive', & - avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & - upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_npp_froot_understory_si_scls ) - - call this%set_history_var(vname='NPP_BSW_UNDERSTORY_SCLS', units = 'kg C / ha / yr', & - long='NPP_BSW for understory plants by size class', use_default='inactive', & - avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & - upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_npp_bsw_understory_si_scls ) - - call this%set_history_var(vname='NPP_BDEAD_UNDERSTORY_SCLS', units = 'kg C / ha / yr', & - long='NPP_BDEAD for understory plants by size class', use_default='inactive', & - avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & - upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_npp_bdead_understory_si_scls ) - - call this%set_history_var(vname='NPP_BSEED_UNDERSTORY_SCLS', units = 'kg C / ha / yr', & - long='NPP_BSEED for understory plants by size class', use_default='inactive', & - avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & - upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_npp_bseed_understory_si_scls ) - - call this%set_history_var(vname='NPP_STORE_UNDERSTORY_SCLS', units = 'kg C / ha / yr', & - long='NPP_STORE for understory plants by size class', use_default='inactive', & - avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & - upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_npp_store_understory_si_scls ) + call this%set_history_var(vname='NPP_LEAF_UNDERSTORY_SCLS', units = 'kg C / ha / yr', & + long='NPP_LEAF for understory plants by size class', use_default='inactive', & + avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=-999.9_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_npp_leaf_understory_si_scls ) + + call this%set_history_var(vname='NPP_FROOT_UNDERSTORY_SCLS', units = 'kg C / ha / yr', & + long='NPP_FROOT for understory plants by size class', use_default='inactive', & + avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=-999.9_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_npp_froot_understory_si_scls ) + + call this%set_history_var(vname='NPP_BSW_UNDERSTORY_SCLS', units = 'kg C / ha / yr', & + long='NPP_BSW for understory plants by size class', use_default='inactive', & + avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=-999.9_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_npp_bsw_understory_si_scls ) + + call this%set_history_var(vname='NPP_BDEAD_UNDERSTORY_SCLS', units = 'kg C / ha / yr', & + long='NPP_BDEAD for understory plants by size class', use_default='inactive', & + avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=-999.9_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_npp_bdead_understory_si_scls ) + + call this%set_history_var(vname='NPP_BSEED_UNDERSTORY_SCLS', units = 'kg C / ha / yr', & + long='NPP_BSEED for understory plants by size class', use_default='inactive', & + avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=-999.9_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_npp_bseed_understory_si_scls ) + + call this%set_history_var(vname='NPP_STORE_UNDERSTORY_SCLS', units = 'kg C / ha / yr', & + long='NPP_STORE for understory plants by size class', use_default='inactive', & + avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=-999.9_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_npp_store_understory_si_scls ) call this%set_history_var(vname='RDARK_UNDERSTORY_SCLS', units = 'kg C / ha / yr', & long='RDARK for understory plants by size class', use_default='inactive', & diff --git a/main/FatesInventoryInitMod.F90 b/main/FatesInventoryInitMod.F90 index d01d9c3d0d..ff8010a281 100644 --- a/main/FatesInventoryInitMod.F90 +++ b/main/FatesInventoryInitMod.F90 @@ -748,9 +748,15 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & ! avgRG (cm/yr?) Average Radial Growth (NOT USED) ! -------------------------------------------------------------------------------------------- - use EDGrowthFunctionsMod, only : hite - use EDGrowthFunctionsMod, only : bleaf - use EDGrowthFunctionsMod, only : bdead + use FatesAllometryMod , only : h_allom + use FatesAllometryMod , only : h2d_allom + use FatesAllometryMod , only : bag_allom + use FatesAllometryMod , only : bcr_allom + use FatesAllometryMod , only : bleaf + use FatesAllometryMod , only : bfineroot + use FatesAllometryMod , only : bsap_allom + use FatesAllometryMod , only : bdead_allom + use EDCohortDynamicsMod , only : create_cohort use FatesInterfaceMod , only : numpft @@ -779,6 +785,11 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & type(ed_cohort_type), pointer :: temp_cohort ! temporary patch (needed for allom funcs) integer :: ipa ! patch idex logical :: matched_patch ! check if cohort was matched w/ patch + real(r8) :: b_ag ! biomass above ground [kgC] + real(r8) :: b_cr ! biomass in coarse roots [kgC] + real(r8) :: b_leaf ! biomass in leaves [kgC] + real(r8) :: b_fineroot ! biomass in fine roots [kgC] + real(r8) :: b_sapwood ! biomass in sapwood [kgC] character(len=128),parameter :: wr_fmt = & '(F7.1,2X,A20,2X,A20,2X,F5.2,2X,F5.2,2X,I4,2X,F5.2,2X,F5.2,2X,F5.2,2X,F5.2)' @@ -857,31 +868,47 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & allocate(temp_cohort) ! A temporary cohort is needed because we want to make ! use of the allometry functions - - temp_cohort%pft = c_pft temp_cohort%n = c_nplant * cpatch%area - temp_cohort%hite = Hite(temp_cohort) temp_cohort%dbh = c_dbh + call h_allom(c_dbh,c_pft,temp_cohort%hite) temp_cohort%canopy_trim = 1.0_r8 - temp_cohort%bdead = Bdead(temp_cohort) - temp_cohort%balive = Bleaf(temp_cohort)*(1.0_r8 + EDPftvarcon_inst%allom_l2fr(c_pft) & - + EDpftvarcon_inst%allom_latosa_int(c_pft)*temp_cohort%hite) + + ! Calculate total above-ground biomass from allometry + + call bag_allom(temp_cohort%dbh,temp_cohort%hite,c_pft,b_ag) + ! Calculate coarse root biomass from allometry + call bcr_allom(temp_cohort%dbh,temp_cohort%hite,c_pft,b_cr) + + ! Calculate the leaf biomass (calculates a maximum first, then applies canopy trim + ! and sla scaling factors) + call bleaf(temp_cohort%dbh,temp_cohort%hite,c_pft,temp_cohort%canopy_trim,b_leaf) + + ! Calculate fine root biomass + call bfineroot(temp_cohort%dbh,temp_cohort%hite,c_pft,temp_cohort%canopy_trim,b_fineroot) + + ! Calculate sapwood biomass + call bsap_allom(temp_cohort%dbh,c_pft,temp_cohort%canopy_trim,b_sapwood) + + temp_cohort%balive = b_leaf + b_fineroot + b_sapwood + + call bdead_allom( b_ag, b_cr, b_sapwood, c_pft, temp_cohort%bdead ) + temp_cohort%b = temp_cohort%balive + temp_cohort%bdead if( EDPftvarcon_inst%evergreen(c_pft) == 1) then - temp_cohort%bstore = Bleaf(temp_cohort) * EDPftvarcon_inst%cushion(c_pft) + temp_cohort%bstore = b_leaf * EDPftvarcon_inst%cushion(c_pft) temp_cohort%laimemory = 0._r8 cstatus = 2 endif if( EDPftvarcon_inst%season_decid(c_pft) == 1 ) then !for dorment places - temp_cohort%bstore = Bleaf(temp_cohort) * EDPftvarcon_inst%cushion(c_pft) !stored carbon in new seedlings. + temp_cohort%bstore = b_leaf * EDPftvarcon_inst%cushion(c_pft) !stored carbon in new seedlings. if(csite%status == 2)then temp_cohort%laimemory = 0.0_r8 else - temp_cohort%laimemory = Bleaf(temp_cohort) + temp_cohort%laimemory = b_leaf endif ! reduce biomass according to size of store, this will be recovered when elaves com on. temp_cohort%balive = temp_cohort%balive - temp_cohort%laimemory @@ -889,8 +916,8 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & endif if ( EDPftvarcon_inst%stress_decid(c_pft) == 1 ) then - temp_cohort%bstore = Bleaf(temp_cohort) * EDPftvarcon_inst%cushion(c_pft) - temp_cohort%laimemory = Bleaf(temp_cohort) + temp_cohort%bstore = b_leaf * EDPftvarcon_inst%cushion(c_pft) + temp_cohort%laimemory = b_leaf temp_cohort%balive = temp_cohort%balive - temp_cohort%laimemory cstatus = csite%dstatus endif diff --git a/main/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index 5f2247614f..7ccbe0bdab 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -1347,12 +1347,12 @@ subroutine create_patchcohort_structure(this, nc, nsites, sites, bc_in) use EDTypesMod, only : maxpft use EDTypesMod, only : area use EDPatchDynamicsMod, only : zero_patch - use EDGrowthFunctionsMod, only : Dbh use EDCohortDynamicsMod, only : create_cohort use EDInitMod, only : zero_site use EDInitMod, only : init_site_vars use EDPatchDynamicsMod, only : create_patch - use EDPftvarcon, only : EDPftvarcon_inst + use EDPftvarcon, only : EDPftvarcon_inst + use FatesAllometryMod, only : h2d_allom ! !ARGUMENTS: class(fates_restart_interface_type) , intent(inout) :: this @@ -1450,7 +1450,7 @@ subroutine create_patchcohort_structure(this, nc, nsites, sites, bc_in) temp_cohort%bdead = 0.0_r8 temp_cohort%bstore = 0.0_r8 temp_cohort%laimemory = 0.0_r8 - temp_cohort%canopy_trim = 0.0_r8 + temp_cohort%canopy_trim = 1.0_r8 temp_cohort%canopy_layer = 1.0_r8 temp_cohort%canopy_layer_yesterday = 1.0_r8 @@ -1469,10 +1469,9 @@ subroutine create_patchcohort_structure(this, nc, nsites, sites, bc_in) endif temp_cohort%hite = 1.25_r8 - ! the dbh function should only take as an argument, the one - ! item it needs, not the entire cohort...refactor - temp_cohort%dbh = Dbh(temp_cohort) + 0.0001_r8*ft - + ! Solve for diameter from height + call h2d_allom(temp_cohort%hite,ft,temp_cohort%dbh) + if (DEBUG) then write(fates_log(),*) 'EDRestVectorMod.F90::createPatchCohortStructure call create_cohort ' end if diff --git a/tools/FatesPFTIndexSwapper.py b/tools/FatesPFTIndexSwapper.py new file mode 100755 index 0000000000..af901b3d18 --- /dev/null +++ b/tools/FatesPFTIndexSwapper.py @@ -0,0 +1,225 @@ +#!/usr/bin/env python + +# ======================================================================================= +# +# This python script will open an input FATES parameter file, and given a list of PFT +# indices supplied by the user, will create a new parameter file with PFTs entries cloned +# from the original file as-per the list of indices supplied by the user. +# +# First Added, Ryan Knox: Thu Jan 11 13:36:14 PST 2018 +# ======================================================================================= + +import numpy as np +import sys +import getopt +import code # For development: code.interact(local=locals()) +from datetime import datetime +from scipy.io import netcdf +import matplotlib.pyplot as plt + +# ======================================================================================= +# Parameters +# ======================================================================================= + +pft_dim_name = 'fates_pft' + + +class timetype: + + # This is time, like the thing that always goes forward and cant be seen + # or touched, insert creative riddle here + + def __init__(self,ntimes): + + self.year = -9*np.ones((ntimes)) + self.month = -9*np.ones((ntimes)) + # This is a floating point decimal day + self.day = -9.0*np.ones((ntimes)) + + # This is a decimal datenumber + self.datenum = -9.0*np.ones((ntimes)) + + +def usage(): + print('') + print('=======================================================================') + print('') + print(' python FatesPFTIndexSwapper.py -h --pft-indices ') + print(' --fin ') + print(' --fout ') + print('') + print('') + print(' -h --help ') + print(' print this help message') + print('') + print('') + print(' --pft-indices ') + print(' This is a comma delimited list of integer positions of the PFTs') + print(' to be copied into the new file. Note that first pft position') + print(' is treated as 1 (not C or python like), and any order or multiples') + print(' of indices can be chosen') + print('') + print('') + print(' --fin ') + print(' This is the full path to the netcdf file you are basing off of') + print('') + print('') + print(' --fout ') + print(' This is the full path to the netcdf file you are writing to.') + print('') + print('') + print('=======================================================================') + + +def interp_args(argv): + + argv.pop(0) # The script itself is the first argument, forget it + + # Name of the conversion file + + input_fname = "none" + output_fname = "none" + donor_pft_indices = -9 + donot_pft_indices_str = '' + try: + opts, args = getopt.getopt(argv, 'h',["fin=","fout=","pft-indices="]) + except getopt.GetoptError as err: + print('Argument error, see usage') + usage() + sys.exit(2) + + for o, a in opts: + if o in ("-h", "--help"): + usage() + sys.exit(0) + elif o in ("--fin"): + input_fname = a + elif o in ("--fout"): + output_fname = a + elif o in ("--pft-indices"): + donor_pft_indices_str = a.strip() + else: + assert False, "unhandled option" + + + if (input_fname == "none"): + print("You must specify an input file:\n\n") + usage() + sys.exit(2) + + if (output_fname == "none"): + print("You must specify an output file:\n\n") + usage() + sys.exit(2) + + if (donor_pft_indices_str == ''): + print("You must specify at least one donor pft index!\n\n") + usage() + sys.exit(2) + else: + donor_pft_indices = [] + for strpft in donor_pft_indices_str.split(','): + donor_pft_indices.append(int(strpft)) + + + return (input_fname,output_fname,donor_pft_indices) + + +# ======================================================================================== +# ======================================================================================== +# Main +# ======================================================================================== +# ======================================================================================== + +def main(argv): + + # Interpret the arguments to the script + [input_fname,output_fname,donor_pft_indices] = interp_args(argv) + + num_pft_out = len(donor_pft_indices) + + # Open the netcdf files + fp_out = netcdf.netcdf_file(output_fname, 'w') + + fp_in = netcdf.netcdf_file(input_fname, 'r') + + for key, value in sorted(fp_in.dimensions.iteritems()): + if(key==pft_dim_name): + fp_out.createDimension(key,int(num_pft_out)) + print('Creating Dimension: {}={}'.format(key,num_pft_out)) + else: + fp_out.createDimension(key,int(value)) + print('Creating Dimension: {}={}'.format(key,value)) + + for key, value in sorted(fp_in.variables.iteritems()): + print('Creating Variable: ',key) + # code.interact(local=locals()) + + out_var = fp_out.createVariable(key,'f',(fp_in.variables.get(key).dimensions)) + in_var = fp_in.variables.get(key) + out_var.units = in_var.units + out_var.long_name = in_var.long_name + + # Idenfity if this variable has pft dimension + pft_dim_found = -1 + pft_dim_len = len(fp_in.variables.get(key).dimensions) + + for idim, name in enumerate(fp_in.variables.get(key).dimensions): + # Manipulate data + if(name==pft_dim_name): + pft_dim_found = idim + + # Copy over the input data + # Tedious, but I have to permute through all combinations of dimension position + if( pft_dim_len == 0 ): + out_var.assignValue(float(fp_in.variables.get(key).data)) + elif(pft_dim_found==-1): + out_var[:] = in_var[:] + elif( (pft_dim_found==0) & (pft_dim_len==1) ): # 1D fates_pft + tmp_out = np.zeros([num_pft_out]) + for id,ipft in enumerate(donor_pft_indices): + tmp_out[id] = fp_in.variables.get(key).data[ipft-1] + out_var[:] = tmp_out + + + elif( (pft_dim_found==1) & (pft_dim_len==2) ): # 2D hdyro_organ - fate_pft + dim2_len = fp_in.dimensions.get(fp_in.variables.get(key).dimensions[0]) + tmp_out = np.zeros([dim2_len,num_pft_out]) + for id,ipft in enumerate(donor_pft_indices): + for idim in range(0,dim2_len): + tmp_out[idim,id] = fp_in.variables.get(key).data[idim,ipft-1] + out_var[:] = tmp_out + else: + print('This variable has a dimensioning that we have not considered yet.') + print('Please add this condition to the logic above this statement.') + print('Aborting') + exit(2) + + fp_out.history = "This file was made from FatesPFTIndexSwapper.py \n Input File = {} \n Indices = {}"\ + .format(input_fname,donor_pft_indices) + + #var_out.mode = var.mode + #fp.flush() + + fp_in.close() + fp_out.close() + + print('Cloneing complete!') + exit(0) + + + + +# ======================================================================================= +# This is the actual call to main + +if __name__ == "__main__": + main(sys.argv) + + + + + + + + diff --git a/tools/modify_fates_paramfile.py b/tools/modify_fates_paramfile.py new file mode 100755 index 0000000000..92ab14dab1 --- /dev/null +++ b/tools/modify_fates_paramfile.py @@ -0,0 +1,102 @@ +#!/usr/bin/env python + +#### this script modifies a FATES parameter file. It accepts the following flags +# --var or --variable: variable. +# --pft or --PFT: PFT number. If this is missing, script will assume that its a global variable that is being modified. +# --input or --fin: input filename. +# --output or --fout: output filename. If missing, will assume its directly modifying the input file, and will prompt unless -O is specified +# --O or --overwrite: overwrite output file without asking. +# --value or --val: value to put in variable +# --s or --silent: don't write anything on successful execution. +#### + + +# ======================================================================================= +# ======================================================================================= + +import numpy as np +import os +from scipy.io import netcdf as nc +import argparse +import shutil + +# ======================================================================================== +# ======================================================================================== +# Main +# ======================================================================================== +# ======================================================================================== + +def main(): + parser = argparse.ArgumentParser(description='Parse command line arguments to this script.') + # + parser.add_argument('--var','--variable', dest='varname', type=str, help="What variable to modify? Required.", required=True) + parser.add_argument('--pft','--PFT', dest='pftnum', type=int, help="PFT number to modify. If this is missing, will assume a global variable.") + parser.add_argument('--fin', '--input', dest='inputfname', type=str, help="Input filename. Required.", required=True) + parser.add_argument('--fout','--output', dest='outputfname', type=str, help="Output filename. Required.", required=True) + parser.add_argument('--val', '--value', dest='val', type=float, help="New value of PFT variable. Required.", required=True) + parser.add_argument('--O','--overwrite', dest='overwrite', help="If present, automatically overwrite the output file.", action="store_true") + parser.add_argument('--silent', '--s', dest='silent', help="prevent writing of output.", action="store_true") + # + args = parser.parse_args() + # print(args.varname, args.pftnum, args.inputfname, args.outputfname, args.val, args.overwrite) + # + # check to see if output file exists + if os.path.isfile(args.outputfname): + if args.overwrite: + if not args.silent: + print('replacing file: '+args.outputfname) + os.remove(args.outputfname) + else: + raise ValueError('Output file already exists and overwrite flag not specified for filename: '+args.outputfname) + # + shutil.copyfile(args.inputfname, args.outputfname) + # + ncfile = nc.netcdf_file(args.outputfname, 'a') + # + var = ncfile.variables[args.varname] + # + ### check to make sure that, if a PFT is specified, the variable has a PFT dimension, and if not, then it doesn't. and also that shape is reasonable. + ndim_file = len(var.dimensions) + ispftvar = False + for i in range(ndim_file): + if var.dimensions[i] == 'fates_pft': + ispftvar = True + npft_file = var.shape[i] + pftdim = 0 + else: + npft_file = None + pftdim = None + if args.pftnum == None and ispftvar: + raise ValueError('pft value is missing but variable has pft dimension.') + if args.pftnum != None and not ispftvar: + raise ValueError('pft value is present but variable does not have pft dimension.') + if ndim_file > 1: + raise ValueError('variable dimensionality is too high for this script') + if ndim_file == 1 and not ispftvar: + raise ValueError('variable dimensionality is too high for this script') + if args.pftnum != None and ispftvar: + if args.pftnum > npft_file: + raise ValueError('PFT specified ('+str(args.pftnum)+') is larger than the number of PFTs in the file ('+str(npft_file)+').') + if pftdim == 0: + if not args.silent: + print('replacing prior value of variable '+args.varname+', for PFT '+str(args.pftnum)+', which was '+str(var[args.pftnum-1])+', with new value of '+str(args.val)) + var[args.pftnum-1] = args.val + elif args.pftnum == None and not ispftvar: + if not args.silent: + print('replacing prior value of variable '+args.varname+', which was '+str(var[:])+', with new value of '+str(args.val)) + var[:] = args.val + else: + raise ValueError('Nothing happened somehow.') + # + # + ncfile.close() + + + + +# ======================================================================================= +# This is the actual call to main + +if __name__ == "__main__": + main() +