From ee4344f118a1ed733c42b4f0962a3c722abb6940 Mon Sep 17 00:00:00 2001 From: Charlie Koven Date: Fri, 2 Nov 2018 17:46:37 -0700 Subject: [PATCH 01/34] initial attempt to allow equal promotion/demotion across equal-statured cohorts --- biogeochem/EDCanopyStructureMod.F90 | 178 ++++++++++++++++++++++++++-- 1 file changed, 171 insertions(+), 7 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index caf6861577..180b6a7452 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -345,6 +345,11 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) real(r8) :: arealayer ! the area of the current canopy layer integer :: exceedance_counter ! when seeking to rebalance demotion exceedance ! keep a loop counter to check for hangs + logical :: tied_size_with_neighbor, has_taller_equalsized_neighbor + logical :: found_shortest_equal_neighbor, found_tallest_equal_neighbor + type(ed_cohort_type), pointer :: cohort_tosearch_relative_to, cohort_tocompare_to + real(r8) :: total_crownarea_of_tied_cohorts + real(r8) :: sumweights_equalsizebuffer ! First, determine how much total canopy area we have in this layer @@ -359,7 +364,8 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) ! In that case, we need to work out which cohorts to demote. ! We go in order from shortest to tallest for ranked demotion - sumweights = 0.0_r8 + sumweights = 0.0_r8 + sumweights_equalsizebuffer = 0.0_r8 currentCohort => currentPatch%shortest do while (associated(currentCohort)) @@ -374,10 +380,85 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) currentCohort%n/(currentCohort%dbh**ED_val_comp_excln) else ! Rank ordered deterministic method - currentCohort%excl_weight = & - max(min(currentCohort%c_area, demote_area - sumweights ), 0._r8) + + ! check to make sure there are no cohorts of equal size + tied_size_with_neighbor = .false. + if (associated(currentCohort%shorter)) then + if (currentCohort%shorter%dbh .eq. currentCohort%dbh ) then + tied_size_with_neighbor = .true. + endif + endif + if (associated(currentCohort%taller)) then + if (currentCohort%taller%dbh .eq. currentCohort%dbh ) then + tied_size_with_neighbor = .true. + endif + endif + + if ( tied_size_with_neighbor ) then + ! now we need to go through and figure out how many equal-size cohorts there are. + ! then we need to go through, add up the collective corwn areas of all equal-sized and equal-canopy-layer cohorts, + ! and then demote from each as if they were a single group + ! + total_crownarea_of_tied_cohorts = currentCohort%c_area + ! + ! first the "shorter" cohorts (scare-quotes because they aren't actually shorter) + found_shortest_equal_neighbor = .false. + cohort_tosearch_relative_to = currentCohort + do while ( .not. found_shortest_equal_neighbor) + if (associated(cohort_tosearch_relative_to%shorter)) then + cohort_tocompare_to = cohort_tosearch_relative_to%shorter + if (cohort_tocompare_to%dbh .eq. currentCohort%dbh ) then + if (cohort_tocompare_to%canopy_layer .eq. currentCohort%canopy_layer ) then + total_crownarea_of_tied_cohorts = total_crownarea_of_tied_cohorts + cohort_tocompare_to%c_area + endif + cohort_tosearch_relative_to = cohort_tocompare_to + else + found_shortest_equal_neighbor = .true. + end if + else + found_shortest_equal_neighbor = .true. + endif + end do + ! + ! then the "taller" cohorts (scare-quotes because they aren't actually taller) + has_taller_equalsized_neighbor = .false. ! init this as false + found_tallest_equal_neighbor = .false. + cohort_tosearch_relative_to = currentCohort + do while ( .not. found_tallest_equal_neighbor) + if (associated(cohort_tosearch_relative_to%taller)) then + cohort_tocompare_to = cohort_tosearch_relative_to%taller + if (cohort_tocompare_to%dbh .eq. currentCohort%dbh ) then + if (cohort_tocompare_to%canopy_layer .eq. currentCohort%canopy_layer ) then + total_crownarea_of_tied_cohorts = total_crownarea_of_tied_cohorts + cohort_tocompare_to%c_area + has_taller_equalsized_neighbor = .true. + endif + cohort_tosearch_relative_to = cohort_tocompare_to + else + found_tallest_equal_neighbor = .true. + end if + else + found_tallest_equal_neighbor = .true. + endif + end do + ! + ! now we know the total crown area of all equal-sized, equal-canopy-layer cohorts + currentCohort%excl_weight = & + max(min(currentCohort%c_area, (currentCohort%c_area/total_crownarea_of_tied_cohorts) * & + (demote_area - sumweights) ), 0._r8) + else + currentCohort%excl_weight = & + max(min(currentCohort%c_area, demote_area - sumweights ), 0._r8) + endif + endif + if ((ED_val_comp_excln .lt. 0.0_r8) .and. tied_size_with_neighbor .and. & + has_taller_equalsized_neighbor) then + sumweights_equalsizebuffer = sumweights_equalsizebuffer + currentCohort%excl_weight + else if (ED_val_comp_excln .lt. 0.0_r8) .and. tied_size_with_neighbor) then + sumweights = sumweights + sumweights_equalsizebuffer + sumweights_equalsizebuffer = 0._r8 + else + sumweights = sumweights + currentCohort%excl_weight endif - sumweights = sumweights + currentCohort%excl_weight endif currentCohort => currentCohort%taller enddo @@ -669,6 +750,13 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) real(r8) :: sapw_c ! sapwood carbon [kg] real(r8) :: store_c ! storage carbon [kg] real(r8) :: struct_c ! structure carbon [kg] + + logical :: tied_size_with_neighbor, has_shorter_equalsized_neighbor + logical :: found_shortest_equal_neighbor, found_tallest_equal_neighbor + type(ed_cohort_type), pointer :: cohort_tosearch_relative_to, cohort_tocompare_to + real(r8) :: total_crownarea_of_tied_cohorts + real(r8) :: sumweights_equalsizebuffer + call CanopyLayerArea(currentPatch,currentSite%spread,i_lyr,arealayer_current) call CanopyLayerArea(currentPatch,currentSite%spread,i_lyr+1,arealayer_below) @@ -722,6 +810,7 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) ! This is the opposite of the demotion weighting... sumweights = 0.0_r8 + sumweights_equalsizebuffer = 0.0_r8 currentCohort => currentPatch%tallest do while (associated(currentCohort)) call carea_allom(currentCohort%dbh,currentCohort%n,currentSite%spread, & @@ -731,10 +820,85 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) ! normal (stochastic) case, as above. currentCohort%prom_weight = currentCohort%n*currentCohort%dbh**ED_val_comp_excln else - currentCohort%prom_weight = max(min(currentCohort%c_area, & - promote_area - sumweights ), 0._r8) + ! Rank ordered deterministic method + + ! check to make sure there are no cohorts of equal size + tied_size_with_neighbor = .false. + if (associated(currentCohort%shorter)) then + if (currentCohort%shorter%dbh .eq. currentCohort%dbh ) then + tied_size_with_neighbor = .true. + endif + endif + if (associated(currentCohort%taller)) then + if (currentCohort%taller%dbh .eq. currentCohort%dbh ) then + tied_size_with_neighbor = .true. + endif + endif + + if ( tied_size_with_neighbor ) then + ! now we need to go through and figure out how many equal-size cohorts there are. + ! then we need to go through, add up the collective corwn areas of all equal-sized and equal-canopy-layer cohorts, + ! and then demote from each as if they were a single group + ! + total_crownarea_of_tied_cohorts = currentCohort%c_area + ! + ! first the "shorter" cohorts (scare-quotes because they aren't actually shorter) + found_shortest_equal_neighbor = .false. + cohort_tosearch_relative_to = currentCohort + do while ( .not. found_shortest_equal_neighbor) + if (associated(cohort_tosearch_relative_to%shorter)) then + cohort_tocompare_to = cohort_tosearch_relative_to%shorter + if (cohort_tocompare_to%dbh .eq. currentCohort%dbh ) then + if (cohort_tocompare_to%canopy_layer .eq. currentCohort%canopy_layer ) then + total_crownarea_of_tied_cohorts = total_crownarea_of_tied_cohorts + cohort_tocompare_to%c_area + has_shorter_equalsized_neighbor = .true. + endif + cohort_tosearch_relative_to = cohort_tocompare_to + else + found_shortest_equal_neighbor = .true. + end if + else + found_shortest_equal_neighbor = .true. + endif + end do + ! + ! then the "taller" cohorts (scare-quotes because they aren't actually taller) + has_shorter_equalsized_neighbor = .false. ! init this as false + found_tallest_equal_neighbor = .false. + cohort_tosearch_relative_to = currentCohort + do while ( .not. found_tallest_equal_neighbor) + if (associated(cohort_tosearch_relative_to%taller)) then + cohort_tocompare_to = cohort_tosearch_relative_to%taller + if (cohort_tocompare_to%dbh .eq. currentCohort%dbh ) then + if (cohort_tocompare_to%canopy_layer .eq. currentCohort%canopy_layer ) then + total_crownarea_of_tied_cohorts = total_crownarea_of_tied_cohorts + cohort_tocompare_to%c_area + endif + cohort_tosearch_relative_to = cohort_tocompare_to + else + found_tallest_equal_neighbor = .true. + end if + else + found_tallest_equal_neighbor = .true. + endif + end do + ! + ! now we know the total crown area of all equal-sized, equal-canopy-layer cohorts + currentCohort%prom_weight = max(min(currentCohort%c_area, & + (currentCohort%c_area/total_crownarea_of_tied_cohorts) * (promote_area - sumweights) ), 0._r8) + else + currentCohort%prom_weight = max(min(currentCohort%c_area, & + promote_area - sumweights ), 0._r8) + endif + endif + if ((ED_val_comp_excln .lt. 0.0_r8) .and. tied_size_with_neighbor .and. & + has_shorter_equalsized_neighbor) then + sumweights_equalsizebuffer = sumweights_equalsizebuffer + currentCohort%prom_weight + else if (ED_val_comp_excln .lt. 0.0_r8) .and. tied_size_with_neighbor) then + sumweights = sumweights + sumweights_equalsizebuffer + sumweights_equalsizebuffer = 0._r8 + else + sumweights = sumweights + currentCohort%prom_weight endif - sumweights = sumweights + currentCohort%prom_weight endif currentCohort => currentCohort%shorter enddo !currentCohort From d91fd6e61de1f33da066ea3c5fddd92391905d75 Mon Sep 17 00:00:00 2001 From: Charlie Koven Date: Fri, 2 Nov 2018 17:53:38 -0700 Subject: [PATCH 02/34] some bug and typo fixes on prior --- biogeochem/EDCanopyStructureMod.F90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 180b6a7452..65efb87a71 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -396,7 +396,7 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) if ( tied_size_with_neighbor ) then ! now we need to go through and figure out how many equal-size cohorts there are. - ! then we need to go through, add up the collective corwn areas of all equal-sized and equal-canopy-layer cohorts, + ! then we need to go through, add up the collective crown areas of all equal-sized and equal-canopy-layer cohorts, ! and then demote from each as if they were a single group ! total_crownarea_of_tied_cohorts = currentCohort%c_area @@ -454,7 +454,7 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) has_taller_equalsized_neighbor) then sumweights_equalsizebuffer = sumweights_equalsizebuffer + currentCohort%excl_weight else if (ED_val_comp_excln .lt. 0.0_r8) .and. tied_size_with_neighbor) then - sumweights = sumweights + sumweights_equalsizebuffer + sumweights = sumweights + currentCohort%excl_weight + sumweights_equalsizebuffer sumweights_equalsizebuffer = 0._r8 else sumweights = sumweights + currentCohort%excl_weight @@ -837,8 +837,8 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) if ( tied_size_with_neighbor ) then ! now we need to go through and figure out how many equal-size cohorts there are. - ! then we need to go through, add up the collective corwn areas of all equal-sized and equal-canopy-layer cohorts, - ! and then demote from each as if they were a single group + ! then we need to go through, add up the collective crown areas of all equal-sized and equal-canopy-layer cohorts, + ! and then promote from each as if they were a single group ! total_crownarea_of_tied_cohorts = currentCohort%c_area ! @@ -894,7 +894,7 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) has_shorter_equalsized_neighbor) then sumweights_equalsizebuffer = sumweights_equalsizebuffer + currentCohort%prom_weight else if (ED_val_comp_excln .lt. 0.0_r8) .and. tied_size_with_neighbor) then - sumweights = sumweights + sumweights_equalsizebuffer + sumweights = sumweights + currentCohort%prom_weight + sumweights_equalsizebuffer sumweights_equalsizebuffer = 0._r8 else sumweights = sumweights + currentCohort%prom_weight From 53bc103a3163a59e3bc3be021965b259d7860510 Mon Sep 17 00:00:00 2001 From: Charles Koven Date: Sat, 3 Nov 2018 19:05:39 -0700 Subject: [PATCH 03/34] fixed some compile-time and runtime bugs --- biogeochem/EDCanopyStructureMod.F90 | 65 ++++++++++++++++++++++------- 1 file changed, 49 insertions(+), 16 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 65efb87a71..b9d978543c 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -26,6 +26,7 @@ module EDCanopyStructureMod use FatesInterfaceMod , only : hlm_use_planthydro use FatesInterfaceMod , only : numpft use FatesPlantHydraulicsMod, only : UpdateH2OVeg,InitHydrCohort + use EDTypesMod , only : maxCohortsPerPatch use PRTGenericMod, only : leaf_organ use PRTGenericMod, only : all_carbon_elements @@ -350,6 +351,7 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) type(ed_cohort_type), pointer :: cohort_tosearch_relative_to, cohort_tocompare_to real(r8) :: total_crownarea_of_tied_cohorts real(r8) :: sumweights_equalsizebuffer + integer :: whileloop_counter ! First, determine how much total canopy area we have in this layer @@ -395,6 +397,9 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) endif if ( tied_size_with_neighbor ) then + if ( DEBUG ) then + write(fates_log(),*) 'tied_size_with_neighbor eq true in demotion phase' + endif ! now we need to go through and figure out how many equal-size cohorts there are. ! then we need to go through, add up the collective crown areas of all equal-sized and equal-canopy-layer cohorts, ! and then demote from each as if they were a single group @@ -403,42 +408,54 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) ! ! first the "shorter" cohorts (scare-quotes because they aren't actually shorter) found_shortest_equal_neighbor = .false. - cohort_tosearch_relative_to = currentCohort + cohort_tosearch_relative_to => currentCohort + whileloop_counter = 0 do while ( .not. found_shortest_equal_neighbor) + whileloop_counter = whileloop_counter + 1 if (associated(cohort_tosearch_relative_to%shorter)) then - cohort_tocompare_to = cohort_tosearch_relative_to%shorter + cohort_tocompare_to => cohort_tosearch_relative_to%shorter if (cohort_tocompare_to%dbh .eq. currentCohort%dbh ) then if (cohort_tocompare_to%canopy_layer .eq. currentCohort%canopy_layer ) then total_crownarea_of_tied_cohorts = total_crownarea_of_tied_cohorts + cohort_tocompare_to%c_area endif - cohort_tosearch_relative_to = cohort_tocompare_to + cohort_tosearch_relative_to => cohort_tocompare_to else found_shortest_equal_neighbor = .true. end if else found_shortest_equal_neighbor = .true. endif + if ( whileloop_counter .ge. maxCohortsPerPatch ) then + ! something has gone horribly wrong and we are in an infite loop. + call endrun(msg=errMsg(sourcefile, __LINE__)) + endif end do ! ! then the "taller" cohorts (scare-quotes because they aren't actually taller) has_taller_equalsized_neighbor = .false. ! init this as false found_tallest_equal_neighbor = .false. - cohort_tosearch_relative_to = currentCohort + cohort_tosearch_relative_to => currentCohort + whileloop_counter = 0 do while ( .not. found_tallest_equal_neighbor) + whileloop_counter = whileloop_counter + 1 if (associated(cohort_tosearch_relative_to%taller)) then - cohort_tocompare_to = cohort_tosearch_relative_to%taller + cohort_tocompare_to => cohort_tosearch_relative_to%taller if (cohort_tocompare_to%dbh .eq. currentCohort%dbh ) then if (cohort_tocompare_to%canopy_layer .eq. currentCohort%canopy_layer ) then total_crownarea_of_tied_cohorts = total_crownarea_of_tied_cohorts + cohort_tocompare_to%c_area has_taller_equalsized_neighbor = .true. endif - cohort_tosearch_relative_to = cohort_tocompare_to + cohort_tosearch_relative_to => cohort_tocompare_to else found_tallest_equal_neighbor = .true. end if else found_tallest_equal_neighbor = .true. endif + if ( whileloop_counter .ge. maxCohortsPerPatch ) then + ! something has gone horribly wrong and we are in an infite loop. + call endrun(msg=errMsg(sourcefile, __LINE__)) + endif end do ! ! now we know the total crown area of all equal-sized, equal-canopy-layer cohorts @@ -453,7 +470,7 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) if ((ED_val_comp_excln .lt. 0.0_r8) .and. tied_size_with_neighbor .and. & has_taller_equalsized_neighbor) then sumweights_equalsizebuffer = sumweights_equalsizebuffer + currentCohort%excl_weight - else if (ED_val_comp_excln .lt. 0.0_r8) .and. tied_size_with_neighbor) then + else if ( (ED_val_comp_excln .lt. 0.0_r8) .and. tied_size_with_neighbor) then sumweights = sumweights + currentCohort%excl_weight + sumweights_equalsizebuffer sumweights_equalsizebuffer = 0._r8 else @@ -756,7 +773,7 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) type(ed_cohort_type), pointer :: cohort_tosearch_relative_to, cohort_tocompare_to real(r8) :: total_crownarea_of_tied_cohorts real(r8) :: sumweights_equalsizebuffer - + integer :: whileloop_counter call CanopyLayerArea(currentPatch,currentSite%spread,i_lyr,arealayer_current) call CanopyLayerArea(currentPatch,currentSite%spread,i_lyr+1,arealayer_below) @@ -836,6 +853,10 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) endif if ( tied_size_with_neighbor ) then + if ( DEBUG ) then + write(fates_log(),*) 'tied_size_with_neighbor eq true in promotion phase' + endif + ! now we need to go through and figure out how many equal-size cohorts there are. ! then we need to go through, add up the collective crown areas of all equal-sized and equal-canopy-layer cohorts, ! and then promote from each as if they were a single group @@ -843,43 +864,55 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) total_crownarea_of_tied_cohorts = currentCohort%c_area ! ! first the "shorter" cohorts (scare-quotes because they aren't actually shorter) + has_shorter_equalsized_neighbor = .false. ! init this as false found_shortest_equal_neighbor = .false. - cohort_tosearch_relative_to = currentCohort + cohort_tosearch_relative_to => currentCohort + whileloop_counter = 0 do while ( .not. found_shortest_equal_neighbor) + whileloop_counter = whileloop_counter + 1 if (associated(cohort_tosearch_relative_to%shorter)) then - cohort_tocompare_to = cohort_tosearch_relative_to%shorter + cohort_tocompare_to => cohort_tosearch_relative_to%shorter if (cohort_tocompare_to%dbh .eq. currentCohort%dbh ) then if (cohort_tocompare_to%canopy_layer .eq. currentCohort%canopy_layer ) then total_crownarea_of_tied_cohorts = total_crownarea_of_tied_cohorts + cohort_tocompare_to%c_area has_shorter_equalsized_neighbor = .true. endif - cohort_tosearch_relative_to = cohort_tocompare_to + cohort_tosearch_relative_to => cohort_tocompare_to else found_shortest_equal_neighbor = .true. end if else found_shortest_equal_neighbor = .true. endif + if ( whileloop_counter .ge. maxCohortsPerPatch ) then + ! something has gone horribly wrong and we are in an infite loop. + call endrun(msg=errMsg(sourcefile, __LINE__)) + endif end do ! ! then the "taller" cohorts (scare-quotes because they aren't actually taller) - has_shorter_equalsized_neighbor = .false. ! init this as false found_tallest_equal_neighbor = .false. - cohort_tosearch_relative_to = currentCohort + cohort_tosearch_relative_to => currentCohort + whileloop_counter = 0 do while ( .not. found_tallest_equal_neighbor) + whileloop_counter = whileloop_counter + 1 if (associated(cohort_tosearch_relative_to%taller)) then - cohort_tocompare_to = cohort_tosearch_relative_to%taller + cohort_tocompare_to => cohort_tosearch_relative_to%taller if (cohort_tocompare_to%dbh .eq. currentCohort%dbh ) then if (cohort_tocompare_to%canopy_layer .eq. currentCohort%canopy_layer ) then total_crownarea_of_tied_cohorts = total_crownarea_of_tied_cohorts + cohort_tocompare_to%c_area endif - cohort_tosearch_relative_to = cohort_tocompare_to + cohort_tosearch_relative_to => cohort_tocompare_to else found_tallest_equal_neighbor = .true. end if else found_tallest_equal_neighbor = .true. endif + if ( whileloop_counter .ge. maxCohortsPerPatch ) then + ! something has gone horribly wrong and we are in an infite loop. + call endrun(msg=errMsg(sourcefile, __LINE__)) + endif end do ! ! now we know the total crown area of all equal-sized, equal-canopy-layer cohorts @@ -893,7 +926,7 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) if ((ED_val_comp_excln .lt. 0.0_r8) .and. tied_size_with_neighbor .and. & has_shorter_equalsized_neighbor) then sumweights_equalsizebuffer = sumweights_equalsizebuffer + currentCohort%prom_weight - else if (ED_val_comp_excln .lt. 0.0_r8) .and. tied_size_with_neighbor) then + else if ( (ED_val_comp_excln .lt. 0.0_r8) .and. tied_size_with_neighbor) then sumweights = sumweights + currentCohort%prom_weight + sumweights_equalsizebuffer sumweights_equalsizebuffer = 0._r8 else From 5a639751a41856e64ae4d96ddaffb55de54a8b9d Mon Sep 17 00:00:00 2001 From: Charlie Koven Date: Wed, 7 Nov 2018 16:36:41 -0500 Subject: [PATCH 04/34] swapped dbh for height in cohort sorting and canopy structure logic --- biogeochem/EDCanopyStructureMod.F90 | 20 ++++++++++---------- biogeochem/EDCohortDynamicsMod.F90 | 4 ++-- 2 files changed, 12 insertions(+), 12 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index b9d978543c..a703e7d725 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -379,19 +379,19 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) ! normal (stochastic) case. weight cohort demotion by ! inverse size to a constant power currentCohort%excl_weight = & - currentCohort%n/(currentCohort%dbh**ED_val_comp_excln) + currentCohort%n/(currentCohort%hite**ED_val_comp_excln) else ! Rank ordered deterministic method ! check to make sure there are no cohorts of equal size tied_size_with_neighbor = .false. if (associated(currentCohort%shorter)) then - if (currentCohort%shorter%dbh .eq. currentCohort%dbh ) then + if (currentCohort%shorter%hite .eq. currentCohort%hite ) then tied_size_with_neighbor = .true. endif endif if (associated(currentCohort%taller)) then - if (currentCohort%taller%dbh .eq. currentCohort%dbh ) then + if (currentCohort%taller%hite .eq. currentCohort%hite ) then tied_size_with_neighbor = .true. endif endif @@ -414,7 +414,7 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) whileloop_counter = whileloop_counter + 1 if (associated(cohort_tosearch_relative_to%shorter)) then cohort_tocompare_to => cohort_tosearch_relative_to%shorter - if (cohort_tocompare_to%dbh .eq. currentCohort%dbh ) then + if (cohort_tocompare_to%hite .eq. currentCohort%hite ) then if (cohort_tocompare_to%canopy_layer .eq. currentCohort%canopy_layer ) then total_crownarea_of_tied_cohorts = total_crownarea_of_tied_cohorts + cohort_tocompare_to%c_area endif @@ -440,7 +440,7 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) whileloop_counter = whileloop_counter + 1 if (associated(cohort_tosearch_relative_to%taller)) then cohort_tocompare_to => cohort_tosearch_relative_to%taller - if (cohort_tocompare_to%dbh .eq. currentCohort%dbh ) then + if (cohort_tocompare_to%hite .eq. currentCohort%hite ) then if (cohort_tocompare_to%canopy_layer .eq. currentCohort%canopy_layer ) then total_crownarea_of_tied_cohorts = total_crownarea_of_tied_cohorts + cohort_tocompare_to%c_area has_taller_equalsized_neighbor = .true. @@ -835,19 +835,19 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) 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. - currentCohort%prom_weight = currentCohort%n*currentCohort%dbh**ED_val_comp_excln + currentCohort%prom_weight = currentCohort%n*currentCohort%hite**ED_val_comp_excln else ! Rank ordered deterministic method ! check to make sure there are no cohorts of equal size tied_size_with_neighbor = .false. if (associated(currentCohort%shorter)) then - if (currentCohort%shorter%dbh .eq. currentCohort%dbh ) then + if (currentCohort%shorter%hite .eq. currentCohort%hite ) then tied_size_with_neighbor = .true. endif endif if (associated(currentCohort%taller)) then - if (currentCohort%taller%dbh .eq. currentCohort%dbh ) then + if (currentCohort%taller%hite .eq. currentCohort%hite ) then tied_size_with_neighbor = .true. endif endif @@ -872,7 +872,7 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) whileloop_counter = whileloop_counter + 1 if (associated(cohort_tosearch_relative_to%shorter)) then cohort_tocompare_to => cohort_tosearch_relative_to%shorter - if (cohort_tocompare_to%dbh .eq. currentCohort%dbh ) then + if (cohort_tocompare_to%hite .eq. currentCohort%hite ) then if (cohort_tocompare_to%canopy_layer .eq. currentCohort%canopy_layer ) then total_crownarea_of_tied_cohorts = total_crownarea_of_tied_cohorts + cohort_tocompare_to%c_area has_shorter_equalsized_neighbor = .true. @@ -898,7 +898,7 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) whileloop_counter = whileloop_counter + 1 if (associated(cohort_tosearch_relative_to%taller)) then cohort_tocompare_to => cohort_tosearch_relative_to%taller - if (cohort_tocompare_to%dbh .eq. currentCohort%dbh ) then + if (cohort_tocompare_to%hite .eq. currentCohort%hite ) then if (cohort_tocompare_to%canopy_layer .eq. currentCohort%canopy_layer ) then total_crownarea_of_tied_cohorts = total_crownarea_of_tied_cohorts + cohort_tocompare_to%c_area endif diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 71a1416cfc..0e014b4b5d 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -1146,7 +1146,7 @@ subroutine insert_cohort(pcc, ptall, pshort, tnull, snull, storebigcohort, store icohort => pcc ! assign address to icohort local name !place in the correct place in the linked list of heights !begin by finding cohort that is just taller than the new cohort - tsp = icohort%dbh + tsp = icohort%hite current => pshortest exitloop = 0 @@ -1154,7 +1154,7 @@ subroutine insert_cohort(pcc, ptall, pshort, tnull, snull, storebigcohort, store !taller than tree being considered and return its pointer if (associated(current)) then do while (associated(current).and.exitloop == 0) - if (current%dbh < tsp) then + if (current%hite < tsp) then current => current%taller else exitloop = 1 From 143ab36ade8d60263411f913f5a41bbe7ca3fc2c Mon Sep 17 00:00:00 2001 From: Charlie Koven Date: Wed, 7 Nov 2018 16:44:52 -0500 Subject: [PATCH 05/34] added some comments --- biogeochem/EDCanopyStructureMod.F90 | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index a703e7d725..e30d339bb7 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -462,11 +462,14 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) currentCohort%excl_weight = & max(min(currentCohort%c_area, (currentCohort%c_area/total_crownarea_of_tied_cohorts) * & (demote_area - sumweights) ), 0._r8) - else + else ! i.e. tied_size_with_neighbor = .false. currentCohort%excl_weight = & max(min(currentCohort%c_area, demote_area - sumweights ), 0._r8) endif endif + ! + ! when two or more cohorts have the same size, we need to keep track of their cumulative demoted crown area + ! in a separate buffer and add it once we reach the last of the equal-sized cohorts if ((ED_val_comp_excln .lt. 0.0_r8) .and. tied_size_with_neighbor .and. & has_taller_equalsized_neighbor) then sumweights_equalsizebuffer = sumweights_equalsizebuffer + currentCohort%excl_weight @@ -918,11 +921,14 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) ! now we know the total crown area of all equal-sized, equal-canopy-layer cohorts currentCohort%prom_weight = max(min(currentCohort%c_area, & (currentCohort%c_area/total_crownarea_of_tied_cohorts) * (promote_area - sumweights) ), 0._r8) - else + else ! i.e. tied_size_with_neighbor = .false. currentCohort%prom_weight = max(min(currentCohort%c_area, & promote_area - sumweights ), 0._r8) endif endif + ! + ! when two or more cohorts have the same size, we need to keep track of their cumulative demoted crown area + ! in a separate buffer and add it once we reach the last of the equal-sized cohorts if ((ED_val_comp_excln .lt. 0.0_r8) .and. tied_size_with_neighbor .and. & has_shorter_equalsized_neighbor) then sumweights_equalsizebuffer = sumweights_equalsizebuffer + currentCohort%prom_weight From 2613b5b3d7f2c61342dbc0e64864d8a5c361a783 Mon Sep 17 00:00:00 2001 From: chonggang xu Date: Sat, 15 Dec 2018 14:17:18 -0800 Subject: [PATCH 06/34] Adding types and parameter structure for leaf age discretization --- main/EDInitMod.F90 | 6 +- main/EDPftvarcon.F90 | 86 +++++++++++++++++++----- main/EDTypesMod.F90 | 17 +++++ main/FatesRestartInterfaceMod.F90 | 7 +- parameter_files/fates_params_14pfts.cdl | 5 +- parameter_files/fates_params_default.cdl | 15 +++-- 6 files changed, 112 insertions(+), 24 deletions(-) diff --git a/main/EDInitMod.F90 b/main/EDInitMod.F90 index 696ba70a72..516a1ac00f 100644 --- a/main/EDInitMod.F90 +++ b/main/EDInitMod.F90 @@ -21,10 +21,12 @@ module EDInitMod use EDTypesMod , only : AREA use EDTypesMod , only : init_spread_near_bare_ground use EDTypesMod , only : init_spread_inventory + use EDTypesMod , only : first_leaf_aclass use FatesInterfaceMod , only : bc_in_type use FatesInterfaceMod , only : hlm_use_planthydro use FatesInterfaceMod , only : hlm_use_inventory_init use FatesInterfaceMod , only : numpft + use FatesInterfaceMod , only : nleafage use ChecksBalancesMod , only : SiteCarbonStock use FatesInterfaceMod , only : nlevsclass use FatesAllometryMod , only : h2d_allom @@ -447,8 +449,8 @@ subroutine init_cohorts( site_in, patch_in, bc_in) call create_cohort(site_in, patch_in, pft, temp_cohort%n, temp_cohort%hite, temp_cohort%dbh, & b_leaf, b_fineroot, b_sapwood, b_dead, b_store, & - temp_cohort%laimemory, cstatus, rstatus, temp_cohort%canopy_trim, 1, site_in%spread, bc_in) - + temp_cohort%laimemory, cstatus, rstatus, temp_cohort%canopy_trim, 1, & + site_in%spread, first_leaf_aclass, bc_in) deallocate(temp_cohort) ! get rid of temporary cohort diff --git a/main/EDPftvarcon.F90 b/main/EDPftvarcon.F90 index b32532f86e..434a28bbf1 100644 --- a/main/EDPftvarcon.F90 +++ b/main/EDPftvarcon.F90 @@ -81,8 +81,7 @@ module EDPftvarcon ! decreases light interception real(r8), allocatable :: c3psn(:) ! index defining the photosynthetic ! pathway C4 = 0, C3 = 1 - real(r8), allocatable :: vcmax25top(:) ! maximum carboxylation rate of Rub. at 25C, - ! canopy top [umol CO2/m^2/s] + real(r8), allocatable :: smpso(:) ! Soil water potential at full stomatal opening ! (non-HYDRO mode only) [mm] real(r8), allocatable :: smpsc(:) ! Soil water potential at full stomatal closure @@ -202,7 +201,7 @@ module EDPftvarcon real(r8), allocatable :: phenflush_fraction(:) ! Maximum fraction of storage carbon used to flush leaves ! on bud-burst [kgC/kgC] - real(r8), allocatable :: leaf_long(:) ! Leaf turnover time (longevity) (pft) [yr] + real(r8), allocatable :: root_long(:) ! root turnover time (longevity) (pft) [yr] real(r8), allocatable :: branch_turnover(:) ! Turnover time for branchfall on live trees (pft) [yr] real(r8), allocatable :: turnover_retrans_mode(:) ! Retranslocation method (pft) @@ -211,7 +210,14 @@ module EDPftvarcon real(r8), allocatable :: turnover_nitr_retrans(:,:) ! nitrogen re-translocation fraction (pft x organ) real(r8), allocatable :: turnover_phos_retrans(:,:) ! phosphorous re-translocation fraction (pft x organ) + ! Parameters dimensioned by PFT and leaf age + real(r8), allocatable :: leaf_long(:,:) ! Leaf turnover time (longevity) (pft x age-class) + ! If there is >1 class, it is the longevity from + ! one class to the next [yr] + real(r8), allocatable :: vcmax25top(:,:) ! maximum carboxylation rate of Rub. at 25C, + ! canopy top [umol CO2/m^2/s]. Dimensioned by + ! leaf age-class ! Plant Hydraulic Parameters ! --------------------------------------------------------------------------------------------- @@ -246,6 +252,8 @@ module EDPftvarcon procedure, private :: Receive_PFT_hydr_organs procedure, private :: Register_PFT_prt_organs procedure, private :: Receive_PFT_prt_organs + procedure, private :: Register_PFT_leafage + procedure, private :: Receive_PFT_leafage procedure, private :: Register_PFT_numrad procedure, private :: Receive_PFT_numrad end type EDPftvarcon_type @@ -288,6 +296,7 @@ subroutine Register(this, fates_params) call this%Register_PFT_nvariants(fates_params) call this%Register_PFT_hydr_organs(fates_params) call this%Register_PFT_prt_organs(fates_params) + call this%Register_PFT_leafage(fates_params) end subroutine Register @@ -306,6 +315,7 @@ subroutine Receive(this, fates_params) call this%Receive_PFT_nvariants(fates_params) call this%Receive_PFT_hydr_organs(fates_params) call this%Receive_PFT_prt_organs(fates_params) + call this%Receive_PFT_leafage(fates_params) end subroutine Receive @@ -426,10 +436,6 @@ subroutine Register_PFT(this, fates_params) call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) - name = 'fates_leaf_long' - call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & - dimension_names=dim_names, lower_bounds=dim_lower_bound) - name = 'fates_roota_par' call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) @@ -475,10 +481,6 @@ subroutine Register_PFT(this, fates_params) call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) - name = 'fates_leaf_vcmax25top' - call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & - dimension_names=dim_names, lower_bounds=dim_lower_bound) - name = 'fates_smpso' call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) @@ -906,10 +908,6 @@ subroutine Receive_PFT(this, fates_params) call fates_params%RetreiveParameterAllocate(name=name, & data=this%c3psn) - name = 'fates_leaf_vcmax25top' - call fates_params%RetreiveParameterAllocate(name=name, & - data=this%vcmax25top) - name = 'fates_smpso' call fates_params%RetreiveParameterAllocate(name=name, & data=this%smpso) @@ -1403,6 +1401,64 @@ subroutine Receive_PFT_nvariants(this, fates_params) end subroutine Receive_PFT_nvariants ! ----------------------------------------------------------------------- + + subroutine Register_PFT_leafage(this, fates_params) + + use FatesParametersInterface, only : fates_parameters_type, param_string_length + use FatesParametersInterface, only : max_dimensions, dimension_name_leaf_age + use FatesParametersInterface, only : dimension_name_pft, dimension_shape_2d + + implicit none + + class(EDPftvarcon_type), intent(inout) :: this + class(fates_parameters_type), intent(inout) :: fates_params + + integer, parameter :: dim_lower_bound(2) = (/ lower_bound_pft, lower_bound_general /) + character(len=param_string_length) :: dim_names(2) + character(len=param_string_length) :: name + + dim_names(1) = dimension_name_pft + dim_names(2) = dimension_name_leaf_age + + name = 'fates_leaf_long' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_2d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) + + name = 'fates_leaf_vcmax25top' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_2d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) + + + return + end subroutine Register_PFT_leafage + + ! ===================================================================================== + + subroutine Receive_PFT_leafage(this, fates_params) + + use FatesParametersInterface, only : fates_parameters_type + use FatesParametersInterface, only : param_string_length + + implicit none + + class(EDPftvarcon_type), intent(inout) :: this + class(fates_parameters_type), intent(inout) :: fates_params + + character(len=param_string_length) :: name + + name = 'fates_leaf_long' + call fates_params%RetreiveParameterAllocate(name=name, & + data=this%leaf_long) + + name = 'fates_leaf_vcmax25top' + call fates_params%RetreiveParameterAllocate(name=name, & + data=this%vcmax25top) + + return + end subroutine Receive_PFT_leafage + + ! ===================================================================================== + subroutine Register_PFT_prt_organs(this, fates_params) diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 16bd7db647..497c872d8a 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -156,6 +156,16 @@ module EDTypesMod ! special mode to cause PFTs to create seed mass of all currently-existing PFTs logical, parameter :: homogenize_seed_pfts = .false. + + ! Leaf age class initialization schemes + integer, parameter :: nan_leaf_aclass = 0 ! initialize leaf age classes as undefined + ! (used when copying) + integer, parameter :: equal_leaf_aclass = 1 ! initialize leaf age classes equal + ! (used for inventory initialization) + integer, parameter :: first_leaf_aclass = 2 ! initialize leaf age classes as all in + ! youngest class (used for recruitment) + + !************************************ !** COHORT type structure ** !************************************ @@ -237,6 +247,12 @@ module EDTypesMod real(r8) :: resp_acc real(r8) :: resp_acc_hold + real(r8),allocatable :: frac_leaf_aclass(:) ! This array's size is defined by the + ! number of leaf age classes. + ! and this defines the fraction of + ! leaf mass in different age classes + ! This should ALWAYS sum to 1 !! + real(r8) :: ts_net_uptake(nlevleaf) ! Net uptake of leaf layers: kgC/m2/timestep real(r8) :: year_net_uptake(nlevleaf) ! Net uptake of leaf layers: kgC/m2/year @@ -797,6 +813,7 @@ subroutine dump_cohort(ccohort) write(fates_log(),*) 'co%resp_acc = ', ccohort%resp_acc write(fates_log(),*) 'co%resp_acc_hold = ', ccohort%resp_acc_hold write(fates_log(),*) 'co%rdark = ', ccohort%rdark + write(fates_log(),*) 'co%frac_leaf_aclass = ', ccohort%frac_leaf_aclass(:) write(fates_log(),*) 'co%resp_m = ', ccohort%resp_m write(fates_log(),*) 'co%resp_g = ', ccohort%resp_g write(fates_log(),*) 'co%livestem_mr = ', ccohort%livestem_mr diff --git a/main/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index a710fd922e..afdb49d65d 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -1417,6 +1417,7 @@ subroutine create_patchcohort_structure(this, nc, nsites, sites, bc_in) use EDTypesMod, only : ed_patch_type use EDTypesMod, only : ncwd use EDTypesMod, only : maxSWb + use EDTypesMod, only : nan_leaf_aclass use FatesInterfaceMod, only : fates_maxElementsPerPatch use EDTypesMod, only : maxpft use EDTypesMod, only : area @@ -1553,7 +1554,7 @@ subroutine create_patchcohort_structure(this, nc, nsites, sites, bc_in) call create_cohort(sites(s),newp, temp_cohort%pft, temp_cohort%n, temp_cohort%hite, temp_cohort%dbh, & b_leaf, b_fineroot, b_sapwood, b_dead, b_store, & temp_cohort%laimemory, cohortstatus,recruitstatus, temp_cohort%canopy_trim, newp%NCL_p, & - site_spread, bc_in(s)) + site_spread, nan_leaf_aclass, bc_in(s)) deallocate(temp_cohort) @@ -1819,6 +1820,10 @@ subroutine get_restart_vectors(this, nc, nsites, sites) ccohort%fmort = rio_fmort_co(io_idx_co) ccohort%frmort = rio_frmort_co(io_idx_co) + ! THE ACTUAL LEAF AGES MUST BE READ IN AFTER NEW HYDRO PR + ! (THAT PR HAS THE VECTOR READ-IN MACHINERY) + ccohort%frac_leaf_aclass(1:nleafage) = 1._r8 / real(nleafage,r8) + !Logging ccohort%lmort_direct = rio_lmort_direct_co(io_idx_co) ccohort%lmort_collateral = rio_lmort_collateral_co(io_idx_co) diff --git a/parameter_files/fates_params_14pfts.cdl b/parameter_files/fates_params_14pfts.cdl index 39f1d6ee40..2a9396da80 100644 --- a/parameter_files/fates_params_14pfts.cdl +++ b/parameter_files/fates_params_14pfts.cdl @@ -5,6 +5,7 @@ dimensions: fates_history_size_bins = 13 ; fates_history_height_bins = 6 ; fates_hydr_organs = 4 ; + fates_leafage_class = 1; fates_NCWD = 4 ; fates_litterclass = 6 ; fates_scalar = 1 ; @@ -344,7 +345,7 @@ variables: float fates_leaf_jmaxse(fates_pft) ; fates_leaf_jmaxse:units = "J/mol/K" ; fates_leaf_jmaxse:long_name = "entropy term for jmax" ; - float fates_leaf_long(fates_pft) ; + float fates_leaf_long(fates_leafage_class,fates_pft) ; fates_leaf_long:units = "yr" ; fates_leaf_long:long_name = "Leaf longevity (ie turnover timescale)" ; float fates_leaf_slamax(fates_pft) ; @@ -365,7 +366,7 @@ variables: float fates_leaf_tpuse(fates_pft) ; fates_leaf_tpuse:units = "J/mol/K" ; fates_leaf_tpuse:long_name = "entropy term for tpu" ; - float fates_leaf_vcmax25top(fates_pft) ; + float fates_leaf_vcmax25top(fates_leafage_class,fates_pft) ; fates_leaf_vcmax25top:units = "umol CO2/m^2/s" ; fates_leaf_vcmax25top:long_name = "maximum carboxylation rate of Rub. at 25C, canopy top" ; float fates_leaf_vcmaxha(fates_pft) ; diff --git a/parameter_files/fates_params_default.cdl b/parameter_files/fates_params_default.cdl index 179eefa42d..47c7e55900 100644 --- a/parameter_files/fates_params_default.cdl +++ b/parameter_files/fates_params_default.cdl @@ -6,6 +6,7 @@ dimensions: fates_history_size_bins = 13 ; fates_hydr_organs = 4 ; fates_prt_organs = 6 ; + fates_leafage_class = 4; fates_litterclass = 6 ; fates_pft = 2 ; fates_scalar = 1 ; @@ -287,7 +288,7 @@ variables: float fates_leaf_jmaxse(fates_pft) ; fates_leaf_jmaxse:units = "J/mol/K" ; fates_leaf_jmaxse:long_name = "entropy term for jmax" ; - float fates_leaf_long(fates_pft) ; + float fates_leaf_long(fates_leafage_class, fates_pft) ; fates_leaf_long:units = "yr" ; fates_leaf_long:long_name = "Leaf longevity (ie turnover timescale)" ; float fates_leaf_slamax(fates_pft) ; @@ -308,7 +309,7 @@ variables: float fates_leaf_tpuse(fates_pft) ; fates_leaf_tpuse:units = "J/mol/K" ; fates_leaf_tpuse:long_name = "entropy term for tpu" ; - float fates_leaf_vcmax25top(fates_pft) ; + float fates_leaf_vcmax25top(fates_leafage_class, fates_pft) ; fates_leaf_vcmax25top:units = "umol CO2/m^2/s" ; fates_leaf_vcmax25top:long_name = "maximum carboxylation rate of Rub. at 25C, canopy top" ; float fates_leaf_vcmaxha(fates_pft) ; @@ -878,7 +879,10 @@ data: fates_leaf_jmaxse = 495, 495 ; - fates_leaf_long = 1.5, 1.5 ; + fates_leaf_long = 0.4, 0.4, + 0.4, 0.4, + 0.4, 0.4, + 0.4, 0.4 ; fates_leaf_slamax = 0.0954, 0.0954 ; @@ -892,7 +896,10 @@ data: fates_leaf_tpuse = 490, 490 ; - fates_leaf_vcmax25top = 50, 50 ; + fates_leaf_vcmax25top = 62, 62, + 54, 54, + 46, 46, + 38, 38; fates_leaf_vcmaxha = 65330, 65330 ; From cbdc5e87181f5a15061507804ad4bb5b90e95c36 Mon Sep 17 00:00:00 2001 From: chonggang xu Date: Sat, 15 Dec 2018 14:22:21 -0800 Subject: [PATCH 07/34] Adding componenents of discrete age bins (transcribed by rknox) --- biogeochem/EDCanopyStructureMod.F90 | 3 + biogeochem/EDCohortDynamicsMod.F90 | 108 +++++++++++++++------ biogeochem/EDPatchDynamicsMod.F90 | 4 + biogeochem/EDPhysiologyMod.F90 | 6 +- biogeochem/FatesAllometryMod.F90 | 2 +- biogeophys/FatesPlantRespPhotosynthMod.F90 | 30 +++--- main/FatesInterfaceMod.F90 | 22 +++-- main/FatesInventoryInitMod.F90 | 5 +- main/FatesParameterDerivedMod.F90 | 55 +++++------ main/FatesParametersInterface.F90 | 1 + parteh/PRTLossFluxesMod.F90 | 4 +- 11 files changed, 158 insertions(+), 82 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index caf6861577..50b034026d 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -25,6 +25,7 @@ module EDCanopyStructureMod use FatesInterfaceMod , only : hlm_days_per_year use FatesInterfaceMod , only : hlm_use_planthydro use FatesInterfaceMod , only : numpft + use FatesInterfaceMod , only : nleafage use FatesPlantHydraulicsMod, only : UpdateH2OVeg,InitHydrCohort use PRTGenericMod, only : leaf_organ @@ -503,6 +504,7 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) ! demoted to the understory allocate(copyc) + allocate(copyc%frac_leaf_aclass(nleafage)) call InitPRTCohort(copyc) if( hlm_use_planthydro.eq.itrue ) then call InitHydrCohort(currentSite,copyc) @@ -850,6 +852,7 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) elseif ( cc_gain > nearzero .and. cc_gain < currentCohort%c_area) then allocate(copyc) + allocate(copyc%frac_leaf_aclass(nleafage)) call InitPRTCohort(copyc) if( hlm_use_planthydro.eq.itrue ) then call InitHydrCohort(CurrentSite,copyc) diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 71a1416cfc..2a7284d323 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -12,6 +12,7 @@ module EDCohortDynamicsMod use FatesConstantsMod , only : fates_unset_int use FatesConstantsMod , only : itrue,ifalse use FatesInterfaceMod , only : hlm_days_per_year + use FatesInterfaceMod , only : nleafage use EDPftvarcon , only : EDPftvarcon_inst use EDTypesMod , only : ed_site_type, ed_patch_type, ed_cohort_type use EDTypesMod , only : nclmax @@ -60,7 +61,7 @@ module EDCohortDynamicsMod use PRTAllometricCarbonMod, only : ac_bc_in_id_ctrim use PRTAllometricCarbonMod, only : ac_bc_inout_id_dbh - + use shr_infnan_mod, only : nan => shr_infnan_nan, assignment(=) ! CIME globals use shr_log_mod , only : errMsg => shr_log_errMsg @@ -91,34 +92,38 @@ module EDCohortDynamicsMod !-------------------------------------------------------------------------------------! - subroutine create_cohort(currentSite, patchptr, pft, nn, hite, dbh, bleaf, bfineroot, bsap, & - bdead, bstore, laimemory, status, recruitstatus,ctrim, clayer, spread, bc_in) + subroutine create_cohort(currentSite, patchptr, pft, nn, hite, dbh, bleaf, bfineroot, & + bsap, bdead, bstore, laimemory, status, recruitstatus,ctrim, & + clayer, spread, leaf_aclass_init, bc_in) + - ! - ! !DESCRIPTION: - ! create new cohort - ! - ! !USES: - ! - ! !ARGUMENTS type(ed_site_type), intent(inout), target :: currentSite type(ed_patch_type), intent(inout), pointer :: patchptr - integer, intent(in) :: pft ! Cohort Plant Functional Type - integer, intent(in) :: clayer ! canopy status of cohort (1 = canopy, 2 = understorey, etc.) - integer, intent(in) :: status ! growth status of plant (2 = leaves on , 1 = leaves off) - integer, intent(in) :: recruitstatus ! recruit status of plant (1 = recruitment , 0 = other) - real(r8), intent(in) :: nn ! number of individuals in cohort per 'area' (10000m2 default) - real(r8), intent(in) :: hite ! height: meters - real(r8), intent(in) :: dbh ! dbh: cm - real(r8), intent(in) :: bleaf ! biomass in leaves: kgC - real(r8), intent(in) :: bfineroot ! biomass in fineroots: kgC - real(r8), intent(in) :: bsap ! biomass in sapwood: kgC - real(r8), intent(in) :: bdead ! total dead biomass: kGC per indiv - real(r8), intent(in) :: bstore ! stored carbon: kGC per indiv - real(r8), intent(in) :: laimemory ! target leaf biomass- set from previous year: kGC per indiv - real(r8), intent(in) :: ctrim ! What is the fraction of the maximum leaf biomass that we are targeting? :- - real(r8), intent(in) :: spread ! The community assembly effects how spread crowns are in horizontal space - type(bc_in_type), intent(in) :: bc_in ! External boundary conditions + integer, intent(in) :: pft ! Cohort Plant Functional Type + integer, intent(in) :: clayer ! canopy status of cohort + ! (1 = canopy, 2 = understorey, etc.) + integer, intent(in) :: status ! growth status of plant + ! (2 = leaves on , 1 = leaves off) + integer, intent(in) :: recruitstatus ! recruit status of plant + ! (1 = recruitment , 0 = other) + real(r8), intent(in) :: nn ! number of individuals in cohort + ! per 'area' (10000m2 default) + real(r8), intent(in) :: hite ! height: meters + real(r8), intent(in) :: dbh ! dbh: cm + real(r8), intent(in) :: bleaf ! biomass in leaves: kgC + real(r8), intent(in) :: bfineroot ! biomass in fineroots: kgC + real(r8), intent(in) :: bsap ! biomass in sapwood: kgC + real(r8), intent(in) :: bdead ! total dead biomass: kGC per indiv + real(r8), intent(in) :: bstore ! stored carbon: kGC per indiv + real(r8), intent(in) :: laimemory ! target leaf biomass- set from + ! previous year: kGC per indiv + real(r8), intent(in) :: ctrim ! What is the fraction of the maximum + ! leaf biomass that we are targeting? + real(r8), intent(in) :: spread ! The community assembly effects how + ! spread crowns are in horizontal space + integer, intent(in) :: leaf_aclass_init ! how to initialized the leaf age class + ! distribution + type(bc_in_type), intent(in) :: bc_in ! External boundary conditions ! ! !LOCAL VARIABLES: @@ -129,6 +134,7 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, hite, dbh, bleaf, bfine !---------------------------------------------------------------------- allocate(new_cohort) + allocate(new_cohort%frac_leaf_aclass(nleafage)) call nan_cohort(new_cohort) ! Make everything in the cohort not-a-number call zero_cohort(new_cohort) ! Zero things that need to be zeroed. @@ -152,6 +158,22 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, hite, dbh, bleaf, bfine new_cohort%canopy_layer_yesterday = real(clayer, r8) new_cohort%laimemory = laimemory + + ! All newly initialized cohorts start off with all leaf biomass + ! as classified in the newest age class + if(leaf_aclass_init .eq. equal_leaf_aclass) then + new_cohort%frac_leaf_aclass(1:nleafage) = 1._r8 / real(nleafage,r8) + elseif(leaf_aclass_init .eq. first_leaf_aclass) then + new_cohort%frac_leaf_aclass(1:nleafage) = 0._r8 + new_cohort%frac_leaf_aclass(1) = 1._r8 + elseif(leaf_aclass_init .eq. nan_leaf_aclass) then + new_cohort%frac_leaf_aclass(1:nleafage) = nan + else + write(fates_log(),*) 'An unknown leaf age distribution was' + write(fates_log(),*) 'requested during create cohort' + write(fates_log(),*) 'leaf_aclass_init: ',leaf_aclass_init + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if ! Initialize the Plant allocative Reactive Transport (PaRT) module ! Choose from one of the extensible hypotheses (EH) @@ -343,7 +365,7 @@ subroutine nan_cohort(cc_p) ! Make all the cohort variables NaN so they aren't used before defined. ! ! !USES: - use shr_infnan_mod, only : nan => shr_infnan_nan, assignment(=) + use FatesConstantsMod, only : fates_unset_int ! @@ -389,7 +411,8 @@ subroutine nan_cohort(cc_p) currentCohort%c_area = nan ! areal extent of canopy (m2) currentCohort%treelai = nan ! lai of tree (total leaf area (m2) / canopy area (m2) currentCohort%treesai = nan ! stem area index of tree (total stem area (m2) / canopy area (m2) - + currentCohort%frac_leaf_aclass(:)= nan ! leaf age classes + ! CARBON FLUXES currentCohort%gpp_acc_hold = nan ! GPP: kgC/indiv/year currentCohort%gpp_tstep = nan ! GPP: kgC/indiv/timestep @@ -686,6 +709,7 @@ subroutine terminate_cohorts( currentSite, currentPatch, level ) ! Deallocate the cohort's PRT structure call currentCohort%prt%DeallocatePRTVartypes() deallocate(currentCohort%prt) + deallocate(currentCohort%frac_leaf_aclass) deallocate(currentCohort) nullify(currentCohort) @@ -706,7 +730,6 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) ! ! !USES: use EDParamsMod , only : ED_val_cohort_fusion_tol - use shr_infnan_mod, only : nan => shr_infnan_nan, assignment(=) ! ! !ARGUMENTS type (ed_site_type), intent(inout), target :: currentSite @@ -728,6 +751,8 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) integer :: nocohorts real(r8) :: newn real(r8) :: diff + real(r8) :: leaf_c_next ! Leaf carbon * plant density of current (for weighting) + real(r8) :: leaf_c_curr ! Leaf carbon * plant density of next (for weighting) real(r8) :: dynamic_fusion_tolerance integer :: largersc, smallersc, sc_i ! indices for tracking the growth flux caused by fusion @@ -851,6 +876,26 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) ! recent canopy history currentCohort%canopy_layer_yesterday = (currentCohort%n*currentCohort%canopy_layer_yesterday + & nextc%n*nextc%canopy_layer_yesterday)/newn + + ! Leaf age class fractions + ! ----------------------------------------------------------------- + leaf_c_next = nextc%prt%GetState(leaf_organ, all_carbon_elements)*nextc%n + leaf_c_curr = currentCohort%prt%GetState(leaf_organ, all_carbon_elements)*currentCohort%n + + if( (leaf_c_next + leaf_c_curr) > nearzero ) then + currentCohort%frac_leaf_aclass(1:nleafage) = & + (leaf_c_curr*currentCohort%frac_leaf_aclass(1:nleafage) + & + leaf_c_next*nextc%frac_leaf_aclass(1:nleafage))/(leaf_c_next+leaf_c_curr) + + ! Force leaf age class to sum to unity + currentCohort%frac_leaf_aclass(1:nleafage) = & + currentCohort%frac_leaf_aclass(1:nleafage) / & + sum(currentCohort%frac_leaf_aclass(1:nleafage)) + else + currentCohort%frac_leaf_aclass(1:nleafage) = 0._r8 + currentCohort%frac_leaf_aclass(1) = 1._r8 + end if + ! keep track of the size class bins so that we can monitor growth fluxes ! compare the values. if they are the same, then nothing needs to be done. if not, track the diagnostic flux @@ -981,7 +1026,7 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) ! Deallocate the cohort's PRT structure call nextc%prt%DeallocatePRTVartypes() deallocate(nextc%prt) - + deallocate(nextc%frac_leaf_aclass) deallocate(nextc) nullify(nextc) @@ -1262,7 +1307,8 @@ subroutine copy_cohort( currentCohort,copyc ) ! This transfers the PRT objects over. call n%prt%CopyPRTVartypes(o%prt) - + + n%frac_leaf_aclass(1:nleafage) = o%frac_leaf_aclass(1:nleafage) ! CARBON FLUXES n%gpp_acc_hold = o%gpp_acc_hold diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 98546557e3..78abbef13c 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -22,6 +22,7 @@ module EDPatchDynamicsMod use FatesInterfaceMod , only : bc_in_type use FatesInterfaceMod , only : hlm_days_per_year use FatesInterfaceMod , only : numpft + use FatesInterfaceMod , only : nleafage use FatesGlobals , only : endrun => fates_endrun use FatesConstantsMod , only : r8 => fates_r8 use FatesConstantsMod , only : itrue @@ -401,6 +402,7 @@ subroutine spawn_patches( currentSite, bc_in) do while(associated(currentCohort)) allocate(nc) + allocate(nc%frac_leaf_aclass(nleafage)) if(hlm_use_planthydro.eq.itrue) call InitHydrCohort(CurrentSite,nc) call InitPRTCohort(nc) call zero_cohort(nc) @@ -674,6 +676,7 @@ subroutine spawn_patches( currentSite, bc_in) if(hlm_use_planthydro.eq.itrue) call DeallocateHydrCohort(nc) call nc%prt%DeallocatePRTVartypes() deallocate(nc%prt) + deallocate(nc%frac_leaf_aclass) deallocate(nc) endif @@ -1926,6 +1929,7 @@ subroutine dealloc_patch(cpatch) if(hlm_use_planthydro.eq.itrue) call DeallocateHydrCohort(ccohort) call ccohort%prt%DeallocatePRTVartypes() deallocate(ccohort%prt) + deallocate(ccohort%frac_leaf_aclass) deallocate(ccohort) ccohort => ncohort diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index ad639950aa..4d2a885ff2 100644 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -12,6 +12,7 @@ module EDPhysiologyMod use FatesInterfaceMod, only : hlm_freq_day use FatesInterfaceMod, only : hlm_day_of_year use FatesInterfaceMod, only : numpft + use FatesInterfaceMod, only : nleafage use FatesInterfaceMod, only : hlm_use_planthydro use FatesConstantsMod, only : r8 => fates_r8 use FatesConstantsMod, only : nearzero @@ -32,6 +33,7 @@ module EDPhysiologyMod use EDTypesMod , only : maxpft use EDTypesMod , only : ed_site_type, ed_patch_type, ed_cohort_type use EDTypesMod , only : dump_cohort + use EDTypesMod , only : first_leaf_aclass use shr_log_mod , only : errMsg => shr_log_errMsg use FatesGlobals , only : fates_log @@ -326,7 +328,7 @@ subroutine trim_canopy( currentSite ) else !evergreen costs ! Leaf cost at leaf level z accounting for sla profile currentCohort%leaf_cost = 1.0_r8/(sla_levleaf* & - EDPftvarcon_inst%leaf_long(ipft)*1000.0_r8) !convert from sla in m2g-1 to m2kg-1 + EDPftvarcon_inst%leaf_long(ipft,0)*1000.0_r8) !convert from sla in m2g-1 to m2kg-1 if ( int(EDPftvarcon_inst%allom_fmode(ipft)) .eq. 1 ) then ! if using trimmed leaf for fine root biomass allometry, add the cost of the root increment ! to the leaf increment; otherwise do not. @@ -1027,7 +1029,7 @@ subroutine recruitment( currentSite, currentPatch, bc_in ) call create_cohort(currentSite,currentPatch, temp_cohort%pft, temp_cohort%n, temp_cohort%hite, temp_cohort%dbh, & b_leaf, b_fineroot, b_sapwood, b_dead, b_store, & temp_cohort%laimemory, cohortstatus,recruitstatus, temp_cohort%canopy_trim, currentPatch%NCL_p, & - currentSite%spread, bc_in) + currentSite%spread, first_leaf_aclass, bc_in) ! keep track of how many individuals were recruited for passing to history diff --git a/biogeochem/FatesAllometryMod.F90 b/biogeochem/FatesAllometryMod.F90 index 7bb2f0ec9c..3c03524f3f 100644 --- a/biogeochem/FatesAllometryMod.F90 +++ b/biogeochem/FatesAllometryMod.F90 @@ -2103,7 +2103,7 @@ real(r8) function decay_coeff_kn(pft) ! kn = 0.11. Here, we derive kn from vcmax25 as in Lloyd et al ! (2010) Biogeosciences, 7, 1833-1859 - decay_coeff_kn = exp(0.00963_r8 * EDPftvarcon_inst%vcmax25top(pft) - 2.43_r8) + decay_coeff_kn = exp(0.00963_r8 * EDPftvarcon_inst%vcmax25top(pft,0) - 2.43_r8) return end function decay_coeff_kn diff --git a/biogeophys/FatesPlantRespPhotosynthMod.F90 b/biogeophys/FatesPlantRespPhotosynthMod.F90 index 2fae4eceb2..c9b022ee30 100644 --- a/biogeophys/FatesPlantRespPhotosynthMod.F90 +++ b/biogeophys/FatesPlantRespPhotosynthMod.F90 @@ -385,14 +385,17 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) ! cohort-layer combo of interest. ! but in the vanilla case, we only re-calculate if it has ! not been done yet. + ! Other cases where we need to solve for every cohort + ! in every leaf layer: nutrient dynamic mode, multiple leaf + ! age classes ! ------------------------------------------------------------ if ( .not.rate_mask_z(iv,ft,cl) .or. & (hlm_use_planthydro.eq.itrue) .or. & + (nleafage > 1) .or. & (hlm_parteh_mode .ne. prt_carbon_allom_hyp ) ) then - if (hlm_use_planthydro.eq.itrue) then - + if (hlm_use_planthydro.eq.itrue ) then bbb = max (bbbopt(nint(c3psn(ft)))*currentCohort%co_hydr%btran(1), 1._r8) btran_eff = currentCohort%co_hydr%btran(1) @@ -476,10 +479,11 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) call LeafLayerBiophysicalRates(currentPatch%ed_parsun_z(cl,ft,iv), & ! in ft, & ! in - EDPftvarcon_inst%vcmax25top(ft), & ! in - param_derived%jmax25top(ft), & ! in - param_derived%tpu25top(ft), & ! in - param_derived%kp25top(ft), & ! in + currentCohort%frac_leaf_aclass(:), & ! in + EDPftvarcon_inst%vcmax25top(ft,:), & ! in + param_derived%jmax25top(ft,:), & ! in + param_derived%tpu25top(ft,:), & ! in + param_derived%kp25top(ft,:), & ! in nscaler, & ! in bc_in(s)%t_veg_pa(ifp), & ! in btran_eff, & ! in @@ -1701,6 +1705,7 @@ end subroutine LeafLayerMaintenanceRespiration subroutine LeafLayerBiophysicalRates( parsun_lsl, & ft, & + frac_leaf_aclass, & vcmax25top_ft, & jmax25top_ft, & tpu25top_ft, & @@ -1736,11 +1741,12 @@ subroutine LeafLayerBiophysicalRates( parsun_lsl, & real(r8), intent(in) :: parsun_lsl ! PAR absorbed in sunlit leaves for this layer integer, intent(in) :: ft ! (plant) Functional Type Index real(r8), intent(in) :: nscaler ! Scale for leaf nitrogen profile - real(r8), intent(in) :: vcmax25top_ft ! canopy top maximum rate of carboxylation at 25C + real(r8), intent(in) :: frac_leaf_aclass(nleafage) ! Fraction of leaves in each age-class + real(r8), intent(in) :: vcmax25top_ft(nleafage) ! canopy top maximum rate of carboxylation at 25C ! for this pft (umol CO2/m**2/s) - real(r8), intent(in) :: jmax25top_ft ! canopy top maximum electron transport rate at 25C + real(r8), intent(in) :: jmax25top_ft(nleafage) ! canopy top maximum electron transport rate at 25C ! for this pft (umol electrons/m**2/s) - real(r8), intent(in) :: tpu25top_ft ! canopy top triose phosphate utilization rate at 25C + real(r8), intent(in) :: tpu25top_ft(nleafage) ! canopy top triose phosphate utilization rate at 25C ! for this pft (umol CO2/m**2/s) real(r8), intent(in) :: co2_rcurve_islope25top_ft ! initial slope of CO2 response curve ! (C4 plants) at 25C, canopy top, this pft @@ -1803,9 +1809,9 @@ subroutine LeafLayerBiophysicalRates( parsun_lsl, & tpu = 0._r8 co2_rcurve_islope = 0._r8 else ! day time - vcmax25 = vcmax25top_ft * nscaler - jmax25 = jmax25top_ft * nscaler - tpu25 = tpu25top_ft * nscaler + vcmax25 = sum(vcmax25top_ft(1:nleafage) * frac_leaf_aclass(1:nleafage)) * nscaler + jmax25 = sum(jmax25top_ft(1:nleafage) * frac_leaf_aclass(1:nleafage)) * nscaler + tpu25 = sum(tpu25top_ft(1:nleafage) * frac_leaf_aclass(1:nleafage)) * nscaler co2_rcurve_islope25 = co2_rcurve_islope25top_ft * nscaler ! Adjust for temperature diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index fb5002d2ac..5c7586db17 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -249,11 +249,11 @@ module FatesInterfaceMod ! ! ------------------------------------------------------------------------------------- - integer, protected :: numpft ! The total number of PFTs defined in the simulation - integer, protected :: nlevsclass ! The total number of cohort size class bins output to history - integer, protected :: nlevage ! The total number of patch age bins output to history - integer, protected :: nlevheight ! The total number of height bins output to history - + integer, protected :: numpft ! The total number of PFTs defined in the simulation + integer, protected :: nlevsclass ! The total number of cohort size class bins output to history + integer, protected :: nlevage ! The total number of patch age bins output to history + integer, protected :: nlevheight ! The total number of height bins output to history + integer, protected :: nleafage ! The total number of leaf age classes ! ------------------------------------------------------------------------------------- ! Structured Boundary Conditions (SITE/PATCH SCALE) @@ -965,8 +965,18 @@ subroutine set_fates_global_elements(use_fates) write(fates_log(), *) 'FatesInterfaceMod.F90:maxpft accordingly' call endrun(msg=errMsg(sourcefile, __LINE__)) end if - + ! Identify the number of leaf age-classes + + if( (lbound(EDPftvarcon_inst%leaf_long(:,:),dim=2) .eq. 0) .or. & + (ubound(EDPftvarcon_inst%leaf_long(:,:),dim=2) .eq. 0) ) then + write(fates_log(), *) 'While assessing the number of FATES leaf age classes,' + write(fates_log(), *) 'The second dimension of leaf_long was 0?' + call endrun(msg=errMsg(sourcefile, __LINE__)) + else + nleafage = size(EDPftvarcon_inst%leaf_long,dim=2) + end if + ! These values are used to define the restart file allocations and general structure ! of memory for the cohort arrays diff --git a/main/FatesInventoryInitMod.F90 b/main/FatesInventoryInitMod.F90 index 5c682a9f02..657d00beab 100644 --- a/main/FatesInventoryInitMod.F90 +++ b/main/FatesInventoryInitMod.F90 @@ -33,6 +33,7 @@ module FatesInventoryInitMod use EDTypesMod , only : ed_patch_type use EDTypesMod , only : ed_cohort_type use EDTypesMod , only : area + use EDTypesMod , only : equal_leaf_aclass use EDPftvarcon , only : EDPftvarcon_inst @@ -880,6 +881,7 @@ 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 + ! Don't need to allocate leaf age classes (not used) temp_cohort%pft = c_pft temp_cohort%n = c_nplant * cpatch%area @@ -934,7 +936,8 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & ! Since spread is a canopy level calculation, we need to provide an initial guess here. call create_cohort(csite, cpatch, c_pft, temp_cohort%n, temp_cohort%hite, temp_cohort%dbh, & b_leaf, b_fineroot, b_sapwood, b_dead, b_store, & - temp_cohort%laimemory, cstatus, rstatus, temp_cohort%canopy_trim, 1, csite%spread, bc_in) + temp_cohort%laimemory, cstatus, rstatus, temp_cohort%canopy_trim, & + 1, csite%spread, equal_leaf_aclass, bc_in) deallocate(temp_cohort) ! get rid of temporary cohort diff --git a/main/FatesParameterDerivedMod.F90 b/main/FatesParameterDerivedMod.F90 index 9da1eeb01e..f175114acc 100644 --- a/main/FatesParameterDerivedMod.F90 +++ b/main/FatesParameterDerivedMod.F90 @@ -12,17 +12,16 @@ module FatesParameterDerivedMod use FatesConstantsMod, only : r8 => fates_r8 use FatesConstantsMod, only : umolC_to_kgC use FatesConstantsMod, only : g_per_kg + use FatesInterfaceMod, only : nleafage type param_derived_type - real(r8), allocatable :: vcmax25top(:) ! canopy top: maximum rate of carboxylation - ! at 25C (umol CO2/m**2/s) - real(r8), allocatable :: jmax25top(:) ! canopy top: maximum electron transport - ! rate at 25C (umol electrons/m**2/s) - real(r8), allocatable :: tpu25top(:) ! canopy top: triose phosphate utilization - ! rate at 25C (umol CO2/m**2/s) - real(r8), allocatable :: kp25top(:) ! canopy top: initial slope of CO2 response - ! curve (C4 plants) at 25C + real(r8), allocatable :: jmax25top(:,:) ! canopy top: maximum electron transport + ! rate at 25C (umol electrons/m**2/s) + real(r8), allocatable :: tpu25top(:,:) ! canopy top: triose phosphate utilization + ! rate at 25C (umol CO2/m**2/s) + real(r8), allocatable :: kp25top(:,:) ! canopy top: initial slope of CO2 response + ! curve (C4 plants) at 25C contains procedure :: Init @@ -39,10 +38,9 @@ subroutine InitAllocate(this,numpft) class(param_derived_type), intent(inout) :: this integer, intent(in) :: numpft - allocate(this%vcmax25top(numpft)) - allocate(this%jmax25top(numpft)) - allocate(this%tpu25top(numpft)) - allocate(this%kp25top(numpft)) + allocate(this%jmax25top(numpft,nleafage)) + allocate(this%tpu25top(numpft,nleafage)) + allocate(this%kp25top(numpft,nleafage)) return end subroutine InitAllocate @@ -58,6 +56,7 @@ subroutine Init(this,numpft) ! local variables integer :: ft ! pft index + integer :: iage ! leaf age class index real(r8) :: lnc ! leaf N concentration (gN leaf/m^2) associate( vcmax25top => EDPftvarcon_inst%vcmax25top ) @@ -66,22 +65,24 @@ subroutine Init(this,numpft) do ft = 1,numpft - ! Parameters derived from vcmax25top. - ! Bonan et al (2011) JGR, 116, doi:10.1029/2010JG001593 - ! used jmax25 = 1.97 vcmax25, from Wullschleger (1993) Journal of - ! Experimental Botany 44:907-920. Here use a factor "1.67", from - ! Medlyn et al (2002) Plant, Cell and Environment 25:1167-1179 - - ! RF - copied this from the CLM trunk code, but where did it come from, - ! and how can we make these consistant? - ! jmax25top(ft) = & - ! (2.59_r8 - 0.035_r8*min(max((t10(p)-tfrzc),11._r8),35._r8)) * vcmax25top(ft) - - this%jmax25top(ft) = 1.67_r8 * vcmax25top(ft) - this%tpu25top(ft) = 0.167_r8 * vcmax25top(ft) - this%kp25top(ft) = 20000._r8 * vcmax25top(ft) - + do iage = 1, nleafage + + ! Parameters derived from vcmax25top. + ! Bonan et al (2011) JGR, 116, doi:10.1029/2010JG001593 + ! used jmax25 = 1.97 vcmax25, from Wullschleger (1993) Journal of + ! Experimental Botany 44:907-920. Here use a factor "1.67", from + ! Medlyn et al (2002) Plant, Cell and Environment 25:1167-1179 + + ! RF - copied this from the CLM trunk code, but where did it come from, + ! and how can we make these consistant? + ! jmax25top(ft) = & + ! (2.59_r8 - 0.035_r8*min(max((t10(p)-tfrzc),11._r8),35._r8)) * vcmax25top(ft) + + this%jmax25top(ft,iage) = 1.67_r8 * vcmax25top(ft,iage) + this%tpu25top(ft,iage) = 0.167_r8 * vcmax25top(ft,iage) + this%kp25top(ft,iage) = 20000._r8 * vcmax25top(ft,iage) + end do end do !ft diff --git a/main/FatesParametersInterface.F90 b/main/FatesParametersInterface.F90 index 8af57330dc..c863542808 100644 --- a/main/FatesParametersInterface.F90 +++ b/main/FatesParametersInterface.F90 @@ -30,6 +30,7 @@ module FatesParametersInterface character(len=*), parameter, public :: dimension_name_variants = 'fates_variants' character(len=*), parameter, public :: dimension_name_hydr_organs = 'fates_hydr_organs' character(len=*), parameter, public :: dimension_name_prt_organs = 'fates_prt_organs' + character(len=*), parameter, public :: dimension_name_leaf_age = 'fates_leafage_class' character(len=*), parameter, public :: dimension_name_history_size_bins = 'fates_history_size_bins' character(len=*), parameter, public :: dimension_name_history_age_bins = 'fates_history_age_bins' character(len=*), parameter, public :: dimension_name_history_height_bins = 'fates_history_height_bins' diff --git a/parteh/PRTLossFluxesMod.F90 b/parteh/PRTLossFluxesMod.F90 index e91ee09bf2..57d727206b 100644 --- a/parteh/PRTLossFluxesMod.F90 +++ b/parteh/PRTLossFluxesMod.F90 @@ -614,9 +614,9 @@ subroutine MaintTurnoverSimpleRetranslocation(prt,ipft) ! Only EVERGREENS HAVE MAINTENANCE LEAF TURNOVER ! ------------------------------------------------------------------------------------- - if ( (EDPftvarcon_inst%leaf_long(ipft) > nearzero ) .and. & + if ( (EDPftvarcon_inst%leaf_long(ipft,0) > nearzero ) .and. & (EDPftvarcon_inst%evergreen(ipft) == 1) ) then - base_turnover(leaf_organ) = hlm_freq_day / EDPftvarcon_inst%leaf_long(ipft) + base_turnover(leaf_organ) = hlm_freq_day / EDPftvarcon_inst%leaf_long(ipft,0) else base_turnover(leaf_organ) = 0.0_r8 endif From a3767195290b03204bdff4147c0a59dd06ecfa65 Mon Sep 17 00:00:00 2001 From: chonggang xu Date: Sun, 16 Dec 2018 16:23:59 -0800 Subject: [PATCH 08/34] leaf age discretization --- biogeochem/EDCohortDynamicsMod.F90 | 11 +++- biogeochem/EDPhysiologyMod.F90 | 98 ++++++++++++++++++------------ parteh/PRTLossFluxesMod.F90 | 28 ++++++--- 3 files changed, 86 insertions(+), 51 deletions(-) diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 2a7284d323..0129fce997 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -159,8 +159,10 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, hite, dbh, bleaf, bfine new_cohort%laimemory = laimemory - ! All newly initialized cohorts start off with all leaf biomass - ! as classified in the newest age class + ! All newly initialized cohorts start off with an assumption + ! about leaf age (depending on what is calling the initialization + ! of this cohort + if(leaf_aclass_init .eq. equal_leaf_aclass) then new_cohort%frac_leaf_aclass(1:nleafage) = 1._r8 / real(nleafage,r8) elseif(leaf_aclass_init .eq. first_leaf_aclass) then @@ -878,6 +880,8 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) nextc%n*nextc%canopy_layer_yesterday)/newn ! Leaf age class fractions + ! First calculate the total amount of leaf in each to enable + ! a weighted average of the two ! ----------------------------------------------------------------- leaf_c_next = nextc%prt%GetState(leaf_organ, all_carbon_elements)*nextc%n leaf_c_curr = currentCohort%prt%GetState(leaf_organ, all_carbon_elements)*currentCohort%n @@ -892,6 +896,8 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) currentCohort%frac_leaf_aclass(1:nleafage) / & sum(currentCohort%frac_leaf_aclass(1:nleafage)) else + ! If there is no leaf at all... then we set + ! the age fraction to be all in the growing class currentCohort%frac_leaf_aclass(1:nleafage) = 0._r8 currentCohort%frac_leaf_aclass(1) = 1._r8 end if @@ -1308,6 +1314,7 @@ subroutine copy_cohort( currentCohort,copyc ) ! This transfers the PRT objects over. call n%prt%CopyPRTVartypes(o%prt) + ! Leaf age class fractions n%frac_leaf_aclass(1:nleafage) = o%frac_leaf_aclass(1:nleafage) ! CARBON FLUXES diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index 4d2a885ff2..a07d324408 100644 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -34,6 +34,8 @@ module EDPhysiologyMod use EDTypesMod , only : ed_site_type, ed_patch_type, ed_cohort_type use EDTypesMod , only : dump_cohort use EDTypesMod , only : first_leaf_aclass + use EDTypesMod , only : leaves_on + use EDTypesMod , only : leaves_off use shr_log_mod , only : errMsg => shr_log_errMsg use FatesGlobals , only : fates_log @@ -491,11 +493,11 @@ subroutine phenology( currentSite, bc_in ) !2) The leaves should not be on already !3) There should have been at least on chilling day in the counting period. if (currentSite%ED_GDD_site > gdd_threshold)then - if (currentSite%status == 1) then + if (currentSite%status == leaves_off) then if (currentSite%ncd >= 1) then - currentSite%status = 2 !alter status of site to 'leaves on' - ! NOTE(bja, 2015-01) should leafondate = model_day to be consistent with leaf off? - currentSite%leafondate = t !record leaf on date + currentSite%status = leaves_on ! alter status of site to 'leaves on' + ! NOTE(bja, 2015-01) should leafondate = model_day to be consistent with leaf off? + currentSite%leafondate = t !record leaf on date if ( debug ) write(fates_log(),*) 'leaves on' endif !ncd endif !status @@ -513,8 +515,8 @@ subroutine phenology( currentSite, bc_in ) if (ncolddays > ED_val_phen_ncolddayslim)then if (timesinceleafon > ED_val_phen_mindayson)then - if (currentSite%status == 2)then - currentSite%status = 1 !alter status of site to 'leaves on' + if (currentSite%status == leaves_on)then + currentSite%status = leaves_off !alter status of site to 'leaves off' currentSite%leafoffdate = hlm_model_day !record leaf off date if ( debug ) write(fates_log(),*) 'leaves off' endif @@ -523,8 +525,8 @@ subroutine phenology( currentSite, bc_in ) !LEAF OFF: COLD LIFESPAN THRESHOLD if(timesinceleafoff > 400)then !remove leaves after a whole year when there is no 'off' period. - if(currentSite%status == 2)then - currentSite%status = 1 !alter status of site to 'leaves on' + if(currentSite%status == leaves_on)then + currentSite%status = leaves_off !alter status of site to 'leaves off' currentSite%leafoffdate = hlm_model_day !record leaf off date if ( debug ) write(fates_log(),*) 'leaves off' endif @@ -566,7 +568,7 @@ subroutine phenology( currentSite, bc_in ) !In drought phenology, we often need to force the leaves to stay on or off as moisture fluctuates... timesincedleafoff = 0 - if (currentSite%dstatus == 1)then !the leaves are off. How long have they been off? + if (currentSite%dstatus == leaves_off)then !the leaves are off. How long have they been off? !leaves have come on, but last year, so at a later date than now. if (currentSite%dleafoffdate > 0.and.currentSite%dleafoffdate > t)then timesincedleafoff = t + (360 - currentSite%dleafoffdate) @@ -577,7 +579,7 @@ subroutine phenology( currentSite, bc_in ) timesincedleafon = 0 !the leaves are on. How long have they been on? - if (currentSite%dstatus == 2)then + if (currentSite%dstatus == leaves_on)then !leaves have come on, but last year, so at a later date than now. if (currentSite%dleafondate > 0.and.currentSite%dleafondate > t)then timesincedleafon = t + (360 - currentSite%dleafondate) @@ -592,37 +594,37 @@ subroutine phenology( currentSite, bc_in ) currentSite%dleafondate < 15))then ! are we in the window? ! TODO: CHANGE THIS MATH, MOVE THE DENOMENATOR OUTSIDE OF THE SUM (rgk 01-2017) if (sum(currentSite%water_memory(1:numWaterMem)/real(numWaterMem,r8)) & - >= ED_val_phen_drought_threshold.and.currentSite%dstatus == 1.and.t >= 10)then + >= ED_val_phen_drought_threshold.and.currentSite%dstatus == leaves_off.and.t >= 10)then ! leave some minimum time between leaf off and leaf on to prevent 'flickering'. if (timesincedleafoff > ED_val_phen_doff_time)then - currentSite%dstatus = 2 !alter status of site to 'leaves on' - currentSite%dleafondate = t !record leaf on date + currentSite%dstatus = leaves_on !alter status of site to 'leaves on' + currentSite%dleafondate = t !record leaf on date endif endif endif !we still haven't done budburst by end of window - if (t == currentSite%dleafondate+30.and.currentSite%dstatus == 1)then - currentSite%dstatus = 2 ! force budburst! + if (t == currentSite%dleafondate+30.and.currentSite%dstatus == leaves_off)then + currentSite%dstatus = leaves_on ! force budburst! currentSite%dleafondate = t ! record leaf on date endif !LEAF OFF: DROUGHT DECIDUOUS LIFESPAN - if the leaf gets to the end of its useful life. A*, E* - if (currentSite%dstatus == 2.and.t >= 10)then !D* + if (currentSite%dstatus == leaves_on.and.t >= 10)then !D* !Are the leaves at the end of their lives? !FIX(RF,0401014)- this is hardwiring.... !FIX(RGK:changed from hard-coded pft 7 leaf lifespan to labeled constant (1 year) if ( timesincedleafon > canopy_leaf_lifespan )then - currentSite%dstatus = 1 !alter status of site to 'leaves on' - currentSite%dleafoffdate = t !record leaf on date + currentSite%dstatus = leaves_off !alter status of site to 'leaves off' + currentSite%dleafoffdate = t !record leaf on date endif endif !LEAF OFF: DROUGHT DECIDUOUS DRYNESS - if the soil gets too dry, and the leaves have already been on a while... - if (currentSite%dstatus == 2.and.t >= 10)then !D* + if (currentSite%dstatus == leaves_on.and.t >= 10)then !D* if (sum(currentSite%water_memory(1:10)/10._r8) <= ED_val_phen_drought_threshold)then if (timesincedleafon > 100)then !B* Have the leaves been on for some reasonable length of time? To prevent flickering. - currentSite%dstatus = 1 !alter status of site to 'leaves on' + currentSite%dstatus = leaves_off !alter status of site to 'leaves on' currentSite%dleafoffdate = t !record leaf on date endif endif @@ -672,10 +674,10 @@ subroutine phenology_leafonoff(currentSite) !COLD LEAF ON if (EDPftvarcon_inst%season_decid(ipft) == 1)then - if (currentSite%status == 2)then !we have just moved to leaves being on . - if (currentCohort%status_coh == 1)then !Are the leaves currently off? - currentCohort%status_coh = 2 ! Leaves are on, so change status to - ! stop flow of carbon out of bstore. + if (currentSite%status == leaves_on)then !we have just moved to leaves being on . + if (currentCohort%status_coh == leaves_off)then !Are the leaves currently off? + currentCohort%status_coh = leaves_on ! Leaves are on, so change status to + ! stop flow of carbon out of bstore. if(store_c>nearzero) then store_c_transfer_frac = & @@ -690,15 +692,23 @@ subroutine phenology_leafonoff(currentSite) currentCohort%laimemory = 0.0_r8 + ! Set the leaf age class to all be in the youngest since + ! phenology is pushing out new leaves (THIS ASSUMES THAT WE ARE + ! STARTING FROM HAVING NO LEAVES (ie a clean flushing of leaves) + + currentCohort%frac_leaf_aclass(:) = 0._r8 + currentCohort%frac_leaf_aclass(1) = 1._r8 + + endif !pft phenology endif ! growing season !COLD LEAF OFF - if (currentSite%status == 1)then !past leaf drop day? Leaves still on tree? - if (currentCohort%status_coh == 2)then ! leaves have not dropped + if (currentSite%status == leaves_off)then !past leaf drop day? Leaves still on tree? + if (currentCohort%status_coh == leaves_on)then ! leaves have not dropped ! This sets the cohort to the "leaves off" flag - currentCohort%status_coh = 1 + currentCohort%status_coh = leaves_off ! Remember what the lai was (leaf mass actually) was for next year ! the same amount back on in the spring... @@ -711,7 +721,8 @@ subroutine phenology_leafonoff(currentSite) call PRTDeciduousTurnover(currentCohort%prt,ipft, & leaf_organ, leaf_drop_fraction) - + + endif !leaf status endif !currentSite status endif !season_decid @@ -719,15 +730,15 @@ subroutine phenology_leafonoff(currentSite) !DROUGHT LEAF ON if (EDPftvarcon_inst%stress_decid(ipft) == 1)then - if (currentSite%dstatus == 2)then + if ( currentSite%dstatus == leaves_on )then ! we have just moved to leaves being on . - if (currentCohort%status_coh == 1)then + if (currentCohort%status_coh == leaves_off)then !is it the leaf-on day? Are the leaves currently off? - currentCohort%status_coh = 2 ! Leaves are on, so change status to - ! stop flow of carbon out of bstore. + currentCohort%status_coh = leaves_on ! Leaves are on, so change status to + ! stop flow of carbon out of bstore. if(store_c>nearzero) then store_c_transfer_frac = & @@ -743,15 +754,22 @@ subroutine phenology_leafonoff(currentSite) currentCohort%laimemory = 0.0_r8 + ! Set the leaf age class to all be in the youngest since + ! phenology is pushing out new leaves (THIS ASSUMES THAT WE ARE + ! STARTING FROM HAVING NO LEAVES (ie a clean flushing of leaves) + + currentCohort%frac_leaf_aclass(:) = 0._r8 + currentCohort%frac_leaf_aclass(1) = 1._r8 + endif !currentCohort status again? endif !currentSite status !DROUGHT LEAF OFF - if (currentSite%dstatus == 1)then - if (currentCohort%status_coh == 2)then ! leaves have not dropped + if (currentSite%dstatus == leaves_off)then + if (currentCohort%status_coh == leaves_on)then ! leaves have not dropped ! This sets the cohort to the "leaves off" flag - currentCohort%status_coh = 1 + currentCohort%status_coh = leaves_off ! Remember what the lai (leaf mass actually) was for next year currentCohort%laimemory = leaf_c @@ -933,10 +951,10 @@ subroutine seed_germination( currentSite, currentPatch ) currentPatch%seed_germination(p) = min(currentSite%seed_bank(p) * & EDPftvarcon_inst%germination_timescale(p),max_germination) !set the germination only under the growing season...c.xu - if (EDPftvarcon_inst%season_decid(p) == 1.and.currentSite%status == 1)then + if (EDPftvarcon_inst%season_decid(p) == 1.and.currentSite%status == leaves_off)then currentPatch%seed_germination(p) = 0.0_r8 endif - if (EDPftvarcon_inst%stress_decid(p) == 1.and.currentSite%dstatus == 1)then + if (EDPftvarcon_inst%stress_decid(p) == 1.and.currentSite%dstatus == leaves_off)then currentPatch%seed_germination(p) = 0.0_r8 endif enddo @@ -991,11 +1009,11 @@ subroutine recruitment( currentSite, currentPatch, bc_in ) call bstore_allom(temp_cohort%dbh,ft,temp_cohort%canopy_trim,b_store) temp_cohort%laimemory = 0.0_r8 - if (EDPftvarcon_inst%season_decid(temp_cohort%pft) == 1.and.currentSite%status == 1)then + if (EDPftvarcon_inst%season_decid(temp_cohort%pft) == 1.and.currentSite%status == leaves_off)then temp_cohort%laimemory = b_leaf b_leaf = 0.0_r8 endif - if (EDPftvarcon_inst%stress_decid(temp_cohort%pft) == 1.and.currentSite%dstatus == 1)then + if (EDPftvarcon_inst%stress_decid(temp_cohort%pft) == 1.and.currentSite%dstatus == leaves_off)then temp_cohort%laimemory = b_leaf b_leaf = 0.0_r8 endif @@ -1007,7 +1025,7 @@ subroutine recruitment( currentSite, currentPatch, bc_in ) if (EDPftvarcon_inst%evergreen(ft) == 1) then temp_cohort%laimemory = 0._r8 - cohortstatus = 2 + cohortstatus = leaves_on endif if (hlm_use_ed_prescribed_phys .eq. ifalse .or. EDPftvarcon_inst%prescribed_recruitment(ft) .lt. 0. ) then diff --git a/parteh/PRTLossFluxesMod.F90 b/parteh/PRTLossFluxesMod.F90 index 57d727206b..e2dffde6ad 100644 --- a/parteh/PRTLossFluxesMod.F90 +++ b/parteh/PRTLossFluxesMod.F90 @@ -524,7 +524,7 @@ end subroutine DeciduousTurnoverSimpleRetranslocation ! ==================================================================================== - subroutine PRTMaintTurnover(prt,ipft) + subroutine PRTMaintTurnover(prt,ipft,senescent_frac) ! --------------------------------------------------------------------------------- ! Generic subroutine (wrapper) calling specialized routines handling @@ -534,7 +534,7 @@ subroutine PRTMaintTurnover(prt,ipft) integer,intent(in) :: ipft if ( int(EDPftvarcon_inst%turnover_retrans_mode(ipft)) == 1 ) then - call MaintTurnoverSimpleRetranslocation(prt,ipft) + call MaintTurnoverSimpleRetranslocation(prt,ipft,senescent_frac) else write(fates_log(),*) 'A maintenance/retranslocation mode was specified' write(fates_log(),*) 'that is unknown.' @@ -548,7 +548,7 @@ end subroutine PRTMaintTurnover ! =================================================================================== - subroutine MaintTurnoverSimpleRetranslocation(prt,ipft) + subroutine MaintTurnoverSimpleRetranslocation(prt,ipft,senescent_frac) ! --------------------------------------------------------------------------------- ! This subroutine removes biomass from all applicable pools due to @@ -567,14 +567,17 @@ subroutine MaintTurnoverSimpleRetranslocation(prt,ipft) ! 4) There are currently no reaction costs associated with re-translocation ! --------------------------------------------------------------------------------- - class(prt_vartypes) :: prt - integer,intent(in) :: ipft + class(prt_vartypes) :: prt + integer, intent(in) :: ipft + real(r8), intent(in) :: senescent_frac integer :: i_var ! the variable index integer :: element_id ! the element associated w/ each variable integer :: organ_id ! the organ associated w/ each variable integer :: i_pos ! spatial position loop counter - + integer :: aclass_sen_id ! the index of the leaf age class dimension + ! associated with the senescing pool + real(r8) :: turnover ! Actual turnover removed from each ! pool [kg] real(r8) :: retrans ! A temp for the actual re-translocated mass @@ -612,11 +615,18 @@ subroutine MaintTurnoverSimpleRetranslocation(prt,ipft) base_turnover(fnrt_organ) = 0.0_r8 end if + + ! The last index of the leaf longevity array contains the turnover + ! timescale for the senescent pool + aclass_sen_id = size(EDPftvarcon_inst%leaf_long(ipft,:)) + ! Only EVERGREENS HAVE MAINTENANCE LEAF TURNOVER ! ------------------------------------------------------------------------------------- - if ( (EDPftvarcon_inst%leaf_long(ipft,0) > nearzero ) .and. & - (EDPftvarcon_inst%evergreen(ipft) == 1) ) then - base_turnover(leaf_organ) = hlm_freq_day / EDPftvarcon_inst%leaf_long(ipft,0) + if ( (EDPftvarcon_inst%leaf_long(ipft,aclass_sen_id) > nearzero ) .and. & + (EDPftvarcon_inst%evergreen(ipft) == itrue) ) then + + base_turnover(leaf_organ) = senescent_frac * hlm_freq_day / & + EDPftvarcon_inst%leaf_long(ipft,aclass_sen_id) else base_turnover(leaf_organ) = 0.0_r8 endif From a3219e99ca5ca29fd9ce31850473030629a35386 Mon Sep 17 00:00:00 2001 From: chonggang xu Date: Mon, 17 Dec 2018 16:17:45 -0800 Subject: [PATCH 09/34] Added leaf age to parteh framework --- parteh/PRTAllometricCarbonMod.F90 | 103 +++++++++++++++++++++++------- parteh/PRTGenericMod.F90 | 63 +++++++++++++----- parteh/PRTLossFluxesMod.F90 | 25 ++++++-- 3 files changed, 147 insertions(+), 44 deletions(-) diff --git a/parteh/PRTAllometricCarbonMod.F90 b/parteh/PRTAllometricCarbonMod.F90 index 8a2da57bc2..e0c81cc61d 100644 --- a/parteh/PRTAllometricCarbonMod.F90 +++ b/parteh/PRTAllometricCarbonMod.F90 @@ -23,6 +23,7 @@ module PRTAllometricCarbonMod use PRTGenericMod , only : repro_organ use PRTGenericMod , only : struct_organ use PRTGenericMod , only : un_initialized + use PRTGenericMod , only : prt_carbon_allom_hyp use FatesAllometryMod , only : bleaf use FatesAllometryMod , only : bsap_allom @@ -48,6 +49,7 @@ module PRTAllometricCarbonMod use FatesConstantsMod , only : itrue + implicit none private @@ -100,6 +102,9 @@ module PRTAllometricCarbonMod ! ------------------------------------------------------------------------------------- integer, parameter :: icd = 1 ! Only 1 coordinate per variable + + ! This is the maximum number of leaf age pools, used for allocating scratch space + integer, parameter :: max_nleafage = 4 ! ------------------------------------------------------------------------------------- ! This is the core type that holds this specific @@ -171,6 +176,8 @@ subroutine InitPRTGlobalAllometricCarbon() ! ! ----------------------------------------------------------------------------------- + integer :: nleafage + allocate(prt_global_ac) ! The "state descriptor" object holds things like the names, the symbols, the units @@ -181,13 +188,32 @@ subroutine InitPRTGlobalAllometricCarbon() prt_global_ac%hyp_name = 'Allometric Carbon Only' + prt_global_ac%hyp_id = prt_carbon_allom_hyp + ! Set mapping tables to zero call prt_global_ac%ZeroGlobal() + + ! The number of leaf age classes can be determined from the parameter file, + ! notably the size of the leaf-longevity parameter's second dimension. + ! This is the same value in FatesInterfaceMod.F90 + + nleafage = size(EDPftvarcon_inst%leaf_long,dim=2) + + if(nleafage>max_nleafage) then + write(fates_log(),*) 'The allometric carbon PARTEH hypothesis' + write(fates_log(),*) 'sets a maximum number of leaf age classes' + write(fates_log(),*) 'used for scratch space. The model wants' + write(fates_log(),*) 'exceed that. Simply increase max_nleafage' + write(fates_log(),*) 'found in parteh/PRTAllometricCarbonMod.F90' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + ! Register the variables. Each variable must be associated with a global identifier - ! for an organ and species. + ! for an organ and species. Leaves are a little special in that they are discretized + ! further by age class. Although the code works fine if this collapses to 1. - call prt_global_ac%RegisterVarInGlobal(leaf_c_id,"Leaf Carbon","leaf_c",leaf_organ,carbon12_element,icd) + call prt_global_ac%RegisterVarInGlobal(leaf_c_id,"Leaf Carbon","leaf_c",leaf_organ,carbon12_element,nleafage) call prt_global_ac%RegisterVarInGlobal(fnrt_c_id,"Fine Root Carbon","fnrt_c",fnrt_organ,carbon12_element,icd) call prt_global_ac%RegisterVarInGlobal(sapw_c_id,"Sapwood Carbon","sapw_c",sapw_organ,carbon12_element,icd) call prt_global_ac%RegisterVarInGlobal(store_c_id,"Storage Carbon","store_c",store_organ,carbon12_element,icd) @@ -301,7 +327,8 @@ subroutine DailyPRTAllometricCarbon(this) real(r8) :: repro_c_flux ! Transfer into reproduction at the final stage [kgC] real(r8) :: struct_c_flux ! Transfer into structure at various stages [kgC] - real(r8) :: leaf_c0 ! Initial value of carbon used to determine net flux + real(r8),dimension(max_nleafage) :: leaf_c0 + ! Initial value of carbon used to determine net flux real(r8) :: fnrt_c0 ! during this routine real(r8) :: sapw_c0 ! "" real(r8) :: store_c0 ! "" @@ -322,6 +349,8 @@ subroutine DailyPRTAllometricCarbon(this) real(r8) :: hite_out ! dummy height variable integer :: i_var ! local index for iterating state variables + integer :: nleafage ! number of leaf age classifications + real(r8) :: leaf_age_flux ! carbon mass flux between leaf age classification pools ! Integegrator variables c_pool is "mostly" carbon variables, it also includes @@ -333,15 +362,21 @@ subroutine DailyPRTAllometricCarbon(this) logical,dimension(n_integration_vars) :: c_mask ! Mask of active pools during integration integer , parameter :: max_substeps = 300 ! Maximum allowable iterations + real(r8), parameter :: max_trunc_error = 1.0_r8 ! Maximum allowable truncation error + integer, parameter :: ODESolve = 2 ! 1=RKF45, 2=Euler + integer, parameter :: iexp_leaf = 1 ! index 1 is the expanding (i.e. youngest) + ! leaf age class, and therefore + ! all new allocation goes into that pool + real(r8) :: intgr_params(num_bc_in) ! The boundary conditions to this routine, ! are pressed into an array that is also ! passed to the integrators associate( & - leaf_c => this%variables(leaf_c_id)%val(icd), & + leaf_c => this%variables(leaf_c_id)%val, & fnrt_c => this%variables(fnrt_c_id)%val(icd), & sapw_c => this%variables(sapw_c_id)%val(icd), & store_c => this%variables(store_c_id)%val(icd), & @@ -371,12 +406,33 @@ subroutine DailyPRTAllometricCarbon(this) ! transport flux "%net_alloc" at the end. ! ----------------------------------------------------------------------------------- - leaf_c0 = leaf_c ! Set initial leaf carbon - fnrt_c0 = fnrt_c ! Set initial fine-root carbon - sapw_c0 = sapw_c ! Set initial sapwood carbon - store_c0 = store_c ! Set initial storage carbon - repro_c0 = repro_c ! Set initial reproductive carbon - struct_c0 = struct_c ! Set initial structural carbon + nleafage = prt_global%state_descriptor(leaf_c_id)%num_pos + + + leaf_c0(1:nleafage) = leaf_c(1:nleafage) ! Set initial leaf carbon + fnrt_c0 = fnrt_c ! Set initial fine-root carbon + sapw_c0 = sapw_c ! Set initial sapwood carbon + store_c0 = store_c ! Set initial storage carbon + repro_c0 = repro_c ! Set initial reproductive carbon + struct_c0 = struct_c ! Set initial structural carbon + + + ! ----------------------------------------------------------------------------------- + ! If we have more than one leaf age classification, allow + ! some leaf biomass to transition to the older classes. NOTE! This is not handling + ! losses due to turnover (ie. flux from the oldest senescing class). This is only + ! internal. + ! (rgk 12-15-2018: Have Chonggang confirm that aging should not be restricted + ! to evergreens) + ! ----------------------------------------------------------------------------------- + + if(nleafage>1) then + do iage = 1,nleafage-1 + leaf_age_flux = leaf_c0(iage) * hlm_freq_day / EDPftvarcon_inst%leaf_long(ipft,iage) + leaf_c(iage) = leaf_c(iage) - leaf_age_flux + leaf_c(iage+1) = leaf_c(iage+1) + leaf_age_flux + end do + end if ! ----------------------------------------------------------------------------------- @@ -441,7 +497,7 @@ subroutine DailyPRTAllometricCarbon(this) if( EDPftvarcon_inst%evergreen(ipft) ==1 ) then leaf_c_demand = max(0.0_r8, & - EDPftvarcon_inst%leaf_stor_priority(ipft)*this%variables(leaf_c_id)%turnover(icd)) + EDPftvarcon_inst%leaf_stor_priority(ipft)*sum(this%variables(leaf_c_id)%turnover(:))) else leaf_c_demand = 0.0_r8 end if @@ -460,8 +516,9 @@ subroutine DailyPRTAllometricCarbon(this) max(0.0_r8,(store_c+carbon_balance)* & (leaf_c_demand/total_c_demand))) - carbon_balance = carbon_balance - leaf_c_flux - leaf_c = leaf_c + leaf_c_flux + ! Add carbon to the youngest age pool (i.e iexp_leaf = index 1) + carbon_balance = carbon_balance - leaf_c_flux + leaf_c(iexp_leaf) = leaf_c(iexp_leaf) + leaf_c_flux ! If we are testing b4b, then we pay this even if we don't have the carbon fnrt_c_flux = min(fnrt_c_demand, & @@ -503,7 +560,7 @@ subroutine DailyPRTAllometricCarbon(this) ! carbon balance is guaranteed to be >=0 beyond this point ! ----------------------------------------------------------------------------------- - leaf_c_demand = max(0.0_r8,(target_leaf_c - leaf_c)) + leaf_c_demand = max(0.0_r8,(target_leaf_c - sum(leaf_c(1:nleafage)))) fnrt_c_demand = max(0.0_r8,(target_fnrt_c - fnrt_c)) total_c_demand = leaf_c_demand + fnrt_c_demand @@ -513,7 +570,7 @@ subroutine DailyPRTAllometricCarbon(this) leaf_c_flux = min(leaf_c_demand, & carbon_balance*(leaf_c_demand/total_c_demand)) carbon_balance = carbon_balance - leaf_c_flux - leaf_c = leaf_c + leaf_c_flux + leaf_c(iexp_leaf) = leaf_c(iexp_leaf) + leaf_c_flux fnrt_c_flux = min(fnrt_c_demand, & carbon_balance*(fnrt_c_demand/total_c_demand)) @@ -530,7 +587,7 @@ subroutine DailyPRTAllometricCarbon(this) ! ----------------------------------------------------------------------------------- if( carbon_balance > nearzero ) then - leaf_below_target = max(target_leaf_c - leaf_c,0.0_r8) + leaf_below_target = max(target_leaf_c - sum(leaf_c(1:nleafage)),0.0_r8) fnrt_below_target = max(target_fnrt_c - fnrt_c,0.0_r8) sapw_below_target = max(target_sapw_c - sapw_c,0.0_r8) store_below_target = max(target_store_c - store_c,0.0_r8) @@ -553,10 +610,10 @@ subroutine DailyPRTAllometricCarbon(this) end if carbon_balance = carbon_balance - leaf_c_flux - leaf_c = leaf_c + leaf_c_flux + leaf_c(iexp_leaf) = leaf_c(iexp_leaf) + leaf_c_flux carbon_balance = carbon_balance - fnrt_c_flux - fnrt_c = fnrt_c + fnrt_c_flux + fnrt_c = fnrt_c + fnrt_c_flux carbon_balance = carbon_balance - sapw_c_flux sapw_c = sapw_c + sapw_c_flux @@ -622,7 +679,7 @@ subroutine DailyPRTAllometricCarbon(this) end if - call TargetAllometryCheck(leaf_c, fnrt_c, sapw_c, & + call TargetAllometryCheck(sum(leaf_c(1:nleafage)), fnrt_c, sapw_c, & store_c, struct_c, & target_leaf_c, target_fnrt_c, & target_sapw_c, target_store_c, target_struct_c, & @@ -657,7 +714,7 @@ subroutine DailyPRTAllometricCarbon(this) ! large and thus won't accumulate such large ! errors, we always mask as true. - c_pool(leaf_c_id) = leaf_c + c_pool(leaf_c_id) = sum(leaf_c(1:nleafage)) c_pool(fnrt_c_id) = fnrt_c c_pool(sapw_c_id) = sapw_c c_pool(store_c_id) = store_c @@ -732,7 +789,7 @@ subroutine DailyPRTAllometricCarbon(this) write(fates_log(),*) 'carbon_balance',carbon_balance write(fates_log(),*) 'deltaC',deltaC write(fates_log(),*) 'totalC',totalC - write(fates_log(),*) 'leaf:',grow_leaf,target_leaf_c,target_leaf_c - leaf_c + write(fates_log(),*) 'leaf:',grow_leaf,target_leaf_c,target_leaf_c - sum(leaf_c(:)) write(fates_log(),*) 'fnrt:',grow_fnrt,target_fnrt_c,target_fnrt_c - fnrt_c write(fates_log(),*) 'sap:',grow_sapw,target_sapw_c, target_sapw_c - sapw_c write(fates_log(),*) 'store:',grow_store,target_store_c,target_store_c - store_c @@ -751,7 +808,7 @@ subroutine DailyPRTAllometricCarbon(this) if( (totalC < calloc_abs_error) .and. (step_pass) )then ierr = 0 - leaf_c_flux = c_pool(leaf_c_id) - leaf_c + leaf_c_flux = c_pool(leaf_c_id) - sum(leaf_c(1:nleafage)) fnrt_c_flux = c_pool(fnrt_c_id) - fnrt_c sapw_c_flux = c_pool(sapw_c_id) - sapw_c store_c_flux = c_pool(store_c_id) - store_c @@ -771,7 +828,7 @@ subroutine DailyPRTAllometricCarbon(this) repro_c_flux = repro_c_flux*flux_adj carbon_balance = carbon_balance - leaf_c_flux - leaf_c = leaf_c + leaf_c_flux + leaf_c(iexp_leaf) = leaf_c(iexp_leaf) + leaf_c_flux carbon_balance = carbon_balance - fnrt_c_flux fnrt_c = fnrt_c + fnrt_c_flux diff --git a/parteh/PRTGenericMod.F90 b/parteh/PRTGenericMod.F90 index 945283694c..7f9c99dac0 100644 --- a/parteh/PRTGenericMod.F90 +++ b/parteh/PRTGenericMod.F90 @@ -320,6 +320,10 @@ module PRTGenericMod ! Note that index 0 is reserved for "all" or "irrelevant" character(len=maxlen_varname) :: hyp_name + ! This is the hypothesis index, used internally for some non-specific + ! routines where different methods are applied + integer :: hyp_id + ! This will save the specific variable id associated with each organ and element integer, dimension(0:num_organ_types,0:num_element_types) :: sp_organ_map @@ -777,26 +781,55 @@ subroutine WeightedFusePRTVartypes(this,donor_prt_obj, recipient_fuse_weight, po if(present(position_id)) then pos_id = position_id - else - pos_id = 1 - end if - do i_var = 1, prt_global%num_vars + do i_var = 1, prt_global%num_vars + + this%variables(i_var)%val(pos_id) = recipient_fuse_weight * this%variables(i_var)%val(pos_id) + & + (1.0_r8-recipient_fuse_weight) * donor_prt_obj%variables(i_var)%val(pos_id) + + this%variables(i_var)%val0(pos_id) = recipient_fuse_weight * this%variables(i_var)%val0(pos_id) + & + (1.0_r8-recipient_fuse_weight) * donor_prt_obj%variables(i_var)%val0(pos_id) + + this%variables(i_var)%net_alloc(pos_id) = recipient_fuse_weight * this%variables(i_var)%net_alloc(pos_id) + & + (1.0_r8-recipient_fuse_weight) * donor_prt_obj%variables(i_var)%net_alloc(pos_id) + + this%variables(i_var)%turnover(pos_id) = recipient_fuse_weight * this%variables(i_var)%turnover(pos_id) + & + (1.0_r8-recipient_fuse_weight) * donor_prt_obj%variables(i_var)%turnover(pos_id) + + this%variables(i_var)%burned(pos_id) = recipient_fuse_weight * this%variables(i_var)%burned(pos_id) + & + (1.0_r8-recipient_fuse_weight) * donor_prt_obj%variables(i_var)%burned(pos_id) - this%variables(i_var)%val(pos_id) = recipient_fuse_weight * this%variables(i_var)%val(pos_id) + & - (1.0_r8-recipient_fuse_weight) * donor_prt_obj%variables(i_var)%val(pos_id) + end do + + else - this%variables(i_var)%val0(pos_id) = recipient_fuse_weight * this%variables(i_var)%val0(pos_id) + & - (1.0_r8-recipient_fuse_weight) * donor_prt_obj%variables(i_var)%val0(pos_id) + ! If no argument regarding positions is supplied, then fuse all positions + ! sequentially + do i_var = 1, prt_global%num_vars + + do pos_id = 1,prt_global%state_descriptor(leaf_c_id)%num_pos + + this%variables(i_var)%val(pos_id) = recipient_fuse_weight * this%variables(i_var)%val(pos_id) + & + (1.0_r8-recipient_fuse_weight) * donor_prt_obj%variables(i_var)%val(pos_id) + + this%variables(i_var)%val0(pos_id) = recipient_fuse_weight * this%variables(i_var)%val0(pos_id) + & + (1.0_r8-recipient_fuse_weight) * donor_prt_obj%variables(i_var)%val0(pos_id) + + this%variables(i_var)%net_alloc(pos_id) = recipient_fuse_weight * this%variables(i_var)%net_alloc(pos_id) + & + (1.0_r8-recipient_fuse_weight) * donor_prt_obj%variables(i_var)%net_alloc(pos_id) + + this%variables(i_var)%turnover(pos_id) = recipient_fuse_weight * this%variables(i_var)%turnover(pos_id) + & + (1.0_r8-recipient_fuse_weight) * donor_prt_obj%variables(i_var)%turnover(pos_id) + + this%variables(i_var)%burned(pos_id) = recipient_fuse_weight * this%variables(i_var)%burned(pos_id) + & + (1.0_r8-recipient_fuse_weight) * donor_prt_obj%variables(i_var)%burned(pos_id) + + end do + end do - this%variables(i_var)%net_alloc(pos_id) = recipient_fuse_weight * this%variables(i_var)%net_alloc(pos_id) + & - (1.0_r8-recipient_fuse_weight) * donor_prt_obj%variables(i_var)%net_alloc(pos_id) - - this%variables(i_var)%turnover(pos_id) = recipient_fuse_weight * this%variables(i_var)%turnover(pos_id) + & - (1.0_r8-recipient_fuse_weight) * donor_prt_obj%variables(i_var)%turnover(pos_id) + end if - this%variables(i_var)%burned(pos_id) = recipient_fuse_weight * this%variables(i_var)%burned(pos_id) + & - (1.0_r8-recipient_fuse_weight) * donor_prt_obj%variables(i_var)%burned(pos_id) + end do diff --git a/parteh/PRTLossFluxesMod.F90 b/parteh/PRTLossFluxesMod.F90 index e2dffde6ad..f7d0c9b3c0 100644 --- a/parteh/PRTLossFluxesMod.F90 +++ b/parteh/PRTLossFluxesMod.F90 @@ -548,7 +548,7 @@ end subroutine PRTMaintTurnover ! =================================================================================== - subroutine MaintTurnoverSimpleRetranslocation(prt,ipft,senescent_frac) + subroutine MaintTurnoverSimpleRetranslocation(prt,ipft) ! --------------------------------------------------------------------------------- ! This subroutine removes biomass from all applicable pools due to @@ -569,7 +569,6 @@ subroutine MaintTurnoverSimpleRetranslocation(prt,ipft,senescent_frac) class(prt_vartypes) :: prt integer, intent(in) :: ipft - real(r8), intent(in) :: senescent_frac integer :: i_var ! the variable index integer :: element_id ! the element associated w/ each variable @@ -577,6 +576,11 @@ subroutine MaintTurnoverSimpleRetranslocation(prt,ipft,senescent_frac) integer :: i_pos ! spatial position loop counter integer :: aclass_sen_id ! the index of the leaf age class dimension ! associated with the senescing pool + integer :: ipos_1 ! the first index of the "position" + ! loop to cycle. For leaves, we only + ! generate maintenance fluxes from the last + ! senescing class; all other cases this + ! is assumed to be 1. real(r8) :: turnover ! Actual turnover removed from each ! pool [kg] @@ -617,15 +621,15 @@ subroutine MaintTurnoverSimpleRetranslocation(prt,ipft,senescent_frac) ! The last index of the leaf longevity array contains the turnover - ! timescale for the senescent pool + ! timescale for the senescent pool. aclass_sen_id = size(EDPftvarcon_inst%leaf_long(ipft,:)) - + ! Only EVERGREENS HAVE MAINTENANCE LEAF TURNOVER ! ------------------------------------------------------------------------------------- if ( (EDPftvarcon_inst%leaf_long(ipft,aclass_sen_id) > nearzero ) .and. & (EDPftvarcon_inst%evergreen(ipft) == itrue) ) then - base_turnover(leaf_organ) = senescent_frac * hlm_freq_day / & + base_turnover(leaf_organ) = hlm_freq_day / & EDPftvarcon_inst%leaf_long(ipft,aclass_sen_id) else base_turnover(leaf_organ) = 0.0_r8 @@ -670,7 +674,16 @@ subroutine MaintTurnoverSimpleRetranslocation(prt,ipft,senescent_frac) call endrun(msg=errMsg(__FILE__, __LINE__)) end if - do i_pos = 1, prt_global%state_descriptor(i_var)%num_pos + ! Hypotheses 1 & 2 assume that the leaf pools are statified by age + ! We only generate turnover from the last (senescing) position + if( (prt_global%hyp_id .le. 2) .and. & + (organ_id .eq. leaf_organ) ) then + ipos_1 = prt_global%state_descriptor(i_var)%num_pos + else + ipos_1 = 1 + end if + + do i_pos = ipos_1, prt_global%state_descriptor(i_var)%num_pos turnover = (1.0_r8 - retrans) * base_turnover(organ_id) * prt%variables(i_var)%val(i_pos) From 73431d5f84ace710c9200c779f5e6b1f6d630d5d Mon Sep 17 00:00:00 2001 From: chonggang xu Date: Tue, 18 Dec 2018 17:47:22 -0800 Subject: [PATCH 10/34] Finished first pass of getting leaf age partitions into parteh. Added cohort level biophysical rates, and removed leaf age fractions --- biogeochem/EDCanopyStructureMod.F90 | 8 +- biogeochem/EDCohortDynamicsMod.F90 | 142 +++++++++++++++++---- biogeochem/EDPatchDynamicsMod.F90 | 5 +- biogeochem/EDPhysiologyMod.F90 | 24 +--- biogeochem/FatesAllometryMod.F90 | 18 ++- biogeophys/FatesPlantRespPhotosynthMod.F90 | 53 ++++---- main/EDMainMod.F90 | 9 ++ main/EDPftvarcon.F90 | 90 ++++++++++++- main/EDTypesMod.F90 | 24 +++- main/FatesParameterDerivedMod.F90 | 1 - main/FatesRestartInterfaceMod.F90 | 1 + parameter_files/fates_params_default.cdl | 16 +-- parteh/PRTAllometricCarbonMod.F90 | 37 +++--- parteh/PRTGenericMod.F90 | 47 ++----- parteh/PRTLossFluxesMod.F90 | 18 ++- 15 files changed, 326 insertions(+), 167 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 50b034026d..5f91969ebf 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -504,7 +504,6 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) ! demoted to the understory allocate(copyc) - allocate(copyc%frac_leaf_aclass(nleafage)) call InitPRTCohort(copyc) if( hlm_use_planthydro.eq.itrue ) then call InitHydrCohort(currentSite,copyc) @@ -852,7 +851,6 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) elseif ( cc_gain > nearzero .and. cc_gain < currentCohort%c_area) then allocate(copyc) - allocate(copyc%frac_leaf_aclass(nleafage)) call InitPRTCohort(copyc) if( hlm_use_planthydro.eq.itrue ) then call InitHydrCohort(CurrentSite,copyc) @@ -1082,7 +1080,7 @@ subroutine canopy_summarization( nsites, sites, bc_in ) currentCohort%treelai = tree_lai(leaf_c, & currentCohort%pft, currentCohort%c_area, currentCohort%n, & - currentCohort%canopy_layer, currentPatch%canopy_layer_tlai ) + currentCohort%canopy_layer, currentPatch%canopy_layer_tlai,currentCohort%vcmax25top ) canopy_leaf_area = canopy_leaf_area + currentCohort%treelai *currentCohort%c_area @@ -1257,11 +1255,11 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) currentCohort%treelai = tree_lai(leaf_c, currentCohort%pft, currentCohort%c_area, & currentCohort%n, currentCohort%canopy_layer, & - currentPatch%canopy_layer_tlai ) + currentPatch%canopy_layer_tlai,currentCohort%vcmax25top ) currentCohort%treesai = tree_sai(currentCohort%pft, currentCohort%dbh, currentCohort%canopy_trim, & currentCohort%c_area, currentCohort%n, currentCohort%canopy_layer, & - currentPatch%canopy_layer_tlai, currentCohort%treelai ) + currentPatch%canopy_layer_tlai, currentCohort%treelai ,currentCohort%vcmax25top) currentCohort%lai = currentCohort%treelai *currentCohort%c_area/currentPatch%total_canopy_area currentCohort%sai = currentCohort%treesai *currentCohort%c_area/currentPatch%total_canopy_area diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 0129fce997..d1427388f6 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -11,6 +11,7 @@ module EDCohortDynamicsMod use FatesConstantsMod , only : r8 => fates_r8 use FatesConstantsMod , only : fates_unset_int use FatesConstantsMod , only : itrue,ifalse + use FatesConstantsMod , only : nearzero use FatesInterfaceMod , only : hlm_days_per_year use FatesInterfaceMod , only : nleafage use EDPftvarcon , only : EDPftvarcon_inst @@ -22,6 +23,9 @@ module EDCohortDynamicsMod use EDTypesMod , only : min_npm2, min_nppatch use EDTypesMod , only : min_n_safemath use EDTypesMod , only : nlevleaf + use EDTypesMod , only : equal_leaf_aclass + use EDTypesMod , only : first_leaf_aclass + use EDTypesMod , only : nan_leaf_aclass use FatesInterfaceMod , only : hlm_use_planthydro use FatesInterfaceMod , only : hlm_parteh_mode use FatesPlantHydraulicsMod, only : FuseCohortHydraulics @@ -79,6 +83,7 @@ module EDCohortDynamicsMod public :: copy_cohort public :: count_cohorts public :: InitPRTCohort + public :: UpdateCohortBioPhysRates logical, parameter :: debug = .false. ! local debug flag @@ -123,6 +128,7 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, hite, dbh, bleaf, bfine ! spread crowns are in horizontal space integer, intent(in) :: leaf_aclass_init ! how to initialized the leaf age class ! distribution + integer :: iage ! loop counter for leaf age classes type(bc_in_type), intent(in) :: bc_in ! External boundary conditions ! @@ -130,11 +136,11 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, hite, dbh, bleaf, bfine type(ed_cohort_type), pointer :: new_cohort ! Pointer to New Cohort structure. type(ed_cohort_type), pointer :: storesmallcohort type(ed_cohort_type), pointer :: storebigcohort - integer :: tnull,snull ! are the tallest and shortest cohorts allocate + real(r8) :: frac_leaf_aclass(max_nleafage) ! Fraction of leaves in each age-class + integer :: tnull,snull ! are the tallest and shortest cohorts allocate !---------------------------------------------------------------------- allocate(new_cohort) - allocate(new_cohort%frac_leaf_aclass(nleafage)) call nan_cohort(new_cohort) ! Make everything in the cohort not-a-number call zero_cohort(new_cohort) ! Zero things that need to be zeroed. @@ -164,12 +170,12 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, hite, dbh, bleaf, bfine ! of this cohort if(leaf_aclass_init .eq. equal_leaf_aclass) then - new_cohort%frac_leaf_aclass(1:nleafage) = 1._r8 / real(nleafage,r8) + frac_leaf_aclass(1:nleafage) = 1._r8 / real(nleafage,r8) elseif(leaf_aclass_init .eq. first_leaf_aclass) then - new_cohort%frac_leaf_aclass(1:nleafage) = 0._r8 - new_cohort%frac_leaf_aclass(1) = 1._r8 + frac_leaf_aclass(1:nleafage) = 0._r8 + frac_leaf_aclass(1) = 1._r8 elseif(leaf_aclass_init .eq. nan_leaf_aclass) then - new_cohort%frac_leaf_aclass(1:nleafage) = nan + frac_leaf_aclass(1:nleafage) = nan else write(fates_log(),*) 'An unknown leaf age distribution was' write(fates_log(),*) 'requested during create cohort' @@ -192,7 +198,10 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, hite, dbh, bleaf, bfine select case(hlm_parteh_mode) case (prt_carbon_allom_hyp) - call SetState(new_cohort%prt,leaf_organ, carbon12_element, bleaf) + do iage = 1,nleafage + call SetState(new_cohort%prt,leaf_organ, carbon12_element, & + bleaf*frac_leaf_aclass(iage),iage) + end do call SetState(new_cohort%prt,fnrt_organ, carbon12_element, bfineroot) call SetState(new_cohort%prt,sapw_organ, carbon12_element, bsap) call SetState(new_cohort%prt,store_organ, carbon12_element, bstore) @@ -208,6 +217,9 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, hite, dbh, bleaf, bfine call new_cohort%prt%CheckInitialConditions() + ! This sets things like vcmax25top, that depend on the + ! leaf age fractions + call UpdateCohortBioPhysRates(currentCohort) call sizetype_class_index(new_cohort%dbh,new_cohort%pft, & new_cohort%size_class,new_cohort%size_by_pft_class) @@ -413,8 +425,12 @@ subroutine nan_cohort(cc_p) currentCohort%c_area = nan ! areal extent of canopy (m2) currentCohort%treelai = nan ! lai of tree (total leaf area (m2) / canopy area (m2) currentCohort%treesai = nan ! stem area index of tree (total stem area (m2) / canopy area (m2) - currentCohort%frac_leaf_aclass(:)= nan ! leaf age classes - + + currentCohort%vcmax25top = nan + currentCohort%jmax25top = nan + currentCohort%tpu25top = nan + currentCohort%kp25top = nan + ! CARBON FLUXES currentCohort%gpp_acc_hold = nan ! GPP: kgC/indiv/year currentCohort%gpp_tstep = nan ! GPP: kgC/indiv/timestep @@ -711,7 +727,6 @@ subroutine terminate_cohorts( currentSite, currentPatch, level ) ! Deallocate the cohort's PRT structure call currentCohort%prt%DeallocatePRTVartypes() deallocate(currentCohort%prt) - deallocate(currentCohort%frac_leaf_aclass) deallocate(currentCohort) nullify(currentCohort) @@ -879,27 +894,29 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) currentCohort%canopy_layer_yesterday = (currentCohort%n*currentCohort%canopy_layer_yesterday + & nextc%n*nextc%canopy_layer_yesterday)/newn - ! Leaf age class fractions - ! First calculate the total amount of leaf in each to enable - ! a weighted average of the two + ! Leaf biophysical rates (use leaf mass weighting) ! ----------------------------------------------------------------- leaf_c_next = nextc%prt%GetState(leaf_organ, all_carbon_elements)*nextc%n leaf_c_curr = currentCohort%prt%GetState(leaf_organ, all_carbon_elements)*currentCohort%n if( (leaf_c_next + leaf_c_curr) > nearzero ) then - currentCohort%frac_leaf_aclass(1:nleafage) = & - (leaf_c_curr*currentCohort%frac_leaf_aclass(1:nleafage) + & - leaf_c_next*nextc%frac_leaf_aclass(1:nleafage))/(leaf_c_next+leaf_c_curr) - - ! Force leaf age class to sum to unity - currentCohort%frac_leaf_aclass(1:nleafage) = & - currentCohort%frac_leaf_aclass(1:nleafage) / & - sum(currentCohort%frac_leaf_aclass(1:nleafage)) + + currentCohort%vcmax25top = (leaf_c_curr*currentCohort%vcmax25top + & + leaf_c_next*nextc%vcmax25top)/(leaf_c_next+leaf_c_curr) + currentCohort%jmax25top = (leaf_c_curr*currentCohort%jmax25top + & + leaf_c_next*nextc%jmax25top)/(leaf_c_next+leaf_c_curr) + currentCohort%tpu25top = (leaf_c_curr*currentCohort%tpu25top + & + leaf_c_next*nextc%tpu25top)/(leaf_c_next+leaf_c_curr) + currentCohort%kp25top = (leaf_c_curr*currentCohort%kp25top + & + leaf_c_next*nextc%kp25top)/(leaf_c_next+leaf_c_curr) + else - ! If there is no leaf at all... then we set - ! the age fraction to be all in the growing class - currentCohort%frac_leaf_aclass(1:nleafage) = 0._r8 - currentCohort%frac_leaf_aclass(1) = 1._r8 + + currentCohort%vcmax25top = 0._r8 + currentCohort%jmax25top = 0._r8 + currentCohort%tpu25top = 0._r8 + currentCohort%kp25top = 0._r8 + end if @@ -1032,7 +1049,6 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) ! Deallocate the cohort's PRT structure call nextc%prt%DeallocatePRTVartypes() deallocate(nextc%prt) - deallocate(nextc%frac_leaf_aclass) deallocate(nextc) nullify(nextc) @@ -1314,8 +1330,11 @@ subroutine copy_cohort( currentCohort,copyc ) ! This transfers the PRT objects over. call n%prt%CopyPRTVartypes(o%prt) - ! Leaf age class fractions - n%frac_leaf_aclass(1:nleafage) = o%frac_leaf_aclass(1:nleafage) + ! Leaf biophysical rates + n%vcmax25top = o%vcmax25top + n%jmax25top = o%jmax25top + n%tpu25top = o%tpu25top + n%kp25top = o%kp25top ! CARBON FLUXES n%gpp_acc_hold = o%gpp_acc_hold @@ -1434,6 +1453,73 @@ function count_cohorts( currentPatch ) result ( backcount ) endif end function count_cohorts + + ! =================================================================================== + + subroutine UpdateCohortBioPhysRates(currentCohort) + + ! -------------------------------------------------------------------------------- + ! This routine updates the four key biophysical rates of leaves + ! based on the changes in a cohort's leaf age proportions + ! + ! This should be called after growth. Growth occurs + ! after turnover and damage states are applied to the tree. + ! Therefore, following growth, the leaf mass fractions + ! of different age classes are unchanged until the next day. + ! -------------------------------------------------------------------------------- + + type(ed_cohort_type),intent(inout) :: currentCohort + + + real(r8) :: frac_leaf_aclass(max_nleafage) ! Fraction of leaves in each age-class + integer :: iage ! loop index for leaf ages + integer :: ipft ! plant functional type index + + ! First, calculate the fraction of leaves in each age class + ! It is assumed that each class has the same proportion + ! across leaf layers + + do iage = 1, nleafage + frac_leaf_aclass(iage) = & + currentCohort%prt%GetState(leaf_organ, all_carbon_elements,iage) + end do + + ! If there are leaves, then perform proportional weighting on the four rates + ! We assume that leaf age does not effect the specific leaf area, so the mass + ! fractions are applicable to these rates + + if(sum(frac_leaf_aclass(1:nleafage))>nearzero) then + + ipft = currentCohort%pft + + frac_leaf_aclass(1:nleafage) = frac_leaf_aclass(1:nleafage) / & + sum(frac_leaf_aclass(1:nleafage)) + + currentCohort%vcmax25top = sum(EDPftvarcon_inst%vcmax25top(ipft,1:nleafage) * & + frac_leaf_aclass(1:nleafage)) + + currentCohort%jmax25top = sum(param_derived%jmax25top(ipft,1:nleafage) * & + frac_leaf_aclass(1:nleafage)) + + currentCohort%tpu25top = sum(param_derived%tpu25top(ipft,1:nleafage) * & + frac_leaf_aclass(1:nleafage)) + + currentCohort%kp25top = sum(param_derived%kp25top(ipft,1:nleafage) * & + frac_leaf_aclass(1:nleafage)) + + else + + currentCohort%vcmax25top = 0._r8 + currentCohort%jmax25top = 0._r8 + currentCohort%tpu25top = 0._r8 + currentCohort%kp25top = 0._r8 + + end if + + + return + end subroutine UpdateCohortBioPhysRates + ! ============================================================================ diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 78abbef13c..458a949daf 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -401,8 +401,7 @@ subroutine spawn_patches( currentSite, bc_in) currentCohort => currentPatch%shortest do while(associated(currentCohort)) - allocate(nc) - allocate(nc%frac_leaf_aclass(nleafage)) + allocate(nc) if(hlm_use_planthydro.eq.itrue) call InitHydrCohort(CurrentSite,nc) call InitPRTCohort(nc) call zero_cohort(nc) @@ -676,7 +675,6 @@ subroutine spawn_patches( currentSite, bc_in) if(hlm_use_planthydro.eq.itrue) call DeallocateHydrCohort(nc) call nc%prt%DeallocatePRTVartypes() deallocate(nc%prt) - deallocate(nc%frac_leaf_aclass) deallocate(nc) endif @@ -1929,7 +1927,6 @@ subroutine dealloc_patch(cpatch) if(hlm_use_planthydro.eq.itrue) call DeallocateHydrCohort(ccohort) call ccohort%prt%DeallocatePRTVartypes() deallocate(ccohort%prt) - deallocate(ccohort%frac_leaf_aclass) deallocate(ccohort) ccohort => ncohort diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index a07d324408..a1438c406b 100644 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -90,7 +90,7 @@ module EDPhysiologyMod private :: seed_germination public :: flux_into_litter_pools public :: ZeroAllocationRates - + logical, parameter :: debug = .false. ! local debug flag character(len=*), parameter, private :: sourcefile = & @@ -251,11 +251,11 @@ subroutine trim_canopy( currentSite ) currentCohort%treelai = tree_lai(leaf_c, currentCohort%pft, currentCohort%c_area, & currentCohort%n, currentCohort%canopy_layer, & - currentPatch%canopy_layer_tlai ) + currentPatch%canopy_layer_tlai,currentCohort%vcmax25top ) currentCohort%treesai = tree_sai(currentCohort%pft, currentCohort%dbh, currentCohort%canopy_trim, & currentCohort%c_area, currentCohort%n, currentCohort%canopy_layer, & - currentPatch%canopy_layer_tlai, currentCohort%treelai ) + currentPatch%canopy_layer_tlai, currentCohort%treelai,currentCohort%vcmax25top ) currentCohort%nv = ceiling((currentCohort%treelai+currentCohort%treesai)/dinc_ed) @@ -692,14 +692,6 @@ subroutine phenology_leafonoff(currentSite) currentCohort%laimemory = 0.0_r8 - ! Set the leaf age class to all be in the youngest since - ! phenology is pushing out new leaves (THIS ASSUMES THAT WE ARE - ! STARTING FROM HAVING NO LEAVES (ie a clean flushing of leaves) - - currentCohort%frac_leaf_aclass(:) = 0._r8 - currentCohort%frac_leaf_aclass(1) = 1._r8 - - endif !pft phenology endif ! growing season @@ -754,13 +746,6 @@ subroutine phenology_leafonoff(currentSite) currentCohort%laimemory = 0.0_r8 - ! Set the leaf age class to all be in the youngest since - ! phenology is pushing out new leaves (THIS ASSUMES THAT WE ARE - ! STARTING FROM HAVING NO LEAVES (ie a clean flushing of leaves) - - currentCohort%frac_leaf_aclass(:) = 0._r8 - currentCohort%frac_leaf_aclass(1) = 1._r8 - endif !currentCohort status again? endif !currentSite status @@ -1814,9 +1799,6 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) end subroutine flux_into_litter_pools - ! =================================================================================== - - end module EDPhysiologyMod diff --git a/biogeochem/FatesAllometryMod.F90 b/biogeochem/FatesAllometryMod.F90 index 3c03524f3f..920105b8ee 100644 --- a/biogeochem/FatesAllometryMod.F90 +++ b/biogeochem/FatesAllometryMod.F90 @@ -549,7 +549,7 @@ end subroutine storage_fraction_of_target ! ===================================================================================== - real(r8) function tree_lai( leaf_c, pft, c_area, nplant, cl, canopy_lai) + real(r8) function tree_lai( leaf_c, pft, c_area, nplant, cl, canopy_lai, vcmax25top) ! ----------------------------------------------------------------------------------- ! LAI of individual trees is a function of the total leaf area and the total @@ -564,6 +564,8 @@ real(r8) function tree_lai( leaf_c, pft, c_area, nplant, cl, canopy_lai) integer, intent(in) :: cl ! canopy layer index real(r8), intent(in) :: canopy_lai(nclmax) ! total leaf area index of ! each canopy layer + real(r8), intent(in) :: vcmax25top ! maximum carboxylation rate at canopy + ! top, ref 25C ! !LOCAL VARIABLES: real(r8) :: leafc_per_unitarea ! KgC of leaf per m2 area of ground. @@ -600,7 +602,7 @@ real(r8) function tree_lai( leaf_c, pft, c_area, nplant, cl, canopy_lai) end if ! Coefficient for exponential decay of 1/sla with canopy depth: - kn = decay_coeff_kn(pft) + kn = decay_coeff_kn(pft,vcmax25top) ! take PFT-level maximum SLA value, even if under a thick canopy (which has units of m2/gC), ! and put into units of m2/kgC @@ -676,7 +678,7 @@ end function tree_lai ! ============================================================================ - real(r8) function tree_sai( pft, dbh, canopy_trim, c_area, nplant, cl, canopy_lai, treelai ) + real(r8) function tree_sai( pft, dbh, canopy_trim, c_area, nplant, cl, canopy_lai, treelai, vcmax25top ) ! ============================================================================ ! SAI of individual trees is a function of the LAI of individual trees @@ -691,13 +693,15 @@ real(r8) function tree_sai( pft, dbh, canopy_trim, c_area, nplant, cl, canopy_la real(r8), intent(in) :: canopy_lai(nclmax) ! total leaf area index of ! each canopy layer real(r8), intent(in) :: treelai ! tree LAI for checking purposes only + real(r8), intent(in) :: vcmax25top ! maximum carboxylation rate at top of crown + real(r8) :: target_bleaf real(r8) :: target_lai call bleaf(dbh,pft,canopy_trim,target_bleaf) - target_lai = tree_lai( target_bleaf, pft, c_area, nplant, cl, canopy_lai) + target_lai = tree_lai( target_bleaf, pft, c_area, nplant, cl, canopy_lai, vcmax25top) tree_sai = EDPftvarcon_inst%allom_sai_scaler(pft) * target_lai @@ -2081,7 +2085,7 @@ end subroutine jackson_beta_root_profile ! ===================================================================================== - real(r8) function decay_coeff_kn(pft) + real(r8) function decay_coeff_kn(pft,vcmax25top) ! --------------------------------------------------------------------------------- ! This function estimates the decay coefficient used to estimate vertical @@ -2095,6 +2099,8 @@ real(r8) function decay_coeff_kn(pft) !ARGUMENTS integer, intent(in) :: pft + real(r8),intent(in) :: vcmax25top + !LOCAL VARIABLES ! ----------------------------------------------------------------------------------- @@ -2103,7 +2109,7 @@ real(r8) function decay_coeff_kn(pft) ! kn = 0.11. Here, we derive kn from vcmax25 as in Lloyd et al ! (2010) Biogeosciences, 7, 1833-1859 - decay_coeff_kn = exp(0.00963_r8 * EDPftvarcon_inst%vcmax25top(pft,0) - 2.43_r8) + decay_coeff_kn = exp(0.00963_r8 * vcmax25top - 2.43_r8) return end function decay_coeff_kn diff --git a/biogeophys/FatesPlantRespPhotosynthMod.F90 b/biogeophys/FatesPlantRespPhotosynthMod.F90 index c9b022ee30..30d7b0b459 100644 --- a/biogeophys/FatesPlantRespPhotosynthMod.F90 +++ b/biogeophys/FatesPlantRespPhotosynthMod.F90 @@ -24,12 +24,15 @@ module FATESPlantRespPhotosynthMod use FatesGlobals, only : fates_log use FatesConstantsMod, only : r8 => fates_r8 use FatesConstantsMod, only : itrue + use FatesConstantsMod, only : nearzero use FatesInterfaceMod, only : hlm_use_planthydro use FatesInterfaceMod, only : hlm_parteh_mode use FatesInterfaceMod, only : numpft + use FatesInterfaceMod, only : nleafage use EDTypesMod, only : maxpft use EDTypesMod, only : nlevleaf use EDTypesMod, only : nclmax + use EDTypesMod, only : max_nleafage use EDTypesMod, only : do_fates_salinity use PRTGenericMod, only : prt_carbon_allom_hyp use PRTGenericMod, only : prt_cnp_flex_allom_hyp @@ -156,7 +159,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) real(r8) :: co2_cpoint ! CO2 compensation point (Pa) real(r8) :: btran_eff ! effective transpiration wetness factor (0 to 1) real(r8) :: bbb ! Ball-Berry minimum leaf conductance (umol H2O/m**2/s) - real(r8) :: kn(maxpft) ! leaf nitrogen decay coefficient + real(r8) :: kn ! leaf nitrogen decay coefficient real(r8) :: cf ! s m**2/umol -> s/m (ideal gas conversion) [umol/m3] real(r8) :: gb_mol ! leaf boundary layer conductance (molar form: [umol /m**2/s]) real(r8) :: ceair ! vapor pressure of air, constrained (Pa) @@ -202,8 +205,6 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) real(r8) :: lai_current ! the LAI in the current leaf layer real(r8) :: cumulative_lai ! the cumulative LAI, top down, to the leaf layer of interest - - ! ----------------------------------------------------------------------------------- ! Keeping these two definitions in case they need to be added later ! @@ -214,6 +215,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) integer :: cl,s,iv,j,ps,ft,ifp ! indices integer :: nv ! number of leaf layers integer :: NCL_p ! number of canopy layers in patch + integer :: iage ! loop counter for leaf age classes ! Parameters ! ----------------------------------------------------------------------- @@ -313,11 +315,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) do ft = 1,numpft - ! Bonan et al (2011) JGR, 116, doi:10.1029/2010JG001593 used - ! kn = 0.11. Here, derive kn from vcmax25 as in Lloyd et al - ! (2010) Biogeosciences, 7, 1833-1859 - kn(ft) = decay_coeff_kn(ft) ! This is probably unnecessary and already calculated @@ -373,7 +371,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) ! are there any leaves of this pft in this layer? if(currentPatch%canopy_mask(cl,ft) == 1)then - + ! Loop over leaf-layers do iv = 1,currentCohort%nv @@ -430,10 +428,16 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) if(do_fates_salinity)then btran_eff = btran_eff*currentPatch%bstress_sal_ft(ft) endif - + + ! Bonan et al (2011) JGR, 116, doi:10.1029/2010JG001593 used + ! kn = 0.11. Here, derive kn from vcmax25 as in Lloyd et al + ! (2010) Biogeosciences, 7, 1833-1859 + + kn = decay_coeff_kn(ft,currentCohort%vcmax25top) + ! Scale for leaf nitrogen profile - nscaler = exp(-kn(ft) * cumulative_lai) + nscaler = exp(-kn * cumulative_lai) ! Leaf maintenance respiration to match the base rate used in CN ! but with the new temperature functions for C3 and C4 plants. @@ -476,14 +480,15 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) ! These rates are the specific rates used in the actual photosynthesis ! calculations that take localized environmental effects (temperature) ! into consideration. + + call LeafLayerBiophysicalRates(currentPatch%ed_parsun_z(cl,ft,iv), & ! in ft, & ! in - currentCohort%frac_leaf_aclass(:), & ! in - EDPftvarcon_inst%vcmax25top(ft,:), & ! in - param_derived%jmax25top(ft,:), & ! in - param_derived%tpu25top(ft,:), & ! in - param_derived%kp25top(ft,:), & ! in + currentCohort%vcmax25top, & ! in + currentCohort%jmax25top, & ! in + currentCohort%tpu25top, & ! in + currentCohort%kp25top, & ! in nscaler, & ! in bc_in(s)%t_veg_pa(ifp), & ! in btran_eff, & ! in @@ -1705,7 +1710,6 @@ end subroutine LeafLayerMaintenanceRespiration subroutine LeafLayerBiophysicalRates( parsun_lsl, & ft, & - frac_leaf_aclass, & vcmax25top_ft, & jmax25top_ft, & tpu25top_ft, & @@ -1741,14 +1745,13 @@ subroutine LeafLayerBiophysicalRates( parsun_lsl, & real(r8), intent(in) :: parsun_lsl ! PAR absorbed in sunlit leaves for this layer integer, intent(in) :: ft ! (plant) Functional Type Index real(r8), intent(in) :: nscaler ! Scale for leaf nitrogen profile - real(r8), intent(in) :: frac_leaf_aclass(nleafage) ! Fraction of leaves in each age-class - real(r8), intent(in) :: vcmax25top_ft(nleafage) ! canopy top maximum rate of carboxylation at 25C + real(r8), intent(in) :: vcmax25top_ft ! canopy top maximum rate of carboxylation at 25C ! for this pft (umol CO2/m**2/s) - real(r8), intent(in) :: jmax25top_ft(nleafage) ! canopy top maximum electron transport rate at 25C + real(r8), intent(in) :: jmax25top_ft ! canopy top maximum electron transport rate at 25C ! for this pft (umol electrons/m**2/s) - real(r8), intent(in) :: tpu25top_ft(nleafage) ! canopy top triose phosphate utilization rate at 25C + real(r8), intent(in) :: tpu25top_ft ! canopy top triose phosphate utilization rate at 25C ! for this pft (umol CO2/m**2/s) - real(r8), intent(in) :: co2_rcurve_islope25top_ft ! initial slope of CO2 response curve + real(r8), intent(in) :: co2_rcurve_islope25top_ft ! initial slope of CO2 response curve ! (C4 plants) at 25C, canopy top, this pft real(r8), intent(in) :: veg_tempk ! vegetation temperature real(r8), intent(in) :: btran ! transpiration wetness factor (0 to 1) @@ -1809,9 +1812,11 @@ subroutine LeafLayerBiophysicalRates( parsun_lsl, & tpu = 0._r8 co2_rcurve_islope = 0._r8 else ! day time - vcmax25 = sum(vcmax25top_ft(1:nleafage) * frac_leaf_aclass(1:nleafage)) * nscaler - jmax25 = sum(jmax25top_ft(1:nleafage) * frac_leaf_aclass(1:nleafage)) * nscaler - tpu25 = sum(tpu25top_ft(1:nleafage) * frac_leaf_aclass(1:nleafage)) * nscaler + + ! Vcmax25top was already calculated to derive the nscaler function + vcmax25 = vcmax25top_ft * nscaler + jmax25 = jmax25top_ft * nscaler + tpu25 = tpu25top_ft * nscaler co2_rcurve_islope25 = co2_rcurve_islope25top_ft * nscaler ! Adjust for temperature diff --git a/main/EDMainMod.F90 b/main/EDMainMod.F90 index aca5a3df14..353ff0f339 100644 --- a/main/EDMainMod.F90 +++ b/main/EDMainMod.F90 @@ -33,6 +33,7 @@ module EDMainMod use EDPhysiologyMod , only : recruitment use EDPhysiologyMod , only : trim_canopy use EDPhysiologyMod , only : ZeroAllocationRates + use EDCohortDynamicsMod , only : UpdateCohortBioPhysRates use SFMainMod , only : fire_model use FatesSizeAgeTypeIndicesMod, only : get_age_class_index use EDtypesMod , only : ncwd @@ -337,12 +338,20 @@ subroutine ed_integrate_state_variables(currentSite, bc_in ) currentSite%flux_in = currentSite%flux_in + currentCohort%npp_acc * currentCohort%n + ! Conducte Maintenance Turnover (parteh) call currentCohort%prt%CheckMassConservation(ft,3) call PRTMaintTurnover(currentCohort%prt,ft) call currentCohort%prt%CheckMassConservation(ft,4) + + ! Conduct Growth (parteh) call currentCohort%prt%DailyPRT() call currentCohort%prt%CheckMassConservation(ft,5) + ! Update the leaf biophysical rates based on leaf ages + ! leaf ages should be fully updated for each cohort at this point + call UpdateCohortBioPhysRates(currentCohort) + + ! Transfer all reproductive tissues into seed production call PRTReproRelease(currentCohort%prt,repro_organ,carbon12_element, & 1.0_r8, currentCohort%seed_prod) diff --git a/main/EDPftvarcon.F90 b/main/EDPftvarcon.F90 index 434a28bbf1..1825122ad3 100644 --- a/main/EDPftvarcon.F90 +++ b/main/EDPftvarcon.F90 @@ -9,6 +9,8 @@ module EDPftvarcon use EDTypesMod , only : maxSWb, ivis, inir use FatesConstantsMod, only : r8 => fates_r8 use FatesConstantsMod, only : nearzero + use FatesConstantsMod, only : itrue, ifalse + use FatesConstantsMod, only : years_per_day use FatesGlobals, only : fates_log use FatesGlobals, only : endrun => fates_endrun @@ -1852,6 +1854,8 @@ subroutine FatesCheckParams(is_master, parteh_mode) integer :: npft ! number of PFTs integer :: ipft ! pft index + integer :: nleafage ! size of the leaf age class array + integer :: iage ! leaf age class index integer :: norgans ! size of the plant organ dimension npft = size(EDPftvarcon_inst%pft_used,1) @@ -2031,7 +2035,7 @@ subroutine FatesCheckParams(is_master, parteh_mode) write(fates_log(),*) ' Aborting' call endrun(msg=errMsg(sourcefile, __LINE__)) - end if + end if ! Check if photosynthetic pathway is neither C3/C4 ! ---------------------------------------------------------------------------------- @@ -2273,9 +2277,89 @@ subroutine FatesCheckParams(is_master, parteh_mode) end if end if + + ! Check turnover time-scales + + nleafage = size(EDPftvarcon_inst%leaf_long,dim=2) + + do iage = 1, nleafage + + if ( EDPftvarcon_inst%leaf_long(ipft,iage)>nearzero ) then + + ! Check that leaf turnover doesn't exeed 1 day + if ( (years_per_day / EDPftvarcon_inst%leaf_long(ipft,iage)) > 1._r8 ) then + write(fates_log(),*) 'Leaf turnover time-scale is greater than 1 day!' + write(fates_log(),*) 'ipft: ',ipft,' iage: ',iage + write(fates_log(),*) 'leaf_long(ipft,iage): ',EDPftvarcon_inst%leaf_long(ipft,iage),' [years]' + write(fates_log(),*) 'Aborting' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + ! Check to make sure that all other age-classes for this PFT also + ! have non-zero entries, it wouldn't make sense otherwise + if ( any(EDPftvarcon_inst%leaf_long(ipft,:) <= nearzero) ) then + write(fates_log(),*) 'You specified a leaf_long that is zero or' + write(fates_log(),*) 'invalid for a particular age class.' + write(fates_log(),*) 'Yet, other age classes for this PFT are non-zero.' + write(fates_log(),*) 'this doesnt make sense.' + write(fates_log(),*) 'ipft = ',ipft + write(fates_log(),*) 'leaf_long(ipft,:) = ',EDPftvarcon_inst%leaf_long(ipft,:) + write(fates_log(),*) 'Aborting' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + else + if (EDPftvarcon_inst%evergreen(ipft) .eq. itrue) then + write(fates_log(),*) 'You specified zero leaf turnover: ' + write(fates_log(),*) 'ipft: ',ipft,' iage: ',iage + write(fates_log(),*) 'leaf_long(ipft,iage): ',EDPftvarcon_inst%leaf_long(ipft,iage) + write(fates_log(),*) 'yet this is an evergreen PFT, and it only makes sense' + write(fates_log(),*) 'that an evergreen would have leaf maintenance turnover' + write(fates_log(),*) 'disable this error if you are ok with this' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + end if + end do + + if ( EDPftvarcon_inst%root_long(ipft)>nearzero ) then + + ! Check that root turnover doesn't exeed 1 day + if ( (years_per_day / EDPftvarcon_inst%root_long(ipft)) > 1._r8 ) then + write(fates_log(),*) 'Root turnover time-scale is greater than 1 day!' + write(fates_log(),*) 'ipft: ',ipft + write(fates_log(),*) 'root_long(ipft): ',EDPftvarcon_inst%root_long(ipft),' [years]' + write(fates_log(),*) 'Aborting' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + else + if (EDPftvarcon_inst%evergreen(ipft) .eq. itrue) then + write(fates_log(),*) 'You specified zero root turnover: ' + write(fates_log(),*) 'ipft: ',ipft + write(fates_log(),*) 'root_long(ipft): ',EDPftvarcon_inst%root_long(ipft) + write(fates_log(),*) 'yet this is an evergreen PFT, and it only makes sense' + write(fates_log(),*) 'that an evergreen would have root maintenance turnover' + write(fates_log(),*) 'disable this error if you are ok with this' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + end if + + ! Check Branch turnover doesn't exceed one day + if ( EDPftvarcon_inst%branch_turnover(ipft)>nearzero ) then + + ! Check that branch turnover doesn't exeed 1 day + if ( (years_per_day / EDPftvarcon_inst%branch_turnover(ipft)) > 1._r8 ) then + write(fates_log(),*) 'Branch turnover time-scale is greater than 1 day!' + write(fates_log(),*) 'ipft: ',ipft + write(fates_log(),*) 'branch_turnover(ipft): ',EDPftvarcon_inst%branch_turnover(ipft),' [years]' + write(fates_log(),*) 'Aborting' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + end if + + end do - - + return end subroutine FatesCheckParams diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 497c872d8a..a75b86fa45 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -27,6 +27,9 @@ module EDTypesMod ! the parameter file may determine that fewer ! are used, but this helps allocate scratch ! space and output arrays. + + integer, parameter :: max_nleafage = 4 ! This is the maximum number of leaf age pools, + ! used for allocating scratch space ! ------------------------------------------------------------------------------------- ! Radiation parameters @@ -247,11 +250,21 @@ module EDTypesMod real(r8) :: resp_acc real(r8) :: resp_acc_hold - real(r8),allocatable :: frac_leaf_aclass(:) ! This array's size is defined by the - ! number of leaf age classes. - ! and this defines the fraction of - ! leaf mass in different age classes - ! This should ALWAYS sum to 1 !! + ! The following four biophysical rates are assumed to be + ! at the canopy top, at reference temp 25C, and based on the + ! leaf age weighted average of the PFT parameterized values. The last + ! condition is why it is dynamic and tied to the cohort + + real(r8) :: vcmax25top ! Maximum carboxylation at the cohort's top + ! at reference temperature (25C). + real(r8) :: jmax25top ! canopy top: maximum electron transport + ! rate at 25C (umol electrons/m**2/s) + real(r8) :: tpu25top ! canopy top: triose phosphate utilization + ! rate at 25C (umol CO2/m**2/s) + real(r8) :: kp25top ! canopy top: initial slope of CO2 response + ! curve (C4 plants) at 25C + + real(r8) :: ts_net_uptake(nlevleaf) ! Net uptake of leaf layers: kgC/m2/timestep real(r8) :: year_net_uptake(nlevleaf) ! Net uptake of leaf layers: kgC/m2/year @@ -813,7 +826,6 @@ subroutine dump_cohort(ccohort) write(fates_log(),*) 'co%resp_acc = ', ccohort%resp_acc write(fates_log(),*) 'co%resp_acc_hold = ', ccohort%resp_acc_hold write(fates_log(),*) 'co%rdark = ', ccohort%rdark - write(fates_log(),*) 'co%frac_leaf_aclass = ', ccohort%frac_leaf_aclass(:) write(fates_log(),*) 'co%resp_m = ', ccohort%resp_m write(fates_log(),*) 'co%resp_g = ', ccohort%resp_g write(fates_log(),*) 'co%livestem_mr = ', ccohort%livestem_mr diff --git a/main/FatesParameterDerivedMod.F90 b/main/FatesParameterDerivedMod.F90 index f175114acc..15fc287b0c 100644 --- a/main/FatesParameterDerivedMod.F90 +++ b/main/FatesParameterDerivedMod.F90 @@ -89,6 +89,5 @@ subroutine Init(this,numpft) end associate return end subroutine Init - end module FatesParameterDerivedMod diff --git a/main/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index afdb49d65d..30e0d76d73 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -14,6 +14,7 @@ module FatesRestartInterfaceMod use FatesRestartVariableMod, only : fates_restart_variable_type use FatesInterfaceMod, only : bc_in_type use FatesInterfaceMod, only : bc_out_type + use FatesInterfaceMod, only : nleafage use FatesSizeAgeTypeIndicesMod, only : get_sizeage_class_index use PRTGenericMod, only : prt_global diff --git a/parameter_files/fates_params_default.cdl b/parameter_files/fates_params_default.cdl index 47c7e55900..f0917e2e84 100644 --- a/parameter_files/fates_params_default.cdl +++ b/parameter_files/fates_params_default.cdl @@ -879,10 +879,10 @@ data: fates_leaf_jmaxse = 495, 495 ; - fates_leaf_long = 0.4, 0.4, - 0.4, 0.4, - 0.4, 0.4, - 0.4, 0.4 ; + fates_leaf_long = 0.375, 0.375, + 0.375, 0.375, + 0.375, 0.375, + 0.375, 0.375 ; fates_leaf_slamax = 0.0954, 0.0954 ; @@ -896,10 +896,10 @@ data: fates_leaf_tpuse = 490, 490 ; - fates_leaf_vcmax25top = 62, 62, - 54, 54, - 46, 46, - 38, 38; + fates_leaf_vcmax25top = 50, 50, + 50, 50, + 50, 50, + 50, 50; fates_leaf_vcmaxha = 65330, 65330 ; diff --git a/parteh/PRTAllometricCarbonMod.F90 b/parteh/PRTAllometricCarbonMod.F90 index e0c81cc61d..43f250183f 100644 --- a/parteh/PRTAllometricCarbonMod.F90 +++ b/parteh/PRTAllometricCarbonMod.F90 @@ -47,7 +47,7 @@ module PRTAllometricCarbonMod use FatesConstantsMod , only : calloc_abs_error use FatesConstantsMod , only : nearzero use FatesConstantsMod , only : itrue - + use FatesConstantsMod , only : years_per_day implicit none @@ -103,8 +103,8 @@ module PRTAllometricCarbonMod integer, parameter :: icd = 1 ! Only 1 coordinate per variable - ! This is the maximum number of leaf age pools, used for allocating scratch space - integer, parameter :: max_nleafage = 4 + ! This is the maximum number of leaf age pools (used for allocating scratch space) + integer, parameter :: max_nleafage = 10 ! ------------------------------------------------------------------------------------- ! This is the core type that holds this specific @@ -328,6 +328,7 @@ subroutine DailyPRTAllometricCarbon(this) real(r8) :: struct_c_flux ! Transfer into structure at various stages [kgC] real(r8),dimension(max_nleafage) :: leaf_c0 + ! Initial value of carbon used to determine net flux real(r8) :: fnrt_c0 ! during this routine real(r8) :: sapw_c0 ! "" @@ -348,7 +349,8 @@ subroutine DailyPRTAllometricCarbon(this) real(r8) :: totalC ! total carbon allocated over alometric growth step real(r8) :: hite_out ! dummy height variable - integer :: i_var ! local index for iterating state variables + integer :: i_var ! index for iterating state variables + integer :: i_age ! index for iterating leaf ages integer :: nleafage ! number of leaf age classifications real(r8) :: leaf_age_flux ! carbon mass flux between leaf age classification pools @@ -406,8 +408,7 @@ subroutine DailyPRTAllometricCarbon(this) ! transport flux "%net_alloc" at the end. ! ----------------------------------------------------------------------------------- - nleafage = prt_global%state_descriptor(leaf_c_id)%num_pos - + nleafage = prt_global%state_descriptor(leaf_c_id)%num_pos ! Number of leaf age class leaf_c0(1:nleafage) = leaf_c(1:nleafage) ! Set initial leaf carbon fnrt_c0 = fnrt_c ! Set initial fine-root carbon @@ -427,13 +428,15 @@ subroutine DailyPRTAllometricCarbon(this) ! ----------------------------------------------------------------------------------- if(nleafage>1) then - do iage = 1,nleafage-1 - leaf_age_flux = leaf_c0(iage) * hlm_freq_day / EDPftvarcon_inst%leaf_long(ipft,iage) - leaf_c(iage) = leaf_c(iage) - leaf_age_flux - leaf_c(iage+1) = leaf_c(iage+1) + leaf_age_flux + do i_age = 1,nleafage-1 + if (EDPftvarcon_inst%leaf_long(ipft,i_age)>nearzero) then + leaf_age_flux = leaf_c0(i_age) * years_per_day / EDPftvarcon_inst%leaf_long(ipft,i_age) + leaf_c(i_age) = leaf_c(i_age) - leaf_age_flux + leaf_c(i_age+1) = leaf_c(i_age+1) + leaf_age_flux + end if end do end if - + ! ----------------------------------------------------------------------------------- ! II. Calculate target size of the biomass compartment for a given dbh. @@ -827,7 +830,7 @@ subroutine DailyPRTAllometricCarbon(this) struct_c_flux = struct_c_flux*flux_adj repro_c_flux = repro_c_flux*flux_adj - carbon_balance = carbon_balance - leaf_c_flux + carbon_balance = carbon_balance - leaf_c_flux leaf_c(iexp_leaf) = leaf_c(iexp_leaf) + leaf_c_flux carbon_balance = carbon_balance - fnrt_c_flux @@ -866,8 +869,11 @@ subroutine DailyPRTAllometricCarbon(this) ! Track the net allocations and transport from this routine - this%variables(leaf_c_id)%net_alloc(icd) = & - this%variables(leaf_c_id)%net_alloc(icd) + (leaf_c - leaf_c0) + do i_age = 1,nleafage + this%variables(leaf_c_id)%net_alloc(i_age) = & + this%variables(leaf_c_id)%net_alloc(i_age) + & + (leaf_c(i_age) - leaf_c0(i_age)) + end do this%variables(fnrt_c_id)%net_alloc(icd) = & this%variables(fnrt_c_id)%net_alloc(icd) + (fnrt_c - fnrt_c0) @@ -885,9 +891,6 @@ subroutine DailyPRTAllometricCarbon(this) this%variables(struct_c_id)%net_alloc(icd) + (struct_c - struct_c0) - - - end associate diff --git a/parteh/PRTGenericMod.F90 b/parteh/PRTGenericMod.F90 index 7f9c99dac0..e82c93240d 100644 --- a/parteh/PRTGenericMod.F90 +++ b/parteh/PRTGenericMod.F90 @@ -763,7 +763,7 @@ end subroutine CopyPRTVartypes ! ===================================================================================== - subroutine WeightedFusePRTVartypes(this,donor_prt_obj, recipient_fuse_weight, position_id) + subroutine WeightedFusePRTVartypes(this,donor_prt_obj, recipient_fuse_weight) ! This subroutine fuses two PRT objects together based on a fusion weighting ! assigned for the recipient (the class calling this) @@ -773,17 +773,18 @@ subroutine WeightedFusePRTVartypes(this,donor_prt_obj, recipient_fuse_weight, po class(prt_vartypes), intent(in), pointer :: donor_prt_obj real(r8),intent(in) :: recipient_fuse_weight ! This is the weighting ! for the recipient - integer,intent(in),optional :: position_id +! integer,intent(in),optional :: position_id ! Locals integer :: i_var ! Loop iterator over variables integer :: pos_id ! coordinate id (defaults to 1, if not position_id) - if(present(position_id)) then - pos_id = position_id - - do i_var = 1, prt_global%num_vars - + ! If no argument regarding positions is supplied, then fuse all positions + ! sequentially + do i_var = 1, prt_global%num_vars + + do pos_id = 1,prt_global%state_descriptor(i_var)%num_pos + this%variables(i_var)%val(pos_id) = recipient_fuse_weight * this%variables(i_var)%val(pos_id) + & (1.0_r8-recipient_fuse_weight) * donor_prt_obj%variables(i_var)%val(pos_id) @@ -798,41 +799,11 @@ subroutine WeightedFusePRTVartypes(this,donor_prt_obj, recipient_fuse_weight, po this%variables(i_var)%burned(pos_id) = recipient_fuse_weight * this%variables(i_var)%burned(pos_id) + & (1.0_r8-recipient_fuse_weight) * donor_prt_obj%variables(i_var)%burned(pos_id) - - end do - else - - ! If no argument regarding positions is supplied, then fuse all positions - ! sequentially - do i_var = 1, prt_global%num_vars - - do pos_id = 1,prt_global%state_descriptor(leaf_c_id)%num_pos - - this%variables(i_var)%val(pos_id) = recipient_fuse_weight * this%variables(i_var)%val(pos_id) + & - (1.0_r8-recipient_fuse_weight) * donor_prt_obj%variables(i_var)%val(pos_id) - - this%variables(i_var)%val0(pos_id) = recipient_fuse_weight * this%variables(i_var)%val0(pos_id) + & - (1.0_r8-recipient_fuse_weight) * donor_prt_obj%variables(i_var)%val0(pos_id) - - this%variables(i_var)%net_alloc(pos_id) = recipient_fuse_weight * this%variables(i_var)%net_alloc(pos_id) + & - (1.0_r8-recipient_fuse_weight) * donor_prt_obj%variables(i_var)%net_alloc(pos_id) - - this%variables(i_var)%turnover(pos_id) = recipient_fuse_weight * this%variables(i_var)%turnover(pos_id) + & - (1.0_r8-recipient_fuse_weight) * donor_prt_obj%variables(i_var)%turnover(pos_id) - - this%variables(i_var)%burned(pos_id) = recipient_fuse_weight * this%variables(i_var)%burned(pos_id) + & - (1.0_r8-recipient_fuse_weight) * donor_prt_obj%variables(i_var)%burned(pos_id) - - end do end do - - end if - - - end do + this%ode_opt_step = recipient_fuse_weight * this%ode_opt_step + & (1.0_r8-recipient_fuse_weight) * donor_prt_obj%ode_opt_step diff --git a/parteh/PRTLossFluxesMod.F90 b/parteh/PRTLossFluxesMod.F90 index f7d0c9b3c0..1720e5a1a6 100644 --- a/parteh/PRTLossFluxesMod.F90 +++ b/parteh/PRTLossFluxesMod.F90 @@ -24,6 +24,7 @@ module PRTLossFluxesMod use FatesConstantsMod, only : i4 => fates_int use FatesConstantsMod, only : nearzero use FatesConstantsMod, only : calloc_abs_error + use FatesConstantsMod, only : itrue use FatesGlobals , only : endrun => fates_endrun use FatesGlobals , only : fates_log use shr_log_mod , only : errMsg => shr_log_errMsg @@ -524,7 +525,7 @@ end subroutine DeciduousTurnoverSimpleRetranslocation ! ==================================================================================== - subroutine PRTMaintTurnover(prt,ipft,senescent_frac) + subroutine PRTMaintTurnover(prt,ipft) ! --------------------------------------------------------------------------------- ! Generic subroutine (wrapper) calling specialized routines handling @@ -534,7 +535,7 @@ subroutine PRTMaintTurnover(prt,ipft,senescent_frac) integer,intent(in) :: ipft if ( int(EDPftvarcon_inst%turnover_retrans_mode(ipft)) == 1 ) then - call MaintTurnoverSimpleRetranslocation(prt,ipft,senescent_frac) + call MaintTurnoverSimpleRetranslocation(prt,ipft) else write(fates_log(),*) 'A maintenance/retranslocation mode was specified' write(fates_log(),*) 'that is unknown.' @@ -624,7 +625,7 @@ subroutine MaintTurnoverSimpleRetranslocation(prt,ipft) ! timescale for the senescent pool. aclass_sen_id = size(EDPftvarcon_inst%leaf_long(ipft,:)) - ! Only EVERGREENS HAVE MAINTENANCE LEAF TURNOVER + ! Only EVERGREENS HAVE MAINTENANCE LEAF TURNOVER ? ! ------------------------------------------------------------------------------------- if ( (EDPftvarcon_inst%leaf_long(ipft,aclass_sen_id) > nearzero ) .and. & (EDPftvarcon_inst%evergreen(ipft) == itrue) ) then @@ -676,9 +677,14 @@ subroutine MaintTurnoverSimpleRetranslocation(prt,ipft) ! Hypotheses 1 & 2 assume that the leaf pools are statified by age ! We only generate turnover from the last (senescing) position - if( (prt_global%hyp_id .le. 2) .and. & - (organ_id .eq. leaf_organ) ) then - ipos_1 = prt_global%state_descriptor(i_var)%num_pos + if((organ_id .eq. leaf_organ)) then + if (prt_global%hyp_id .le. 2) then + ipos_1 = prt_global%state_descriptor(i_var)%num_pos + else + write(fates_log(),*) 'Unhandled Leaf maintenance turnover condition' + write(fates_log(),*) 'for PARTEH hypothesis id: ',prt_global%hyp_id + call endrun(msg=errMsg(__FILE__, __LINE__)) + end if else ipos_1 = 1 end if From c1ef6794b4e3da7a11c8eec170fec1b18c76077c Mon Sep 17 00:00:00 2001 From: chonggang xu Date: Tue, 18 Dec 2018 17:53:30 -0800 Subject: [PATCH 11/34] Updated leaf_longevity in trimming --- biogeochem/EDPhysiologyMod.F90 | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index a1438c406b..cff0dee414 100644 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -328,9 +328,12 @@ subroutine trim_canopy( currentSite ) currentCohort%leaf_cost = currentCohort%leaf_cost * & (EDPftvarcon_inst%grperc(ipft) + 1._r8) else !evergreen costs + ! Leaf cost at leaf level z accounting for sla profile currentCohort%leaf_cost = 1.0_r8/(sla_levleaf* & - EDPftvarcon_inst%leaf_long(ipft,0)*1000.0_r8) !convert from sla in m2g-1 to m2kg-1 + sum(EDPftvarcon_inst%leaf_long(ipft,:))*1000.0_r8) !convert from sla in m2g-1 to m2kg-1 + + if ( int(EDPftvarcon_inst%allom_fmode(ipft)) .eq. 1 ) then ! if using trimmed leaf for fine root biomass allometry, add the cost of the root increment ! to the leaf increment; otherwise do not. From 71fefa126ee7cd984a1e3dc6bf31c0de05e35535 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Wed, 19 Dec 2018 12:04:24 -0800 Subject: [PATCH 12/34] Cleaning up the leaf-age class implementation --- biogeochem/EDCohortDynamicsMod.F90 | 32 ++++++------------------------ biogeochem/EDPatchDynamicsMod.F90 | 1 - biogeochem/EDPhysiologyMod.F90 | 2 +- main/EDMainMod.F90 | 6 ++++-- main/EDPftvarcon.F90 | 4 ---- main/FatesRestartInterfaceMod.F90 | 16 ++++++++++----- 6 files changed, 22 insertions(+), 39 deletions(-) diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index d1427388f6..6d57594fca 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -15,6 +15,7 @@ module EDCohortDynamicsMod use FatesInterfaceMod , only : hlm_days_per_year use FatesInterfaceMod , only : nleafage use EDPftvarcon , only : EDPftvarcon_inst + use FatesParameterDerivedMod, only : param_derived use EDTypesMod , only : ed_site_type, ed_patch_type, ed_cohort_type use EDTypesMod , only : nclmax use EDTypesMod , only : ncwd @@ -26,6 +27,7 @@ module EDCohortDynamicsMod use EDTypesMod , only : equal_leaf_aclass use EDTypesMod , only : first_leaf_aclass use EDTypesMod , only : nan_leaf_aclass + use EDTypesMod , only : max_nleafage use FatesInterfaceMod , only : hlm_use_planthydro use FatesInterfaceMod , only : hlm_parteh_mode use FatesPlantHydraulicsMod, only : FuseCohortHydraulics @@ -219,7 +221,7 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, hite, dbh, bleaf, bfine ! This sets things like vcmax25top, that depend on the ! leaf age fractions - call UpdateCohortBioPhysRates(currentCohort) + call UpdateCohortBioPhysRates(new_cohort) call sizetype_class_index(new_cohort%dbh,new_cohort%pft, & new_cohort%size_class,new_cohort%size_by_pft_class) @@ -242,11 +244,11 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, hite, dbh, bleaf, bfine new_cohort%treelai = tree_lai(bleaf, new_cohort%pft, new_cohort%c_area, & new_cohort%n, new_cohort%canopy_layer, & - patchptr%canopy_layer_tlai ) + patchptr%canopy_layer_tlai,new_cohort%vcmax25top ) new_cohort%treesai = tree_sai(new_cohort%pft, new_cohort%dbh, new_cohort%canopy_trim, & new_cohort%c_area, new_cohort%n, new_cohort%canopy_layer, & - patchptr%canopy_layer_tlai, new_cohort%treelai ) + patchptr%canopy_layer_tlai, new_cohort%treelai,new_cohort%vcmax25top ) new_cohort%lai = new_cohort%treelai * new_cohort%c_area/patchptr%area @@ -896,30 +898,8 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) ! Leaf biophysical rates (use leaf mass weighting) ! ----------------------------------------------------------------- - leaf_c_next = nextc%prt%GetState(leaf_organ, all_carbon_elements)*nextc%n - leaf_c_curr = currentCohort%prt%GetState(leaf_organ, all_carbon_elements)*currentCohort%n + call UpdateCohortBioPhysRates(currentCohort) - if( (leaf_c_next + leaf_c_curr) > nearzero ) then - - currentCohort%vcmax25top = (leaf_c_curr*currentCohort%vcmax25top + & - leaf_c_next*nextc%vcmax25top)/(leaf_c_next+leaf_c_curr) - currentCohort%jmax25top = (leaf_c_curr*currentCohort%jmax25top + & - leaf_c_next*nextc%jmax25top)/(leaf_c_next+leaf_c_curr) - currentCohort%tpu25top = (leaf_c_curr*currentCohort%tpu25top + & - leaf_c_next*nextc%tpu25top)/(leaf_c_next+leaf_c_curr) - currentCohort%kp25top = (leaf_c_curr*currentCohort%kp25top + & - leaf_c_next*nextc%kp25top)/(leaf_c_next+leaf_c_curr) - - else - - currentCohort%vcmax25top = 0._r8 - currentCohort%jmax25top = 0._r8 - currentCohort%tpu25top = 0._r8 - currentCohort%kp25top = 0._r8 - - end if - - ! keep track of the size class bins so that we can monitor growth fluxes ! compare the values. if they are the same, then nothing needs to be done. if not, track the diagnostic flux if (currentCohort%size_class_lasttimestep .ne. nextc%size_class_lasttimestep ) then diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 458a949daf..5f140a63fe 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -22,7 +22,6 @@ module EDPatchDynamicsMod use FatesInterfaceMod , only : bc_in_type use FatesInterfaceMod , only : hlm_days_per_year use FatesInterfaceMod , only : numpft - use FatesInterfaceMod , only : nleafage use FatesGlobals , only : endrun => fates_endrun use FatesConstantsMod , only : r8 => fates_r8 use FatesConstantsMod , only : itrue diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index cff0dee414..9085777c16 100644 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -299,7 +299,7 @@ subroutine trim_canopy( currentSite ) ! Calculate sla_levleaf following the sla profile with overlying leaf area ! Scale for leaf nitrogen profile - kn = decay_coeff_kn(ipft) + kn = decay_coeff_kn(ipft,currentCohort%vcmax25top) ! Nscaler value at leaf level z nscaler_levleaf = exp(-kn * cumulative_lai) ! Sla value at leaf level z after nitrogen profile scaling (m2/gC) diff --git a/main/EDMainMod.F90 b/main/EDMainMod.F90 index 353ff0f339..e8ae1005ca 100644 --- a/main/EDMainMod.F90 +++ b/main/EDMainMod.F90 @@ -347,8 +347,10 @@ subroutine ed_integrate_state_variables(currentSite, bc_in ) call currentCohort%prt%DailyPRT() call currentCohort%prt%CheckMassConservation(ft,5) - ! Update the leaf biophysical rates based on leaf ages - ! leaf ages should be fully updated for each cohort at this point + ! Update the leaf biophysical rates based on proportion of leaf + ! mass in the different leaf age classes. Following growth + ! and turnover, these proportions won't change again. This + ! routine is also called following fusion call UpdateCohortBioPhysRates(currentCohort) diff --git a/main/EDPftvarcon.F90 b/main/EDPftvarcon.F90 index 1825122ad3..1c3b806d57 100644 --- a/main/EDPftvarcon.F90 +++ b/main/EDPftvarcon.F90 @@ -862,10 +862,6 @@ subroutine Receive_PFT(this, fates_params) call fates_params%RetreiveParameterAllocate(name=name, & data=this%slatop) - name = 'fates_leaf_long' - call fates_params%RetreiveParameterAllocate(name=name, & - data=this%leaf_long) - name = 'fates_roota_par' call fates_params%RetreiveParameterAllocate(name=name, & data=this%roota_par) diff --git a/main/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index 30e0d76d73..6a603baeb8 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -16,7 +16,7 @@ module FatesRestartInterfaceMod use FatesInterfaceMod, only : bc_out_type use FatesInterfaceMod, only : nleafage use FatesSizeAgeTypeIndicesMod, only : get_sizeage_class_index - + use EDCohortDynamicsMod, only : UpdateCohortBioPhysRates use PRTGenericMod, only : prt_global @@ -1799,7 +1799,17 @@ subroutine get_restart_vectors(this, nc, nsites, sites) this%rvars(ir_prt_var)%r81d(io_idx_co) end do end do + + + call UpdateCohortBioPhysRates(ccohort) + + !ccohort%vcmax25top + !ccohort%jmax25top + !ccohort%tpu25top + !ccohort%kp25top + + ccohort%canopy_layer = rio_canopy_layer_co(io_idx_co) ccohort%canopy_layer_yesterday = rio_canopy_layer_yesterday_co(io_idx_co) ccohort%canopy_trim = rio_canopy_trim_co(io_idx_co) @@ -1821,10 +1831,6 @@ subroutine get_restart_vectors(this, nc, nsites, sites) ccohort%fmort = rio_fmort_co(io_idx_co) ccohort%frmort = rio_frmort_co(io_idx_co) - ! THE ACTUAL LEAF AGES MUST BE READ IN AFTER NEW HYDRO PR - ! (THAT PR HAS THE VECTOR READ-IN MACHINERY) - ccohort%frac_leaf_aclass(1:nleafage) = 1._r8 / real(nleafage,r8) - !Logging ccohort%lmort_direct = rio_lmort_direct_co(io_idx_co) ccohort%lmort_collateral = rio_lmort_collateral_co(io_idx_co) From 8202aa755d5aacc7aae5ef5c325691c808b400a1 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Wed, 19 Dec 2018 14:30:23 -0800 Subject: [PATCH 13/34] Refactors of site level drought and cold status --- biogeochem/EDPhysiologyMod.F90 | 127 +++++++++++++++++------------- main/EDInitMod.F90 | 57 ++++++-------- main/EDTypesMod.F90 | 5 +- main/FatesInventoryInitMod.F90 | 34 ++++---- main/FatesRestartInterfaceMod.F90 | 55 ++++++++++--- 5 files changed, 157 insertions(+), 121 deletions(-) diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index 9085777c16..c7a32e8355 100644 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -496,11 +496,11 @@ subroutine phenology( currentSite, bc_in ) !2) The leaves should not be on already !3) There should have been at least on chilling day in the counting period. if (currentSite%ED_GDD_site > gdd_threshold)then - if (currentSite%status == leaves_off) then + if (currentSite%is_cold) then if (currentSite%ncd >= 1) then - currentSite%status = leaves_on ! alter status of site to 'leaves on' - ! NOTE(bja, 2015-01) should leafondate = model_day to be consistent with leaf off? - currentSite%leafondate = t !record leaf on date + ! NOTE(bja, 2015-01) should leafondate = model_day to be consistent with leaf off? + currentSite%is_cold = .false. + currentSite%leafondate = t ! record leaf on date if ( debug ) write(fates_log(),*) 'leaves on' endif !ncd endif !status @@ -514,13 +514,14 @@ subroutine phenology( currentSite, bc_in ) !1) have exceeded the number of cold days threshold !2) have exceeded the minimum leafon time. !3) The leaves should not be off already - !4) The day of the year should be larger than the counting period. (not sure if we need this/if it will break the restarting) + !4) The day of the year should be larger than the counting period. + ! (not sure if we need this/if it will break the restarting) if (ncolddays > ED_val_phen_ncolddayslim)then if (timesinceleafon > ED_val_phen_mindayson)then - if (currentSite%status == leaves_on)then - currentSite%status = leaves_off !alter status of site to 'leaves off' - currentSite%leafoffdate = hlm_model_day !record leaf off date + if (.not.currentSite%is_cold)then + currentSite%is_cold = .true. ! now cold enough for leaf-off + currentSite%leafoffdate = hlm_model_day ! record leaf off date if ( debug ) write(fates_log(),*) 'leaves off' endif endif @@ -528,8 +529,8 @@ subroutine phenology( currentSite, bc_in ) !LEAF OFF: COLD LIFESPAN THRESHOLD if(timesinceleafoff > 400)then !remove leaves after a whole year when there is no 'off' period. - if(currentSite%status == leaves_on)then - currentSite%status = leaves_off !alter status of site to 'leaves off' + if(.not.currentSite%is_cold)then + currentSite%is_cold = .true. ! cold enough for leaf-off currentSite%leafoffdate = hlm_model_day !record leaf off date if ( debug ) write(fates_log(),*) 'leaves off' endif @@ -537,7 +538,7 @@ subroutine phenology( currentSite, bc_in ) !-----------------Drought Phenology--------------------! ! Principles of drought-deciduos phenology model... - ! The 'dstatus' flag is 2 when leaves are on, and 1 when leaves area off. + ! The 'is_drought' flag is false when leaves are on, and true when leaves area off. ! The following sets those site-level flags, which are acted on in phenology_deciduos. ! A* The leaves live for either the length of time the soil moisture is over the threshold ! or the lifetime of the leaves, whichever is shorter. @@ -569,10 +570,14 @@ subroutine phenology( currentSite, bc_in ) enddo currentSite%water_memory(1) = bc_in%h2o_liqvol_sl(1) !waterstate_inst%h2osoi_vol_col(coli,1) - !In drought phenology, we often need to force the leaves to stay on or off as moisture fluctuates... + ! In drought phenology, we often need to force the leaves to stay on or off as moisture fluctuates... + ! Here we incremend how long the leaves have been off; + ! We set the default assumption that no time has elapsed, but if drought + ! status is true, then we update the time + ! If the leaves are off. How long have they been off? + ! leaves have come on, but last year, so at a later date than now. timesincedleafoff = 0 - if (currentSite%dstatus == leaves_off)then !the leaves are off. How long have they been off? - !leaves have come on, but last year, so at a later date than now. + if ( currentSite%is_drought )then if (currentSite%dleafoffdate > 0.and.currentSite%dleafoffdate > t)then timesincedleafoff = t + (360 - currentSite%dleafoffdate) else @@ -582,7 +587,7 @@ subroutine phenology( currentSite, bc_in ) timesincedleafon = 0 !the leaves are on. How long have they been on? - if (currentSite%dstatus == leaves_on)then + if ( .not.currentSite%is_drought )then !leaves have come on, but last year, so at a later date than now. if (currentSite%dleafondate > 0.and.currentSite%dleafondate > t)then timesincedleafon = t + (360 - currentSite%dleafondate) @@ -595,40 +600,40 @@ subroutine phenology( currentSite, bc_in ) !Here, we used a window of oppurtunity to determine if we are close to the time when then leaves came on last year if ((t >= currentSite%dleafondate - 30.and.t <= currentSite%dleafondate + 30).or.(t > 360 - 15.and. & currentSite%dleafondate < 15))then ! are we in the window? - ! TODO: CHANGE THIS MATH, MOVE THE DENOMENATOR OUTSIDE OF THE SUM (rgk 01-2017) - if (sum(currentSite%water_memory(1:numWaterMem)/real(numWaterMem,r8)) & - >= ED_val_phen_drought_threshold.and.currentSite%dstatus == leaves_off.and.t >= 10)then + + if ( sum(currentSite%water_memory(1:numWaterMem))/real(numWaterMem,r8) & + >= ED_val_phen_drought_threshold .and. & + currentSite%is_drought .and. & + (t >= 10) ) then ! leave some minimum time between leaf off and leaf on to prevent 'flickering'. if (timesincedleafoff > ED_val_phen_doff_time)then - currentSite%dstatus = leaves_on !alter status of site to 'leaves on' - currentSite%dleafondate = t !record leaf on date + currentSite%is_drought = .false. ! end the drought + currentSite%dleafondate = t !record leaf on date endif endif endif - !we still haven't done budburst by end of window - if (t == currentSite%dleafondate+30.and.currentSite%dstatus == leaves_off)then - currentSite%dstatus = leaves_on ! force budburst! - currentSite%dleafondate = t ! record leaf on date + ! we still haven't done budburst by end of window + if (t == currentSite%dleafondate+30 .and. currentSite%is_drought)then + currentSite%is_drought = .false. ! force budburst! + currentSite%dleafondate = t ! record leaf on date endif !LEAF OFF: DROUGHT DECIDUOUS LIFESPAN - if the leaf gets to the end of its useful life. A*, E* - if (currentSite%dstatus == leaves_on.and.t >= 10)then !D* + if ( .not.currentSite%is_drought .and. (t >= 10) ) then !D* !Are the leaves at the end of their lives? - !FIX(RF,0401014)- this is hardwiring.... - !FIX(RGK:changed from hard-coded pft 7 leaf lifespan to labeled constant (1 year) if ( timesincedleafon > canopy_leaf_lifespan )then - currentSite%dstatus = leaves_off !alter status of site to 'leaves off' - currentSite%dleafoffdate = t !record leaf on date + currentSite%is_drought = .true. !alter status of site to 'leaves off' + currentSite%dleafoffdate = t !record leaf off date endif endif !LEAF OFF: DROUGHT DECIDUOUS DRYNESS - if the soil gets too dry, and the leaves have already been on a while... - if (currentSite%dstatus == leaves_on.and.t >= 10)then !D* + if ( .not.currentSite%is_drought .and. (t >= 10) ) then !D* if (sum(currentSite%water_memory(1:10)/10._r8) <= ED_val_phen_drought_threshold)then if (timesincedleafon > 100)then !B* Have the leaves been on for some reasonable length of time? To prevent flickering. - currentSite%dstatus = leaves_off !alter status of site to 'leaves on' - currentSite%dleafoffdate = t !record leaf on date + currentSite%is_drought = .true. !alter status of site to 'leaves on' + currentSite%dleafoffdate = t !record leaf on date endif endif endif @@ -675,9 +680,12 @@ subroutine phenology_leafonoff(currentSite) store_c = currentCohort%prt%GetState(store_organ, all_carbon_elements) leaf_c = currentCohort%prt%GetState(leaf_organ, all_carbon_elements) - !COLD LEAF ON + ! COLD LEAF ON + ! The site level flags signify that it is no-longer too cold + ! for leaves. Time to signal flushing + if (EDPftvarcon_inst%season_decid(ipft) == 1)then - if (currentSite%status == leaves_on)then !we have just moved to leaves being on . + if ( .not.currentSite%is_cold )then !we have just moved to leaves being on . if (currentCohort%status_coh == leaves_off)then !Are the leaves currently off? currentCohort%status_coh = leaves_on ! Leaves are on, so change status to ! stop flow of carbon out of bstore. @@ -698,8 +706,11 @@ subroutine phenology_leafonoff(currentSite) endif !pft phenology endif ! growing season - !COLD LEAF OFF - if (currentSite%status == leaves_off)then !past leaf drop day? Leaves still on tree? + ! COLD LEAF OFF + ! The site level flag signifies that it is now too cold for + ! deciduous leaves. Time to drop if they have not already + + if (currentSite%is_cold) then if (currentCohort%status_coh == leaves_on)then ! leaves have not dropped ! This sets the cohort to the "leaves off" flag @@ -722,10 +733,13 @@ subroutine phenology_leafonoff(currentSite) endif !currentSite status endif !season_decid - !DROUGHT LEAF ON + ! DROUGHT LEAF ON + ! Site level flag indicates it is no longer in drought condition + ! deciduous plants can flush + if (EDPftvarcon_inst%stress_decid(ipft) == 1)then - if ( currentSite%dstatus == leaves_on )then + if ( .not.currentSite%is_drought ) then ! we have just moved to leaves being on . if (currentCohort%status_coh == leaves_off)then @@ -752,8 +766,11 @@ subroutine phenology_leafonoff(currentSite) endif !currentCohort status again? endif !currentSite status - !DROUGHT LEAF OFF - if (currentSite%dstatus == leaves_off)then + ! DROUGHT LEAF OFF + ! Site level flag indicates a drought condition is in effect + ! deciduous plants should drop leaves if have not already + + if ( currentSite%is_drought ) then if (currentCohort%status_coh == leaves_on)then ! leaves have not dropped ! This sets the cohort to the "leaves off" flag @@ -939,16 +956,18 @@ subroutine seed_germination( currentSite, currentPatch ) currentPatch%seed_germination(p) = min(currentSite%seed_bank(p) * & EDPftvarcon_inst%germination_timescale(p),max_germination) !set the germination only under the growing season...c.xu - if (EDPftvarcon_inst%season_decid(p) == 1.and.currentSite%status == leaves_off)then + if (EDPftvarcon_inst%season_decid(p) == itrue .and. currentSite%is_cold)then currentPatch%seed_germination(p) = 0.0_r8 endif - if (EDPftvarcon_inst%stress_decid(p) == 1.and.currentSite%dstatus == leaves_off)then + if (EDPftvarcon_inst%stress_decid(p) == itrue .and. currentSite%is_drought)then currentPatch%seed_germination(p) = 0.0_r8 endif enddo end subroutine seed_germination + ! ===================================================================================== + subroutine recruitment( currentSite, currentPatch, bc_in ) ! ! !DESCRIPTION: @@ -996,25 +1015,27 @@ subroutine recruitment( currentSite, currentPatch, bc_in ) call bdead_allom(b_agw,b_bgw,b_sapwood,ft,b_dead) call bstore_allom(temp_cohort%dbh,ft,temp_cohort%canopy_trim,b_store) + ! Default assumption is that leaves are on + cohortstatus = leaves_on temp_cohort%laimemory = 0.0_r8 - if (EDPftvarcon_inst%season_decid(temp_cohort%pft) == 1.and.currentSite%status == leaves_off)then + + ! But if the plant is seasonally (cold) deciduous, and the site status is flagged + ! as "cold", then set the cohort's status to leaves_off, and remember the leaf biomass + if (EDPftvarcon_inst%season_decid(ft) == itrue .and. currentSite%is_cold)then temp_cohort%laimemory = b_leaf b_leaf = 0.0_r8 + cohortstatus = leaves_off endif - if (EDPftvarcon_inst%stress_decid(temp_cohort%pft) == 1.and.currentSite%dstatus == leaves_off)then + + ! Or.. if the plant is drought deciduous, and the site status is flagged as + ! "in a drought", then likewise, set the cohort's status to leaves_off, and remember leaf + ! biomass + if (EDPftvarcon_inst%stress_decid(ft) == itrue .and. currentSite%is_drought )then temp_cohort%laimemory = b_leaf b_leaf = 0.0_r8 + cohortstatus = leaves_off endif - cohortstatus = currentSite%status - if (EDPftvarcon_inst%stress_decid(ft) == 1)then !drought decidous, override status. - cohortstatus = currentSite%dstatus - endif - - if (EDPftvarcon_inst%evergreen(ft) == 1) then - temp_cohort%laimemory = 0._r8 - cohortstatus = leaves_on - endif 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 & diff --git a/main/EDInitMod.F90 b/main/EDInitMod.F90 index 516a1ac00f..73e2470249 100644 --- a/main/EDInitMod.F90 +++ b/main/EDInitMod.F90 @@ -22,6 +22,8 @@ module EDInitMod use EDTypesMod , only : init_spread_near_bare_ground use EDTypesMod , only : init_spread_inventory use EDTypesMod , only : first_leaf_aclass + use EDTypesMod , only : leaves_on + use EDTypesMod , only : leaves_off use FatesInterfaceMod , only : bc_in_type use FatesInterfaceMod , only : hlm_use_planthydro use FatesInterfaceMod , only : hlm_use_inventory_init @@ -101,8 +103,8 @@ subroutine zero_site( site_in ) site_in%total_burn_flux_to_atm = 0._r8 ! PHENOLOGY - site_in%status = 0 ! are leaves in this pixel on or off? - site_in%dstatus = 0 + site_in%is_cold = .false. ! Is cold deciduous leaf-off triggered? + site_in%is_drought = .false. ! Is drought deciduous leaf-off triggered? site_in%ED_GDD_site = nan ! growing degree days site_in%ncd = nan ! no chilling days site_in%last_n_days(:) = 999 ! record of last 10 days temperature for senescence model. @@ -170,10 +172,10 @@ subroutine set_site_properties( nsites, sites) integer :: s real(r8) :: leafon real(r8) :: leafoff - real(r8) :: stat + logical :: stat real(r8) :: NCD real(r8) :: GDD - real(r8) :: dstat + logical :: dstat real(r8) :: acc_NI real(r8) :: watermem integer :: dleafoff @@ -186,9 +188,9 @@ subroutine set_site_properties( nsites, sites) GDD = 30.0_r8 leafon = 100.0_r8 leafoff = 300.0_r8 - stat = 2 + stat = .false. acc_NI = 0.0_r8 - dstat = 2 + dstat = .false. dleafoff = 300 dleafon = 100 watermem = 0.5_r8 @@ -199,9 +201,9 @@ subroutine set_site_properties( nsites, sites) GDD = 0.0_r8 leafon = 0.0_r8 leafoff = 0.0_r8 - stat = 1 + stat = .false. acc_NI = 0.0_r8 - dstat = 2 + dstat = .false. dleafoff = 300 dleafon = 100 watermem = 0.5_r8 @@ -220,9 +222,8 @@ subroutine set_site_properties( nsites, sites) sites(s)%water_memory(1:numWaterMem) = watermem end if - sites(s)%status = stat - !start off with leaves off to initialise - sites(s)%dstatus= dstat + sites(s)%is_cold = stat + sites(s)%is_drought = dstat sites(s)%acc_NI = acc_NI sites(s)%frac_burnt = 0.0_r8 @@ -420,29 +421,19 @@ subroutine init_cohorts( site_in, patch_in, bc_in) call bstore_allom(temp_cohort%dbh, pft, temp_cohort%canopy_trim, b_store) - - if( EDPftvarcon_inst%evergreen(pft) == 1) then - temp_cohort%laimemory = 0._r8 - cstatus = 2 - endif - - if( EDPftvarcon_inst%season_decid(pft) == 1 ) then !for dorment places - if(site_in%status == 2)then - temp_cohort%laimemory = 0.0_r8 - else - temp_cohort%laimemory = b_leaf - endif - ! reduce biomass according to size of store, this will be recovered when elaves com on. - cstatus = site_in%status + temp_cohort%laimemory = 0._r8 + cstatus = leaves_on + + if( EDPftvarcon_inst%season_decid(pft) == itrue .and. site_in%is_cold ) then + temp_cohort%laimemory = b_leaf + b_leaf = 0._r8 + cstatus = leaves_off endif - - if ( EDPftvarcon_inst%stress_decid(pft) == 1 ) then - if(site_in%dstatus == 2)then - temp_cohort%laimemory = 0.0_r8 - else - temp_cohort%laimemory = b_leaf - endif - cstatus = site_in%dstatus + + if ( EDPftvarcon_inst%stress_decid(pft) == itrue .and. site_in%is_drought ) then + temp_cohort%laimemory = b_leaf + b_leaf = 0._r8 + cstatus = leaves_off endif if ( debug ) write(fates_log(),*) 'EDInitMod.F90 call create_cohort ' diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index a75b86fa45..7efc3f5cd9 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -75,7 +75,6 @@ module EDTypesMod ! its leaves and should not be trying to allocate ! towards any growth. - ! Switches that turn on/off ED dynamics process (names are self explanatory) ! IMPORTANT NOTE!!! THESE SWITCHES ARE EXPERIMENTAL. ! THEY SHOULD CORRECTLY TURN OFF OR ON THE PROCESS, BUT.. THERE ARE VARIOUS @@ -607,8 +606,8 @@ module EDTypesMod ! PHENOLOGY real(r8) :: ED_GDD_site ! ED Phenology growing degree days. - integer :: status ! are leaves in this pixel on or off for cold decid - integer :: dstatus ! are leaves in this pixel on or off for drought decid + logical :: is_cold ! is this site/column in a cold-status where its cohorts drop leaves? + logical :: is_drought ! is this site/column in a drought-status where its cohorts drop leaves? real(r8) :: ncd ! no chilling days:- real(r8) :: last_n_days(senes) ! record of last 10 days temperature for senescence model. deg C integer :: leafondate ! doy of leaf on:- diff --git a/main/FatesInventoryInitMod.F90 b/main/FatesInventoryInitMod.F90 index 657d00beab..338ac8618d 100644 --- a/main/FatesInventoryInitMod.F90 +++ b/main/FatesInventoryInitMod.F90 @@ -25,6 +25,7 @@ module FatesInventoryInitMod ! FATES GLOBALS use FatesConstantsMod, only : r8 => fates_r8 use FatesConstantsMod, only : pi_const + use FatesConstantsMod, only : itrue use FatesGlobals , only : endrun => fates_endrun use FatesGlobals , only : fates_log use FatesInterfaceMod, only : bc_in_type @@ -34,6 +35,8 @@ module FatesInventoryInitMod use EDTypesMod , only : ed_cohort_type use EDTypesMod , only : area use EDTypesMod , only : equal_leaf_aclass + use EDTypesMod , only : leaves_on + use EDTypesMod , only : leaves_off use EDPftvarcon , only : EDPftvarcon_inst @@ -909,28 +912,19 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & call bstore_allom(temp_cohort%dbh, c_pft, temp_cohort%canopy_trim, b_store) - if( EDPftvarcon_inst%evergreen(c_pft) == 1) then - temp_cohort%laimemory = 0._r8 - cstatus = 2 - endif - - if( EDPftvarcon_inst%season_decid(c_pft) == 1 ) then !for dorment places - if(csite%status == 2)then - temp_cohort%laimemory = 0.0_r8 - else - temp_cohort%laimemory = b_leaf - endif - ! reduce biomass according to size of store, this will be recovered when elaves com on. - cstatus = csite%status + temp_cohort%laimemory = 0._r8 + cstatus = leaves_on + + if( EDPftvarcon_inst%season_decid(c_pft) == itrue .and. csite%is_cold ) then + temp_cohort%laimemory = b_leaf + b_leaf = 0._r8 + cstatus = leaves_off endif - if ( EDPftvarcon_inst%stress_decid(c_pft) == 1 ) then - if(csite%dstatus == 2)then - temp_cohort%laimemory = 0.0_r8 - else - temp_cohort%laimemory = b_leaf - endif - cstatus = csite%dstatus + if ( EDPftvarcon_inst%stress_decid(c_pft) == itrue .and. csite%is_drought ) then + temp_cohort%laimemory = b_leaf + b_leaf = 0._r8 + cstatus = leaves_off endif ! Since spread is a canopy level calculation, we need to provide an initial guess here. diff --git a/main/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index 6a603baeb8..933d35d497 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -505,12 +505,12 @@ subroutine define_restart_vars(this, initialize_variables) flushval = flushzero, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_oldstock_si ) - call this%set_restart_var(vname='fates_cold_dec_status', vtype=site_r8, & - long_name='status flag for cold deciduous plants', units='unitless', flushval = flushzero, & + call this%set_restart_var(vname='fates_cold_dec_status', vtype=site_int, & + long_name='status flag for cold deciduous plants', units='unitless', flushval = flushinvalid, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_cd_status_si ) - call this%set_restart_var(vname='fates_drought_dec_status', vtype=site_r8, & - long_name='status flag for drought deciduous plants', units='unitless', flushval = flushzero, & + call this%set_restart_var(vname='fates_drought_dec_status', vtype=site_int, & + long_name='status flag for drought deciduous plants', units='unitless', flushval = flushinvalid, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_dd_status_si ) call this%set_restart_var(vname='fates_chilling_days', vtype=site_r8, & @@ -1086,8 +1086,8 @@ subroutine set_restart_vectors(this,nc,nsites,sites) associate( rio_npatch_si => this%rvars(ir_npatch_si)%int1d, & rio_old_stock_si => this%rvars(ir_oldstock_si)%r81d, & - rio_cd_status_si => this%rvars(ir_cd_status_si)%r81d, & - rio_dd_status_si => this%rvars(ir_dd_status_si)%r81d, & + rio_cd_status_si => this%rvars(ir_cd_status_si)%int1d, & + rio_dd_status_si => this%rvars(ir_dd_status_si)%int1d, & rio_nchill_days_si => this%rvars(ir_nchill_days_si)%r81d, & rio_leafondate_si => this%rvars(ir_leafondate_si)%r81d, & rio_leafoffdate_si => this%rvars(ir_leafoffdate_si)%r81d, & @@ -1358,8 +1358,18 @@ subroutine set_restart_vectors(this,nc,nsites,sites) enddo ! cpatch do while rio_old_stock_si(io_idx_si) = sites(s)%old_stock - rio_cd_status_si(io_idx_si) = sites(s)%status - rio_dd_status_si(io_idx_si) = sites(s)%dstatus + + if(sites(s)%is_cold) then + rio_cd_status_si(io_idx_si) = itrue + else + rio_cd_status_si(io_idx_si) = ifalse + end if + if(sites(s)%is_drought) then + rio_dd_status_si(io_idx_si) = itrue + else + rio_dd_status_si(io_idx_si) = ifalse + end if + rio_nchill_days_si(io_idx_si) = sites(s)%ncd rio_leafondate_si(io_idx_si) = sites(s)%leafondate rio_leafoffdate_si(io_idx_si) = sites(s)%leafoffdate @@ -1663,8 +1673,8 @@ subroutine get_restart_vectors(this, nc, nsites, sites) associate( rio_npatch_si => this%rvars(ir_npatch_si)%int1d, & rio_old_stock_si => this%rvars(ir_oldstock_si)%r81d, & - rio_cd_status_si => this%rvars(ir_cd_status_si)%r81d, & - rio_dd_status_si => this%rvars(ir_dd_status_si)%r81d, & + rio_cd_status_si => this%rvars(ir_cd_status_si)%int1d, & + rio_dd_status_si => this%rvars(ir_dd_status_si)%int1d, & rio_nchill_days_si => this%rvars(ir_nchill_days_si)%r81d, & rio_leafondate_si => this%rvars(ir_leafondate_si)%r81d, & rio_leafoffdate_si => this%rvars(ir_leafoffdate_si)%r81d, & @@ -1938,8 +1948,29 @@ subroutine get_restart_vectors(this, nc, nsites, sites) end do sites(s)%old_stock = rio_old_stock_si(io_idx_si) - sites(s)%status = rio_cd_status_si(io_idx_si) - sites(s)%dstatus = rio_dd_status_si(io_idx_si) + + ! Site level phenology status flags + if(rio_cd_status_si(io_idx_si) .eq. itrue) then + sites(s)%is_cold = .true. + elseif(rio_cd_status_si(io_idx_si) .eq. ifalse) then + sites(s)%is_cold = .false. + else + write(fates_log(),*) 'An invalid site level cold stress status was found' + write(fates_log(),*) 'io_idx_si = ',io_idx_si + write(fates_log(),*) 'rio_cd_status_si(io_idx_si) = ',rio_cd_status_si(io_idx_si) + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + if(rio_dd_status_si(io_idx_si) .eq. itrue)then + sites(s)%is_drought = .true. + elseif(rio_dd_status_si(io_idx_si) .eq. ifalse) then + sites(s)%is_drought = .false. + else + write(fates_log(),*) 'An invalid site level drought stress status was found' + write(fates_log(),*) 'io_idx_si = ',io_idx_si + write(fates_log(),*) 'rio_dd_status_si(io_idx_si) = ',rio_dd_status_si(io_idx_si) + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + sites(s)%ncd = rio_nchill_days_si(io_idx_si) sites(s)%leafondate = rio_leafondate_si(io_idx_si) sites(s)%leafoffdate = rio_leafoffdate_si(io_idx_si) From 4d2250baf9d89662adccbf7b2113b577b4553f0c Mon Sep 17 00:00:00 2001 From: chonggang xu Date: Wed, 19 Dec 2018 15:55:18 -0800 Subject: [PATCH 14/34] Added drought amplified senescent leaf turnover --- main/EDMainMod.F90 | 2 +- main/EDPftvarcon.F90 | 41 +++++++++++++++++++++++- parameter_files/fates_params_14pfts.cdl | 6 ++++ parameter_files/fates_params_default.cdl | 5 +++ parteh/PRTLossFluxesMod.F90 | 23 +++++++++---- 5 files changed, 69 insertions(+), 8 deletions(-) diff --git a/main/EDMainMod.F90 b/main/EDMainMod.F90 index e8ae1005ca..bda0d0dabc 100644 --- a/main/EDMainMod.F90 +++ b/main/EDMainMod.F90 @@ -340,7 +340,7 @@ subroutine ed_integrate_state_variables(currentSite, bc_in ) ! Conducte Maintenance Turnover (parteh) call currentCohort%prt%CheckMassConservation(ft,3) - call PRTMaintTurnover(currentCohort%prt,ft) + call PRTMaintTurnover(currentCohort%prt,ft,currentSite%is_drought) call currentCohort%prt%CheckMassConservation(ft,4) ! Conduct Growth (parteh) diff --git a/main/EDPftvarcon.F90 b/main/EDPftvarcon.F90 index 1c3b806d57..7a84f5de8d 100644 --- a/main/EDPftvarcon.F90 +++ b/main/EDPftvarcon.F90 @@ -203,7 +203,9 @@ module EDPftvarcon real(r8), allocatable :: phenflush_fraction(:) ! Maximum fraction of storage carbon used to flush leaves ! on bud-burst [kgC/kgC] - + real(r8), allocatable :: senleaf_long_fdrought(:) ! Multiplication factor for leaf longevity of senescent + ! leaves during drought( 1.0 indicates no change) + real(r8), allocatable :: root_long(:) ! root turnover time (longevity) (pft) [yr] real(r8), allocatable :: branch_turnover(:) ! Turnover time for branchfall on live trees (pft) [yr] real(r8), allocatable :: turnover_retrans_mode(:) ! Retranslocation method (pft) @@ -394,6 +396,10 @@ subroutine Register_PFT(this, fates_params) call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) + name = 'fates_senleaf_long_fdrought' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) + name = 'fates_root_long' call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) @@ -822,6 +828,10 @@ subroutine Receive_PFT(this, fates_params) call fates_params%RetreiveParameterAllocate(name=name, & data=this%BB_slope) + name = 'fates_senleaf_long_fdrought' + call fates_params%RetreiveParameterAllocate(name=name, & + data=this%senleaf_long_fdrought) + name = 'fates_root_long' call fates_params%RetreiveParameterAllocate(name=name, & data=this%root_long) @@ -1716,6 +1726,7 @@ subroutine FatesReportPFTParams(is_master) write(fates_log(),fmt0) 'seed_rain = ',EDPftvarcon_inst%seed_rain write(fates_log(),fmt0) 'BB_slope = ',EDPftvarcon_inst%BB_slope write(fates_log(),fmt0) 'root_long = ',EDPftvarcon_inst%root_long + write(fates_log(),fmt0) 'senleaf_long_fdrought = ',EDPftvarcon_inst%senleaf_long_fdrought write(fates_log(),fmt0) 'seed_alloc_mature = ',EDPftvarcon_inst%seed_alloc_mature write(fates_log(),fmt0) 'seed_alloc = ',EDPftvarcon_inst%seed_alloc write(fates_log(),fmt0) 'woody = ',EDPftvarcon_inst%woody @@ -2315,8 +2326,36 @@ subroutine FatesCheckParams(is_master, parteh_mode) call endrun(msg=errMsg(sourcefile, __LINE__)) end if end if + end do + ! Check the turnover rates on the senescing leaf pool + if ( EDPftvarcon_inst%leaf_long(ipft,nleafage)>nearzero ) then + + ! Check that leaf turnover doesn't exeed 1 day + if ( (years_per_day / & + (EDPftvarcon_inst%leaf_long(ipft,nleafage) * & + EDPftvarcon_inst%senleaf_long_fdrought(ipft))) > 1._r8 ) then + write(fates_log(),*) 'Drought-senescent turnover time-scale is greater than 1 day!' + write(fates_log(),*) 'ipft: ',ipft + write(fates_log(),*) 'leaf_long(ipft,nleafage)*senleaf_long_fdrought: ', & + EDPftvarcon_inst%leaf_long(ipft,nleafage)*EDPftvarcon_inst%senleaf_long_fdrought(ipft),' [years]' + write(fates_log(),*) 'Aborting' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + end if + + if ( EDPftvarcon_inst%senleaf_long_fdrought(ipft)1._r8 ) then + write(fates_log(),*) 'senleaf_long_fdrought(ipft) must be greater than 0 ' + write(fates_log(),*) 'or less than or equal to 1.' + write(fates_log(),*) 'Set this to 1 if you want no accelerated senescence turnover' + write(fates_log(),*) 'ipft = ',ipft + write(fates_log(),*) 'senleaf_long_fdrought(ipft) = ',EDPftvarcon_inst%senleaf_long_fdrought(ipft) + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + if ( EDPftvarcon_inst%root_long(ipft)>nearzero ) then ! Check that root turnover doesn't exeed 1 day diff --git a/parameter_files/fates_params_14pfts.cdl b/parameter_files/fates_params_14pfts.cdl index 2a9396da80..a6b58fa7e7 100644 --- a/parameter_files/fates_params_14pfts.cdl +++ b/parameter_files/fates_params_14pfts.cdl @@ -462,6 +462,9 @@ variables: float fates_root_long(fates_pft) ; fates_root_long:units = "yr" ; fates_root_long:long_name = "root longevity (alternatively, turnover time)" ; + float fates_senleaf_long_fdrought(fates_pft) ; + fates_senleaf_long_fdrought:units = "unitless[0-1]" ; + fates_senleaf_long_fdrought:long_name = "multiplication factor for leaf longevity of senescent leaves during drought" ; float fates_roota_par(fates_pft) ; fates_roota_par:units = "1/m" ; fates_roota_par:long_name = "CLM rooting distribution parameter" ; @@ -1091,6 +1094,9 @@ data: fates_rhosvis = 0.16, 0.16, 0.16, 0.16, 0.16, 0.16, 0.16, 0.16, 0.16, 0.16, 0.16, 0.31, 0.31, 0.31 ; + fates_senleaf_long_fdrought = 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, + 1.0, 1.0, 1.0, 1.0, 1.0 ; + fates_root_long = 1, 2, 3, 1, 1.5, 1, 1, 1, 1.5, 1, 1, 1, 1, 1 ; fates_roota_par = 7, 7, 7, 7, 7, 6, 6, 6, 7, 7, 7, 11, 11, 11 ; diff --git a/parameter_files/fates_params_default.cdl b/parameter_files/fates_params_default.cdl index f0917e2e84..8ec2b6257c 100644 --- a/parameter_files/fates_params_default.cdl +++ b/parameter_files/fates_params_default.cdl @@ -405,6 +405,9 @@ variables: float fates_root_long(fates_pft) ; fates_root_long:units = "yr" ; fates_root_long:long_name = "root longevity (alternatively, turnover time)" ; + float fates_senleaf_long_fdrought(fates_pft) ; + fates_senleaf_long_fdrought:units = "unitless[0-1]" ; + fates_senleaf_long_fdrought:long_name = "multiplication factor for leaf longevity of senescent leaves during drought" ; float fates_roota_par(fates_pft) ; fates_roota_par:units = "1/m" ; fates_roota_par:long_name = "CLM rooting distribution parameter" ; @@ -961,6 +964,8 @@ data: fates_rhosvis = 0.16, 0.16 ; + fates_senleaf_long_fdrought = 1.0, 1.0 ; + fates_root_long = 1, 1 ; fates_roota_par = 7, 7 ; diff --git a/parteh/PRTLossFluxesMod.F90 b/parteh/PRTLossFluxesMod.F90 index 1720e5a1a6..e3b4a4b87a 100644 --- a/parteh/PRTLossFluxesMod.F90 +++ b/parteh/PRTLossFluxesMod.F90 @@ -525,7 +525,7 @@ end subroutine DeciduousTurnoverSimpleRetranslocation ! ==================================================================================== - subroutine PRTMaintTurnover(prt,ipft) + subroutine PRTMaintTurnover(prt,ipft,is_drought) ! --------------------------------------------------------------------------------- ! Generic subroutine (wrapper) calling specialized routines handling @@ -533,9 +533,11 @@ subroutine PRTMaintTurnover(prt,ipft) ! --------------------------------------------------------------------------------- class(prt_vartypes) :: prt integer,intent(in) :: ipft + logical,intent(in) :: is_drought ! Is this plant/cohort operating in a drought + ! stress context? if ( int(EDPftvarcon_inst%turnover_retrans_mode(ipft)) == 1 ) then - call MaintTurnoverSimpleRetranslocation(prt,ipft) + call MaintTurnoverSimpleRetranslocation(prt,ipft,is_drought) else write(fates_log(),*) 'A maintenance/retranslocation mode was specified' write(fates_log(),*) 'that is unknown.' @@ -549,7 +551,7 @@ end subroutine PRTMaintTurnover ! =================================================================================== - subroutine MaintTurnoverSimpleRetranslocation(prt,ipft) + subroutine MaintTurnoverSimpleRetranslocation(prt,ipft,is_drought) ! --------------------------------------------------------------------------------- ! This subroutine removes biomass from all applicable pools due to @@ -570,6 +572,8 @@ subroutine MaintTurnoverSimpleRetranslocation(prt,ipft) class(prt_vartypes) :: prt integer, intent(in) :: ipft + logical, intent(in) :: is_drought ! Is this plant/cohort operating in a drought + ! stress context? integer :: i_var ! the variable index integer :: element_id ! the element associated w/ each variable @@ -625,13 +629,20 @@ subroutine MaintTurnoverSimpleRetranslocation(prt,ipft) ! timescale for the senescent pool. aclass_sen_id = size(EDPftvarcon_inst%leaf_long(ipft,:)) - ! Only EVERGREENS HAVE MAINTENANCE LEAF TURNOVER ? + ! Only evergreens have maintenance turnover (must also change trimming logic + ! if we want to change this) ! ------------------------------------------------------------------------------------- if ( (EDPftvarcon_inst%leaf_long(ipft,aclass_sen_id) > nearzero ) .and. & (EDPftvarcon_inst%evergreen(ipft) == itrue) ) then - base_turnover(leaf_organ) = hlm_freq_day / & - EDPftvarcon_inst%leaf_long(ipft,aclass_sen_id) + if(is_drought) then + base_turnover(leaf_organ) = hlm_freq_day / & + (EDPftvarcon_inst%leaf_long(ipft,aclass_sen_id) * & + EDPftvarcon_inst%senleaf_long_fdrought(ipft) ) + else + base_turnover(leaf_organ) = hlm_freq_day / & + EDPftvarcon_inst%leaf_long(ipft,aclass_sen_id) + end if else base_turnover(leaf_organ) = 0.0_r8 endif From d2b45db0b817ad6666d3499239a786bde37012a4 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Mon, 24 Dec 2018 13:07:29 -0700 Subject: [PATCH 15/34] Added provisions for muliple leaf ages in loss fluxes. Fixed bug that had orrured in nutrient enabled flushing and turnover --- parteh/PRTLossFluxesMod.F90 | 86 +++++++++++++++++++++++++++---------- 1 file changed, 64 insertions(+), 22 deletions(-) diff --git a/parteh/PRTLossFluxesMod.F90 b/parteh/PRTLossFluxesMod.F90 index e3b4a4b87a..a47b613862 100644 --- a/parteh/PRTLossFluxesMod.F90 +++ b/parteh/PRTLossFluxesMod.F90 @@ -86,9 +86,15 @@ subroutine PRTPhenologyFlush(prt, ipft, organ_id, c_store_transfer_frac) integer :: i_var_of_organ ! index for all variables in ! a given organ (mostly likely ! synonymous with diff elements) - integer :: i_cvar ! carbon variable index + integer :: i_cvar ! carbon variable index for leaves + ! or other potential organ of interest integer :: i_pos ! spatial position index integer :: i_store ! storage variable index + integer :: i_leaf_pos ! Flush carbon into a specific + ! leaf pool (probably 1st?) + integer :: i_store_pos ! position index for net allocation + ! from retranslocatoin in/out + ! of storage integer :: element_id ! global element identifier real(r8) :: mass_transfer ! The actual mass ! removed from storage @@ -110,6 +116,20 @@ subroutine PRTPhenologyFlush(prt, ipft, organ_id, c_store_transfer_frac) call endrun(msg=errMsg(__FILE__, __LINE__)) end if + if(prt_global%hyp_id .le. 2) then + i_leaf_pos = 1 + i_store_pos = 1 ! hypothesis 1/2 only have + ! 1 storage pool + else + write(fates_log(),*) 'You picked a hypothesis that has not defined' + write(fates_log(),*) ' how and where flushing interacts' + write(fates_log(),*) ' with the storage pool. specifically, ' + write(fates_log(),*) ' if this hypothesis has multiple storage pools' + write(fates_log(),*) ' to pull carbon/resources from' + write(fates_log(),*) 'Exiting' + call endrun(msg=errMsg(__FILE__, __LINE__)) + end if + associate(organ_map => prt_global%organ_map) @@ -144,12 +164,13 @@ subroutine PRTPhenologyFlush(prt, ipft, organ_id, c_store_transfer_frac) ! Get the variable id of the storage pool for this element (carbon12) i_store = prt_global%sp_organ_map(store_organ,element_id) - - ! Loop over all of the coordinate ids - do i_pos = 1,prt_global%state_descriptor(i_var)%num_pos + + + do i_pos = 1,i_leaf_pos ! Calculate the mass transferred out of storage into the pool of interest - mass_transfer = prt%variables(i_store)%val(i_pos) * c_store_transfer_frac + mass_transfer = prt%variables(i_store)%val(i_store_pos) * & + c_store_transfer_frac ! Increment the c pool of interest's allocation flux prt%variables(i_var)%net_alloc(i_pos) = & @@ -161,17 +182,25 @@ subroutine PRTPhenologyFlush(prt, ipft, organ_id, c_store_transfer_frac) ! Increment the storage pool's allocation flux prt%variables(i_store)%net_alloc(i_pos) = & - prt%variables(i_store)%net_alloc(i_pos) - mass_transfer + prt%variables(i_store)%net_alloc(i_store_pos) - mass_transfer ! Update the storage c pool prt%variables(i_store)%val(i_pos) = & - prt%variables(i_store)%val(i_pos) - mass_transfer + prt%variables(i_store)%val(i_store_pos) - mass_transfer end do end if end do - + + + ! This is the variable index for leaf carbon + ! used to calculate the targets for nutrient flushing + i_cvar = prt_global%sp_organ_map(organ_id,carbon12_element) + if(i_cvar < 1) then + write(fates_log(),*) 'Could not determine the carbon var id during flushing' + call endrun(msg=errMsg(__FILE__, __LINE__)) + end if ! Transfer in other elements (nutrients) ! -------------------------------------------------------------------------------- @@ -203,9 +232,8 @@ subroutine PRTPhenologyFlush(prt, ipft, organ_id, c_store_transfer_frac) call endrun(msg=errMsg(__FILE__, __LINE__)) end if - ! Loop over all of the coordinate ids - do i_pos = 1,prt_global%state_descriptor(i_var)%num_pos + do i_pos = 1,i_leaf_pos ! The target quanitity for this element is based on the amount ! of carbon @@ -214,23 +242,23 @@ subroutine PRTPhenologyFlush(prt, ipft, organ_id, c_store_transfer_frac) sp_demand = max(0.0_r8,sp_target - prt%variables(i_var)%val(i_pos)) ! Assume that all of the storage is transferrable - mass_transfer = min(sp_demand, prt%variables(i_store)%val(i_pos)) + mass_transfer = min(sp_demand, prt%variables(i_store)%val(i_store_pos)) ! Increment the pool of interest prt%variables(i_var)%net_alloc(i_pos) = & prt%variables(i_var)%net_alloc(i_pos) + mass_transfer - ! Update the c pool + ! Update the pool prt%variables(i_var)%val(i_pos) = & prt%variables(i_var)%val(i_pos) + mass_transfer - ! Increment the c pool of interest - prt%variables(i_store)%net_alloc(i_pos) = & - prt%variables(i_store)%net_alloc(i_pos) - mass_transfer + ! Increment the store pool allocation diagnostic + prt%variables(i_store)%net_alloc(i_store_pos) = & + prt%variables(i_store)%net_alloc(i_store_pos) - mass_transfer - ! Update the c pool - prt%variables(i_store)%val(i_pos) = & - prt%variables(i_store)%val(i_pos) - mass_transfer + ! Update the store pool + prt%variables(i_store)%val(i_store_pos) = & + prt%variables(i_store)%val(i_store_pos) - mass_transfer end do @@ -441,6 +469,7 @@ subroutine DeciduousTurnoverSimpleRetranslocation(prt,ipft,organ_id,mass_fractio integer :: i_var_of_organ ! loop counter for all element in this organ integer :: element_id ! Element id of the turnover pool integer :: store_var_id ! Variable id of the storage pool + integer :: i_store_pos ! Position index for storage integer :: i_pos ! position index (spatial) real(r8) :: retrans ! retranslocated fraction real(r8) :: turnover_mass ! mass sent to turnover (leaves the plant) @@ -461,6 +490,19 @@ subroutine DeciduousTurnoverSimpleRetranslocation(prt,ipft,organ_id,mass_fractio end if + if(prt_global%hyp_id .le. 2) then + i_store_pos = 1 ! hypothesis 1/2 only have + ! 1 storage pool + else + write(fates_log(),*) 'You picked a hypothesis that has not defined' + write(fates_log(),*) ' how and where flushing interacts' + write(fates_log(),*) ' with the storage pool. specifically, ' + write(fates_log(),*) ' if this hypothesis has multiple storage pools' + write(fates_log(),*) ' to pull carbon/resources from' + write(fates_log(),*) 'Exiting' + call endrun(msg=errMsg(__FILE__, __LINE__)) + end if + do i_var_of_organ = 1, organ_map(organ_id)%num_vars i_var = organ_map(organ_id)%var_id(i_var_of_organ) @@ -508,11 +550,11 @@ subroutine DeciduousTurnoverSimpleRetranslocation(prt,ipft,organ_id,mass_fractio ! Now, since re-translocation is handled by the storage pool, ! we add the re-translocated mass to it - prt%variables(store_var_id)%net_alloc(i_pos) = & - prt%variables(store_var_id)%net_alloc(i_pos) + retranslocated_mass + prt%variables(store_var_id)%net_alloc(i_store_pos) = & + prt%variables(store_var_id)%net_alloc(i_store_pos) + retranslocated_mass - prt%variables(store_var_id)%val(i_pos) = & - prt%variables(store_var_id)%val(i_pos) + retranslocated_mass + prt%variables(store_var_id)%val(i_store_pos) = & + prt%variables(store_var_id)%val(i_store_pos) + retranslocated_mass end do From 49fee562b8d5e88253b4d57d15102a7edd46b8c2 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Mon, 24 Dec 2018 13:09:06 -0700 Subject: [PATCH 16/34] Enabled 4 leaf ages with same vcmax in the 14 pft file. --- parameter_files/fates_params_14pfts.cdl | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/parameter_files/fates_params_14pfts.cdl b/parameter_files/fates_params_14pfts.cdl index e6a6f6b5da..cafbcd4762 100644 --- a/parameter_files/fates_params_14pfts.cdl +++ b/parameter_files/fates_params_14pfts.cdl @@ -5,7 +5,7 @@ dimensions: fates_history_size_bins = 13 ; fates_history_height_bins = 6 ; fates_hydr_organs = 4 ; - fates_leafage_class = 1; + fates_leafage_class = 4; fates_NCWD = 4 ; fates_litterclass = 6 ; fates_scalar = 1 ; @@ -995,7 +995,10 @@ data: fates_leaf_jmaxse = 495, 495, 495, 495, 495, 495, 495, 495, 495, 495, 495, 495, 495, 495 ; - fates_leaf_long = 1.5, 4, 6, 1, 1.5, 1, 1, 1, 1.5, 1, 1, 1, 1, 1 ; + fates_leaf_long = 0.375, 1, 1.5, 0.25, 0.375, 0.25, 0.25, 0.25, 0.375, 0.25, 0.25, 0.25, 0.25, 0.25, + 0.375, 1, 1.5, 0.25, 0.375, 0.25, 0.25, 0.25, 0.375, 0.25, 0.25, 0.25, 0.25, 0.25, + 0.375, 1, 1.5, 0.25, 0.375, 0.25, 0.25, 0.25, 0.375, 0.25, 0.25, 0.25, 0.25, 0.25, + 0.375, 1, 1.5, 0.25, 0.375, 0.25, 0.25, 0.25, 0.375, 0.25, 0.25, 0.25, 0.25, 0.25; fates_leaf_slamax = 0.0954, 0.0954, 0.0954, 0.0954, 0.0954, 0.0954, 0.0954, 0.0954, 0.012, 0.03, 0.03, 0.03, 0.03, 0.03 ; @@ -1015,7 +1018,10 @@ data: fates_leaf_tpuse = 490, 490, 490, 490, 490, 490, 490, 490, 490, 490, 490, 490, 490, 490 ; - fates_leaf_vcmax25top = 50, 65, 63, 39, 62, 41, 58, 58, 62, 54, 54, 78, 78, 78 ; + fates_leaf_vcmax25top = 50, 65, 63, 39, 62, 41, 58, 58, 62, 54, 54, 78, 78, 78, + 50, 65, 63, 39, 62, 41, 58, 58, 62, 54, 54, 78, 78, 78, + 50, 65, 63, 39, 62, 41, 58, 58, 62, 54, 54, 78, 78, 78, + 50, 65, 63, 39, 62, 41, 58, 58, 62, 54, 54, 78, 78, 78 ; fates_leaf_vcmaxha = 65330, 65330, 65330, 65330, 65330, 65330, 65330, 65330, 65330, 65330, 65330, 65330, 65330, 65330 ; From 42393bc4a5dfb585d0853c5b6218658892d8f739 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Mon, 24 Dec 2018 13:10:27 -0700 Subject: [PATCH 17/34] Expanded reporting of carbon conservation failure. Updated cohort and area minimums for termination --- main/EDMainMod.F90 | 41 +++++++++++++++++++++++++++-------------- main/EDTypesMod.F90 | 10 +++++----- 2 files changed, 32 insertions(+), 19 deletions(-) diff --git a/main/EDMainMod.F90 b/main/EDMainMod.F90 index bda0d0dabc..0e2cd52696 100644 --- a/main/EDMainMod.F90 +++ b/main/EDMainMod.F90 @@ -549,6 +549,12 @@ subroutine ed_total_balance_check (currentSite, call_index ) real(r8) :: error ! How much carbon did we gain or lose (should be zero!) real(r8) :: error_frac ! Error as a fraction of total biomass real(r8) :: net_flux ! Difference between recorded fluxes in and out. KgC/site + real(r8) :: leaf_c + real(r8) :: fnrt_c + real(r8) :: sapw_c + real(r8) :: store_c + real(r8) :: struct_c + real(r8) :: repro_c ! nb. There is no time associated with these variables ! because this routine can be called between any two @@ -590,7 +596,7 @@ subroutine ed_total_balance_check (currentSite, call_index ) ! burned_litter * new_patch%area !kG/site/day ! ----------------------------------------------------------------------------------- - if ( error_frac > 10e-6 ) then + if ( error_frac > 10e-6_r8 ) then write(fates_log(),*) 'carbon balance error detected' write(fates_log(),*) 'error fraction relative to biomass stock:',error_frac write(fates_log(),*) 'call index: ',call_index @@ -604,13 +610,17 @@ subroutine ed_total_balance_check (currentSite, call_index ) write(fates_log(),*) 'seeds',seed_stock write(fates_log(),*) 'previous total',currentSite%old_stock - if(debug)then + write(fates_log(),*) 'lat lon',currentSite%lat,currentSite%lon + + ! If this is the first day of simulation, carbon balance reports but does not end the run + if( int(hlm_current_year*10000 + hlm_current_month*100 + hlm_current_day).ne.hlm_reference_date ) then + change_in_stock = 0.0_r8 biomass_stock = 0.0_r8 litter_stock = 0.0_r8 seed_stock = sum(currentSite%seed_bank)*AREA - currentPatch => currentSite%oldest_patch + currentPatch => currentSite%oldest_patch do while(associated(currentPatch)) write(fates_log(),*) '---------------------------------------' write(fates_log(),*) currentPatch%area , sum(currentPatch%cwd_ag), sum(currentPatch%cwd_bg) @@ -618,19 +628,22 @@ subroutine ed_total_balance_check (currentSite, call_index ) write(fates_log(),*)'---' currentCohort => currentPatch%tallest do while(associated(currentCohort)) - write(fates_log(),*) 'structure: ',currentCohort%prt%GetState(struct_organ,all_carbon_elements) - write(fates_log(),*) 'storage: ',currentCohort%prt%GetState(store_organ,all_carbon_elements) + write(fates_log(),*) 'pft: ',currentCohort%pft + write(fates_log(),*) 'dbh: ',currentCohort%dbh + leaf_c = currentCohort%prt%GetState(leaf_organ,all_carbon_elements) + struct_c = currentCohort%prt%GetState(struct_organ,all_carbon_elements) + store_c = currentCohort%prt%GetState(store_organ,all_carbon_elements) + fnrt_c = currentCohort%prt%GetState(fnrt_organ,all_carbon_elements) + repro_c = currentCohort%prt%GetState(repro_organ,all_carbon_elements) + sapw_c = currentCohort%prt%GetState(sapw_organ,all_carbon_elements) + + write(fates_log(),*) 'lc: ',leaf_c,' dc: ',struct_c,' stc: ',store_c + write(fates_log(),*) 'fc: ',fnrt_c,' rc: ',repro_c,' sac: ',sapw_c write(fates_log(),*) 'N plant: ',currentCohort%n - currentCohort => currentCohort%shorter; - enddo !end cohort loop + currentCohort => currentCohort%shorter + enddo !end cohort loop currentPatch => currentPatch%younger - enddo !end patch loop - end if - - write(fates_log(),*) 'lat lon',currentSite%lat,currentSite%lon - - ! If this is the first day of simulation, carbon balance reports but does not end the run - if( int(hlm_current_year*10000 + hlm_current_month*100 + hlm_current_day).ne.hlm_reference_date ) then + enddo !end patch loop write(fates_log(),*) 'aborting on date:',hlm_current_year,hlm_current_month,hlm_current_day call endrun(msg=errMsg(sourcefile, __LINE__)) end if diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 7efc3f5cd9..78c932b576 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -147,11 +147,11 @@ module EDTypesMod integer , parameter :: N_HITE_BINS = 60 ! no. of hite bins used to distribute LAI ! COHORT TERMINATION - real(r8), parameter :: min_npm2 = 1.0E-8_r8 ! minimum cohort number density per m2 before termination - real(r8), parameter :: min_patch_area = 0.001_r8 ! smallest allowable patch area before termination - real(r8), parameter :: min_nppatch = 1.0E-11_r8 ! minimum number of cohorts per patch (min_npm2*min_patch_area) - real(r8), parameter :: min_n_safemath = 1.0E-15_r8 ! in some cases, we want to immediately remove super small - ! number densities of cohorts to prevent FPEs + real(r8), parameter :: min_npm2 = 1.0E-7_r8 ! minimum cohort number density per m2 before termination + real(r8), parameter :: min_patch_area = 0.01_r8 ! smallest allowable patch area before termination + real(r8), parameter :: min_nppatch = min_npm2*min_patch_area ! minimum number of cohorts per patch (min_npm2*min_patch_area) + real(r8), parameter :: min_n_safemath = 1.0E-12_r8 ! in some cases, we want to immediately remove super small + ! number densities of cohorts to prevent FPEs character*4 yearchar From baaee64da6b733db68977d969c26c63fc11703dd Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 27 Dec 2018 10:22:28 -0700 Subject: [PATCH 18/34] Reverting patch minimums, will change in different set --- main/EDTypesMod.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 78c932b576..46f039508e 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -147,10 +147,10 @@ module EDTypesMod integer , parameter :: N_HITE_BINS = 60 ! no. of hite bins used to distribute LAI ! COHORT TERMINATION - real(r8), parameter :: min_npm2 = 1.0E-7_r8 ! minimum cohort number density per m2 before termination - real(r8), parameter :: min_patch_area = 0.01_r8 ! smallest allowable patch area before termination + real(r8), parameter :: min_npm2 = 1.0E-8_r8 ! minimum cohort number density per m2 before termination + real(r8), parameter :: min_patch_area = 0.001_r8 ! smallest allowable patch area before termination real(r8), parameter :: min_nppatch = min_npm2*min_patch_area ! minimum number of cohorts per patch (min_npm2*min_patch_area) - real(r8), parameter :: min_n_safemath = 1.0E-12_r8 ! in some cases, we want to immediately remove super small + real(r8), parameter :: min_n_safemath = 1.0E-15_r8 ! in some cases, we want to immediately remove super small ! number densities of cohorts to prevent FPEs character*4 yearchar From 2db8c6abf87457ae083cf5c1499d37b4a1884e8a Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Sun, 16 Dec 2018 16:27:27 -0800 Subject: [PATCH 19/34] Added a leaf layer, reduced cohorts --- main/EDTypesMod.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 16bd7db647..518ab7924b 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -16,9 +16,9 @@ module EDTypesMod save integer, parameter :: maxPatchesPerSite = 10 ! maximum number of patches to live on a site - integer, parameter :: maxCohortsPerPatch = 160 ! maximum number of cohorts per patch + integer, parameter :: maxCohortsPerPatch = 40 ! maximum number of cohorts per patch - integer, parameter :: nclmax = 2 ! Maximum number of canopy layers + integer, parameter :: nclmax = 3 ! Maximum number of canopy layers integer, parameter :: ican_upper = 1 ! Nominal index for the upper canopy integer, parameter :: ican_ustory = 2 ! Nominal index for understory in two-canopy system From 1bebe731c6322c9325204e5165a86cb4916aa948 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Fri, 4 Jan 2019 12:04:41 -0700 Subject: [PATCH 20/34] Changed cohort cap and the fusion tolerances to 80 and 8% respectively --- main/EDTypesMod.F90 | 2 +- parameter_files/fates_params_14pfts.cdl | 2 +- parameter_files/fates_params_coastal_veg.cdl | 2 +- parameter_files/fates_params_default.cdl | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 518ab7924b..7d8266051c 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -16,7 +16,7 @@ module EDTypesMod save integer, parameter :: maxPatchesPerSite = 10 ! maximum number of patches to live on a site - integer, parameter :: maxCohortsPerPatch = 40 ! maximum number of cohorts per patch + integer, parameter :: maxCohortsPerPatch = 80 ! maximum number of cohorts per patch integer, parameter :: nclmax = 3 ! Maximum number of canopy layers integer, parameter :: ican_upper = 1 ! Nominal index for the upper canopy diff --git a/parameter_files/fates_params_14pfts.cdl b/parameter_files/fates_params_14pfts.cdl index cf3d4de880..ab665f3e36 100644 --- a/parameter_files/fates_params_14pfts.cdl +++ b/parameter_files/fates_params_14pfts.cdl @@ -616,7 +616,7 @@ data: fates_canopy_closure_thresh = 0.8 ; - fates_cohort_fusion_tol = 0.05 ; + fates_cohort_fusion_tol = 0.08 ; fates_comp_excln = 3 ; diff --git a/parameter_files/fates_params_coastal_veg.cdl b/parameter_files/fates_params_coastal_veg.cdl index a6b1ac59d3..47f41c811e 100644 --- a/parameter_files/fates_params_coastal_veg.cdl +++ b/parameter_files/fates_params_coastal_veg.cdl @@ -880,7 +880,7 @@ data: fates_canopy_closure_thresh = 0.8 ; - fates_cohort_fusion_tol = 0.05 ; + fates_cohort_fusion_tol = 0.08 ; fates_comp_excln = 3 ; diff --git a/parameter_files/fates_params_default.cdl b/parameter_files/fates_params_default.cdl index bed07aed11..f16971c052 100644 --- a/parameter_files/fates_params_default.cdl +++ b/parameter_files/fates_params_default.cdl @@ -999,7 +999,7 @@ data: fates_canopy_closure_thresh = 0.8 ; - fates_cohort_fusion_tol = 0.05 ; + fates_cohort_fusion_tol = 0.08 ; fates_comp_excln = 3 ; From 5c007050df5234eb54cdee3842a77f45d55f995e Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 8 Jan 2019 12:21:28 -0800 Subject: [PATCH 21/34] Half-way through streamlining equal promotion/demotion of same-height cohorts --- biogeochem/EDCanopyStructureMod.F90 | 165 +++++++++++++--------------- 1 file changed, 77 insertions(+), 88 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index e30d339bb7..ecdd0e9ab9 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -343,6 +343,8 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) real(r8) :: remainder_area_hold real(r8) :: sumweights real(r8) :: sumweights_old + real(r8) :: sumactual ! for rank-ordered same-size cohorts + ! this tallies their excluded area real(r8) :: arealayer ! the area of the current canopy layer integer :: exceedance_counter ! when seeking to rebalance demotion exceedance ! keep a loop counter to check for hangs @@ -375,109 +377,96 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) currentSite%spread,currentCohort%pft,currentCohort%c_area) if( 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 + ! ---------------------------------------------------------- + currentCohort%excl_weight = & currentCohort%n/(currentCohort%hite**ED_val_comp_excln) + sumweights = sumweights + currentCohort%excl_weight + else + + ! ----------------------------------------------------------- ! Rank ordered deterministic method + ! ----------------------------------------------------------- - ! check to make sure there are no cohorts of equal size - tied_size_with_neighbor = .false. - if (associated(currentCohort%shorter)) then - if (currentCohort%shorter%hite .eq. currentCohort%hite ) then - tied_size_with_neighbor = .true. - endif - endif - if (associated(currentCohort%taller)) then - if (currentCohort%taller%hite .eq. currentCohort%hite ) then - tied_size_with_neighbor = .true. + ! If there cohorts that have the exact same height (which is possible, really) + ! we don't want to unilaterally promote/demote one before the others. + ! So we promote them (as a unit) + + ! now we need to go through and figure out how many equal-size cohorts there are. + ! then we need to go through, add up the collective crown areas of all equal-sized + ! and equal-canopy-layer cohorts, + ! and then demote from each as if they were a single group + ! + total_crownarea_of_tied_cohorts = currentCohort%c_area + + tied_size_with_neighbors = .false. + nextc => currentCohort%taller + do while (associated(nextc)) then + if ( abs(nextc%hite - currentCohort%hite) < nearzero ) then + if( nextc%canopy_layer .eq. currentCohort%canopy_layer ) then + tied_size_with_neighbors = .true. + total_crownarea_of_tied_cohorts = & + total_crownarea_of_tied_cohort + currentCohort%c_area + end if + else + break endif + nextc => nextc%taller endif + + if ( tied_size_with_neighbors ) then + + currentCohort%excl_weight = & + max(0.0_r8,min(currentCohort%c_area, & + (currentCohort%c_area/total_crownarea_of_tied_cohorts) * & + (demote_area - sumweights) )) - if ( tied_size_with_neighbor ) then - if ( DEBUG ) then - write(fates_log(),*) 'tied_size_with_neighbor eq true in demotion phase' - endif - ! now we need to go through and figure out how many equal-size cohorts there are. - ! then we need to go through, add up the collective crown areas of all equal-sized and equal-canopy-layer cohorts, - ! and then demote from each as if they were a single group - ! - total_crownarea_of_tied_cohorts = currentCohort%c_area - ! - ! first the "shorter" cohorts (scare-quotes because they aren't actually shorter) - found_shortest_equal_neighbor = .false. - cohort_tosearch_relative_to => currentCohort - whileloop_counter = 0 - do while ( .not. found_shortest_equal_neighbor) - whileloop_counter = whileloop_counter + 1 - if (associated(cohort_tosearch_relative_to%shorter)) then - cohort_tocompare_to => cohort_tosearch_relative_to%shorter - if (cohort_tocompare_to%hite .eq. currentCohort%hite ) then - if (cohort_tocompare_to%canopy_layer .eq. currentCohort%canopy_layer ) then - total_crownarea_of_tied_cohorts = total_crownarea_of_tied_cohorts + cohort_tocompare_to%c_area - endif - cohort_tosearch_relative_to => cohort_tocompare_to - else - found_shortest_equal_neighbor = .true. - end if - else - found_shortest_equal_neighbor = .true. - endif - if ( whileloop_counter .ge. maxCohortsPerPatch ) then - ! something has gone horribly wrong and we are in an infite loop. - call endrun(msg=errMsg(sourcefile, __LINE__)) - endif - end do - ! - ! then the "taller" cohorts (scare-quotes because they aren't actually taller) - has_taller_equalsized_neighbor = .false. ! init this as false - found_tallest_equal_neighbor = .false. - cohort_tosearch_relative_to => currentCohort - whileloop_counter = 0 - do while ( .not. found_tallest_equal_neighbor) - whileloop_counter = whileloop_counter + 1 - if (associated(cohort_tosearch_relative_to%taller)) then - cohort_tocompare_to => cohort_tosearch_relative_to%taller - if (cohort_tocompare_to%hite .eq. currentCohort%hite ) then - if (cohort_tocompare_to%canopy_layer .eq. currentCohort%canopy_layer ) then - total_crownarea_of_tied_cohorts = total_crownarea_of_tied_cohorts + cohort_tocompare_to%c_area - has_taller_equalsized_neighbor = .true. - endif - cohort_tosearch_relative_to => cohort_tocompare_to - else - found_tallest_equal_neighbor = .true. + sumactual = currentCohort%excl_weight + + + nextc => currentCohort%taller + do while (associated(nextc)) then + if ( abs(nextc%hite - currentCohort%hite) < nearzero ) then + if (nextc%canopy_layer .eq. currentCohort%canopy_layer ) then + ! now we know the total crown area of all equal-sized, + ! equal-canopy-layer cohorts + nextc%excl_weight = & + max(0.0_r8,min(nextc%c_area, & + (nextc%c_area/total_crownarea_of_tied_cohorts) * & + (demote_area - sumweights) )) + sumactual = sumactual + currentCohort%excl_weight end if else - found_tallest_equal_neighbor = .true. - endif - if ( whileloop_counter .ge. maxCohortsPerPatch ) then - ! something has gone horribly wrong and we are in an infite loop. - call endrun(msg=errMsg(sourcefile, __LINE__)) + break endif - end do - ! - ! now we know the total crown area of all equal-sized, equal-canopy-layer cohorts - currentCohort%excl_weight = & - max(min(currentCohort%c_area, (currentCohort%c_area/total_crownarea_of_tied_cohorts) * & - (demote_area - sumweights) ), 0._r8) - else ! i.e. tied_size_with_neighbor = .false. + nextc => nextc%taller + endif + + ! Update the current cohort pointer to the last similar cohort + ! Its ok if this is not in the right layer + if(associated(nextc))then + currentCohort => nextc%shorter + else + currentCohort => currentPatch%tallest + end if + sumweights = sumweights + sumactual + + else + currentCohort%excl_weight = & max(min(currentCohort%c_area, demote_area - sumweights ), 0._r8) - endif - endif - ! - ! when two or more cohorts have the same size, we need to keep track of their cumulative demoted crown area - ! in a separate buffer and add it once we reach the last of the equal-sized cohorts - if ((ED_val_comp_excln .lt. 0.0_r8) .and. tied_size_with_neighbor .and. & - has_taller_equalsized_neighbor) then - sumweights_equalsizebuffer = sumweights_equalsizebuffer + currentCohort%excl_weight - else if ( (ED_val_comp_excln .lt. 0.0_r8) .and. tied_size_with_neighbor) then - sumweights = sumweights + currentCohort%excl_weight + sumweights_equalsizebuffer - sumweights_equalsizebuffer = 0._r8 - else - sumweights = sumweights + currentCohort%excl_weight + sumweights = sumweights + currentCohort%excl_weight + + end if + + endif endif currentCohort => currentCohort%taller @@ -856,7 +845,7 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) endif if ( tied_size_with_neighbor ) then - if ( DEBUG ) then + if ( debug ) then write(fates_log(),*) 'tied_size_with_neighbor eq true in promotion phase' endif From b07be09b0aa64066e668bfe12d6b0f9a40c29296 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 8 Jan 2019 14:54:11 -0800 Subject: [PATCH 22/34] changed fusion algorithm to have forward only same-height evaluation, should also work for more than two cohorts with nigh similar height --- biogeochem/EDCanopyStructureMod.F90 | 218 ++++++++++++---------------- 1 file changed, 93 insertions(+), 125 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index ecdd0e9ab9..120286cde3 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -59,6 +59,9 @@ module EDCanopyStructureMod real(r8), parameter :: area_target_precision = 1.0E-11_r8 ! Area conservation must be within this tolerance real(r8), parameter :: area_check_precision = 1.0E-9_r8 ! Area conservation checks must be within this tolerance + real(r8), parameter :: similar_height_tol = 1.0E-3_r8 ! I think trees that differ by 1mm + ! can be roughly considered the same right? + ! 10/30/09: Created by Rosie Fisher ! 2017/2018: Modifications and updates by Ryan Knox @@ -328,7 +331,9 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) integer, intent(in) :: i_lyr ! Current canopy layer of interest ! !LOCAL VARIABLES: - type(ed_cohort_type), pointer :: currentCohort,copyc + type(ed_cohort_type), pointer :: currentCohort + type(ed_cohort_type), pointer :: copyc + type(ed_cohort_type), pointer :: nextc ! The next cohort in line integer :: i_cwd ! Index for CWD pool real(r8) :: cc_loss ! cohort crown area loss in demotion (m2) real(r8) :: leaf_c ! leaf carbon [kg] @@ -343,18 +348,14 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) real(r8) :: remainder_area_hold real(r8) :: sumweights real(r8) :: sumweights_old - real(r8) :: sumactual ! for rank-ordered same-size cohorts + real(r8) :: sumequal ! for rank-ordered same-size cohorts ! this tallies their excluded area real(r8) :: arealayer ! the area of the current canopy layer integer :: exceedance_counter ! when seeking to rebalance demotion exceedance ! keep a loop counter to check for hangs - logical :: tied_size_with_neighbor, has_taller_equalsized_neighbor - logical :: found_shortest_equal_neighbor, found_tallest_equal_neighbor + logical :: tied_size_with_neighbors type(ed_cohort_type), pointer :: cohort_tosearch_relative_to, cohort_tocompare_to real(r8) :: total_crownarea_of_tied_cohorts - real(r8) :: sumweights_equalsizebuffer - integer :: whileloop_counter - ! First, determine how much total canopy area we have in this layer @@ -369,7 +370,6 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) ! We go in order from shortest to tallest for ranked demotion sumweights = 0.0_r8 - sumweights_equalsizebuffer = 0.0_r8 currentCohort => currentPatch%shortest do while (associated(currentCohort)) @@ -394,32 +394,30 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) ! ----------------------------------------------------------- ! Rank ordered deterministic method ! ----------------------------------------------------------- - ! If there cohorts that have the exact same height (which is possible, really) ! we don't want to unilaterally promote/demote one before the others. - ! So we promote them (as a unit) - + ! So we <>mote them as a unit ! now we need to go through and figure out how many equal-size cohorts there are. ! then we need to go through, add up the collective crown areas of all equal-sized ! and equal-canopy-layer cohorts, ! and then demote from each as if they were a single group - ! + total_crownarea_of_tied_cohorts = currentCohort%c_area tied_size_with_neighbors = .false. nextc => currentCohort%taller - do while (associated(nextc)) then - if ( abs(nextc%hite - currentCohort%hite) < nearzero ) then + do while (associated(nextc)) + if ( abs(nextc%hite - currentCohort%hite) < similar_height_tol ) then if( nextc%canopy_layer .eq. currentCohort%canopy_layer ) then tied_size_with_neighbors = .true. total_crownarea_of_tied_cohorts = & - total_crownarea_of_tied_cohort + currentCohort%c_area + total_crownarea_of_tied_cohorts + currentCohort%c_area end if else - break + exit endif nextc => nextc%taller - endif + end do if ( tied_size_with_neighbors ) then @@ -428,12 +426,11 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) (currentCohort%c_area/total_crownarea_of_tied_cohorts) * & (demote_area - sumweights) )) - sumactual = currentCohort%excl_weight - - + sumequal = currentCohort%excl_weight + nextc => currentCohort%taller - do while (associated(nextc)) then - if ( abs(nextc%hite - currentCohort%hite) < nearzero ) then + do while (associated(nextc)) + if ( abs(nextc%hite - currentCohort%hite) < similar_height_tol ) then if (nextc%canopy_layer .eq. currentCohort%canopy_layer ) then ! now we know the total crown area of all equal-sized, ! equal-canopy-layer cohorts @@ -441,13 +438,13 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) max(0.0_r8,min(nextc%c_area, & (nextc%c_area/total_crownarea_of_tied_cohorts) * & (demote_area - sumweights) )) - sumactual = sumactual + currentCohort%excl_weight + sumequal = sumequal + currentCohort%excl_weight end if else - break + exit endif nextc => nextc%taller - endif + end do ! Update the current cohort pointer to the last similar cohort ! Its ok if this is not in the right layer @@ -456,17 +453,14 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) else currentCohort => currentPatch%tallest end if - sumweights = sumweights + sumactual + sumweights = sumweights + sumequal else - currentCohort%excl_weight = & max(min(currentCohort%c_area, demote_area - sumweights ), 0._r8) sumweights = sumweights + currentCohort%excl_weight - end if - endif endif currentCohort => currentCohort%taller @@ -742,12 +736,16 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) ! !LOCAL VARIABLES: type(ed_cohort_type), pointer :: currentCohort type(ed_cohort_type), pointer :: copyc + type(ed_cohort_type), pointer :: nextc ! the next cohort, or used for looping + ! cohorts against the current real(r8) :: promote_area real(r8) :: newarea real(r8) :: sumweights real(r8) :: sumweights_old + real(r8) :: sumequal ! for tied cohorts, the sum of weights in + ! their group integer :: exceedance_counter real(r8) :: remainder_area real(r8) :: remainder_area_hold @@ -760,12 +758,9 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) real(r8) :: store_c ! storage carbon [kg] real(r8) :: struct_c ! structure carbon [kg] - logical :: tied_size_with_neighbor, has_shorter_equalsized_neighbor - logical :: found_shortest_equal_neighbor, found_tallest_equal_neighbor + logical :: tied_size_with_neighbors type(ed_cohort_type), pointer :: cohort_tosearch_relative_to, cohort_tocompare_to real(r8) :: total_crownarea_of_tied_cohorts - real(r8) :: sumweights_equalsizebuffer - integer :: whileloop_counter call CanopyLayerArea(currentPatch,currentSite%spread,i_lyr,arealayer_current) call CanopyLayerArea(currentPatch,currentSite%spread,i_lyr+1,arealayer_below) @@ -819,113 +814,86 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) ! This is the opposite of the demotion weighting... sumweights = 0.0_r8 - sumweights_equalsizebuffer = 0.0_r8 currentCohort => currentPatch%tallest do while (associated(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. currentCohort%prom_weight = currentCohort%n*currentCohort%hite**ED_val_comp_excln else + + ! ------------------------------------------------------------------ ! Rank ordered deterministic method - - ! check to make sure there are no cohorts of equal size - tied_size_with_neighbor = .false. - if (associated(currentCohort%shorter)) then - if (currentCohort%shorter%hite .eq. currentCohort%hite ) then - tied_size_with_neighbor = .true. - endif - endif - if (associated(currentCohort%taller)) then - if (currentCohort%taller%hite .eq. currentCohort%hite ) then - tied_size_with_neighbor = .true. - endif - endif - - if ( tied_size_with_neighbor ) then - if ( debug ) then - write(fates_log(),*) 'tied_size_with_neighbor eq true in promotion phase' + ! If there cohorts that have the exact same height (which is possible, really) + ! we don't want to unilaterally promote/demote one before the others. + ! So we <>mote them as a unit + ! now we need to go through and figure out how many equal-size cohorts there are. + ! then we need to go through, add up the collective crown areas of all equal-sized + ! and equal-canopy-layer cohorts, + ! and then demote from each as if they were a single group + ! ------------------------------------------------------------------ + + total_crownarea_of_tied_cohorts = currentCohort%c_area + tied_size_with_neighbors = .false. + nextc => currentCohort%shorter + do while (associated(nextc)) + if ( abs(nextc%hite - currentCohort%hite) < similar_height_tol ) then + if( nextc%canopy_layer .eq. currentCohort%canopy_layer ) then + tied_size_with_neighbors = .true. + total_crownarea_of_tied_cohorts = & + total_crownarea_of_tied_cohorts + currentCohort%c_area + end if + else + exit endif - - ! now we need to go through and figure out how many equal-size cohorts there are. - ! then we need to go through, add up the collective crown areas of all equal-sized and equal-canopy-layer cohorts, - ! and then promote from each as if they were a single group - ! - total_crownarea_of_tied_cohorts = currentCohort%c_area - ! - ! first the "shorter" cohorts (scare-quotes because they aren't actually shorter) - has_shorter_equalsized_neighbor = .false. ! init this as false - found_shortest_equal_neighbor = .false. - cohort_tosearch_relative_to => currentCohort - whileloop_counter = 0 - do while ( .not. found_shortest_equal_neighbor) - whileloop_counter = whileloop_counter + 1 - if (associated(cohort_tosearch_relative_to%shorter)) then - cohort_tocompare_to => cohort_tosearch_relative_to%shorter - if (cohort_tocompare_to%hite .eq. currentCohort%hite ) then - if (cohort_tocompare_to%canopy_layer .eq. currentCohort%canopy_layer ) then - total_crownarea_of_tied_cohorts = total_crownarea_of_tied_cohorts + cohort_tocompare_to%c_area - has_shorter_equalsized_neighbor = .true. - endif - cohort_tosearch_relative_to => cohort_tocompare_to - else - found_shortest_equal_neighbor = .true. + nextc => nextc%shorter + end do + + if ( tied_size_with_neighbors ) then + + currentCohort%prom_weight = & + max(0.0_r8,min(currentCohort%c_area, & + (currentCohort%c_area/total_crownarea_of_tied_cohorts) * & + (promote_area - sumweights) )) + sumequal = currentCohort%prom_weight + + nextc => currentCohort%shorter + do while (associated(nextc)) + if ( abs(nextc%hite - currentCohort%hite) < similar_height_tol ) then + if (nextc%canopy_layer .eq. currentCohort%canopy_layer ) then + ! now we know the total crown area of all equal-sized, + ! equal-canopy-layer cohorts + nextc%prom_weight = & + max(0.0_r8,min(nextc%c_area, & + (nextc%c_area/total_crownarea_of_tied_cohorts) * & + (promote_area - sumweights) )) + sumequal = sumequal + currentCohort%prom_weight end if else - found_shortest_equal_neighbor = .true. - endif - if ( whileloop_counter .ge. maxCohortsPerPatch ) then - ! something has gone horribly wrong and we are in an infite loop. - call endrun(msg=errMsg(sourcefile, __LINE__)) - endif - end do - ! - ! then the "taller" cohorts (scare-quotes because they aren't actually taller) - found_tallest_equal_neighbor = .false. - cohort_tosearch_relative_to => currentCohort - whileloop_counter = 0 - do while ( .not. found_tallest_equal_neighbor) - whileloop_counter = whileloop_counter + 1 - if (associated(cohort_tosearch_relative_to%taller)) then - cohort_tocompare_to => cohort_tosearch_relative_to%taller - if (cohort_tocompare_to%hite .eq. currentCohort%hite ) then - if (cohort_tocompare_to%canopy_layer .eq. currentCohort%canopy_layer ) then - total_crownarea_of_tied_cohorts = total_crownarea_of_tied_cohorts + cohort_tocompare_to%c_area - endif - cohort_tosearch_relative_to => cohort_tocompare_to - else - found_tallest_equal_neighbor = .true. - end if - else - found_tallest_equal_neighbor = .true. - endif - if ( whileloop_counter .ge. maxCohortsPerPatch ) then - ! something has gone horribly wrong and we are in an infite loop. - call endrun(msg=errMsg(sourcefile, __LINE__)) + exit endif + nextc => nextc%shorter end do - ! - ! now we know the total crown area of all equal-sized, equal-canopy-layer cohorts - currentCohort%prom_weight = max(min(currentCohort%c_area, & - (currentCohort%c_area/total_crownarea_of_tied_cohorts) * (promote_area - sumweights) ), 0._r8) - else ! i.e. tied_size_with_neighbor = .false. - currentCohort%prom_weight = max(min(currentCohort%c_area, & - promote_area - sumweights ), 0._r8) - endif - endif - ! - ! when two or more cohorts have the same size, we need to keep track of their cumulative demoted crown area - ! in a separate buffer and add it once we reach the last of the equal-sized cohorts - if ((ED_val_comp_excln .lt. 0.0_r8) .and. tied_size_with_neighbor .and. & - has_shorter_equalsized_neighbor) then - sumweights_equalsizebuffer = sumweights_equalsizebuffer + currentCohort%prom_weight - else if ( (ED_val_comp_excln .lt. 0.0_r8) .and. tied_size_with_neighbor) then - sumweights = sumweights + currentCohort%prom_weight + sumweights_equalsizebuffer - sumweights_equalsizebuffer = 0._r8 - else - sumweights = sumweights + currentCohort%prom_weight + + ! Update the current cohort pointer to the last similar cohort + ! Its ok if this is not in the right layer + if(associated(nextc))then + currentCohort => nextc%taller + else + currentCohort => currentPatch%shortest + end if + sumweights = sumweights + sumequal + + else + currentCohort%prom_weight = & + max(min(currentCohort%c_area, promote_area - sumweights ), 0._r8) + sumweights = sumweights + currentCohort%prom_weight + + end if + endif endif currentCohort => currentCohort%shorter From bacaf1e7e48c30464f4158ad719a1d3899555708 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 8 Jan 2019 16:08:15 -0800 Subject: [PATCH 23/34] small bug fix on rank-ordered shared area tally --- biogeochem/EDCanopyStructureMod.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 120286cde3..791fbafb69 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -411,7 +411,7 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) if( nextc%canopy_layer .eq. currentCohort%canopy_layer ) then tied_size_with_neighbors = .true. total_crownarea_of_tied_cohorts = & - total_crownarea_of_tied_cohorts + currentCohort%c_area + total_crownarea_of_tied_cohorts + nextc%c_area end if else exit @@ -844,7 +844,7 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) if( nextc%canopy_layer .eq. currentCohort%canopy_layer ) then tied_size_with_neighbors = .true. total_crownarea_of_tied_cohorts = & - total_crownarea_of_tied_cohorts + currentCohort%c_area + total_crownarea_of_tied_cohorts + nextc%c_area end if else exit From b2382baf439c5423eba05007253512f981a2d4d7 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 8 Jan 2019 16:55:16 -0800 Subject: [PATCH 24/34] Another bug fix to same-height rank weighting --- biogeochem/EDCanopyStructureMod.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 791fbafb69..dae2b5bf84 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -438,7 +438,7 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) max(0.0_r8,min(nextc%c_area, & (nextc%c_area/total_crownarea_of_tied_cohorts) * & (demote_area - sumweights) )) - sumequal = sumequal + currentCohort%excl_weight + sumequal = sumequal + nextc%excl_weight end if else exit @@ -870,7 +870,7 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) max(0.0_r8,min(nextc%c_area, & (nextc%c_area/total_crownarea_of_tied_cohorts) * & (promote_area - sumweights) )) - sumequal = sumequal + currentCohort%prom_weight + sumequal = sumequal + nextc%prom_weight end if else exit From 2dc56c78f59210f58b4129af0be2f3106b467ce4 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 10 Jan 2019 18:07:35 -0700 Subject: [PATCH 25/34] Forgot to sumweights for non ranked case during promotion --- biogeochem/EDCanopyStructureMod.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index dae2b5bf84..985f5244a8 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -823,6 +823,7 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) if (ED_val_comp_excln .ge. 0.0_r8 ) then ! normal (stochastic) case, as above. currentCohort%prom_weight = currentCohort%n*currentCohort%hite**ED_val_comp_excln + sumweights = sumweights + currentCohort%prom_weight else ! ------------------------------------------------------------------ From 14bceb85c52e6e7d0d909e56d83897ce1532abfb Mon Sep 17 00:00:00 2001 From: Charlie Koven Date: Thu, 10 Jan 2019 17:56:19 -0800 Subject: [PATCH 26/34] Update biogeochem/EDCanopyStructureMod.F90 Co-Authored-By: rgknox --- biogeochem/EDCanopyStructureMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 985f5244a8..9aefc379f7 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -394,7 +394,7 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) ! ----------------------------------------------------------- ! Rank ordered deterministic method ! ----------------------------------------------------------- - ! If there cohorts that have the exact same height (which is possible, really) + ! If there are cohorts that have the exact same height (which is possible, really) ! we don't want to unilaterally promote/demote one before the others. ! So we <>mote them as a unit ! now we need to go through and figure out how many equal-size cohorts there are. From 5a617ea298931f7be5c47029e63d8f0eca55133e Mon Sep 17 00:00:00 2001 From: Charlie Koven Date: Thu, 10 Jan 2019 17:56:30 -0800 Subject: [PATCH 27/34] Update biogeochem/EDCanopyStructureMod.F90 Co-Authored-By: rgknox --- biogeochem/EDCanopyStructureMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 9aefc379f7..8be0045177 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -828,7 +828,7 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) ! ------------------------------------------------------------------ ! Rank ordered deterministic method - ! If there cohorts that have the exact same height (which is possible, really) + ! If there are cohorts that have the exact same height (which is possible, really) ! we don't want to unilaterally promote/demote one before the others. ! So we <>mote them as a unit ! now we need to go through and figure out how many equal-size cohorts there are. From 6ef2a30be6bb468da8ba9b661195b87bd21330b8 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 31 Jan 2019 13:08:41 -0700 Subject: [PATCH 28/34] Reduced the number of default leaf age bins to 2 --- parameter_files/fates_params_14pfts.cdl | 10 +++------- parameter_files/fates_params_default.cdl | 11 ++++------- 2 files changed, 7 insertions(+), 14 deletions(-) diff --git a/parameter_files/fates_params_14pfts.cdl b/parameter_files/fates_params_14pfts.cdl index 442f6f350f..3e53fa8659 100644 --- a/parameter_files/fates_params_14pfts.cdl +++ b/parameter_files/fates_params_14pfts.cdl @@ -5,7 +5,7 @@ dimensions: fates_history_size_bins = 13 ; fates_history_height_bins = 6 ; fates_hydr_organs = 4 ; - fates_leafage_class = 4; + fates_leafage_class = 2; fates_NCWD = 4 ; fates_litterclass = 6 ; fates_scalar = 1 ; @@ -995,10 +995,8 @@ data: fates_leaf_jmaxse = 495, 495, 495, 495, 495, 495, 495, 495, 495, 495, 495, 495, 495, 495 ; - fates_leaf_long = 0.375, 1, 1.5, 0.25, 0.375, 0.25, 0.25, 0.25, 0.375, 0.25, 0.25, 0.25, 0.25, 0.25, - 0.375, 1, 1.5, 0.25, 0.375, 0.25, 0.25, 0.25, 0.375, 0.25, 0.25, 0.25, 0.25, 0.25, - 0.375, 1, 1.5, 0.25, 0.375, 0.25, 0.25, 0.25, 0.375, 0.25, 0.25, 0.25, 0.25, 0.25, - 0.375, 1, 1.5, 0.25, 0.375, 0.25, 0.25, 0.25, 0.375, 0.25, 0.25, 0.25, 0.25, 0.25; + fates_leaf_long = 0.75, 2, 3.0, 0.5, 0.75, 0.5, 0.5, 0.5, 0.75, 0.5, 0.5, 0.5, 0.5, 0.5, + 0.75, 2, 3.0, 0.5, 0.75, 0.5, 0.5, 0.5, 0.75, 0.5, 0.5, 0.5, 0.5, 0.5; fates_leaf_slamax = 0.0954, 0.0954, 0.0954, 0.0954, 0.0954, 0.0954, 0.0954, 0.0954, 0.012, 0.03, 0.03, 0.03, 0.03, 0.03 ; @@ -1019,8 +1017,6 @@ data: 490, 490, 490 ; fates_leaf_vcmax25top = 50, 65, 63, 39, 62, 41, 58, 58, 62, 54, 54, 78, 78, 78, - 50, 65, 63, 39, 62, 41, 58, 58, 62, 54, 54, 78, 78, 78, - 50, 65, 63, 39, 62, 41, 58, 58, 62, 54, 54, 78, 78, 78, 50, 65, 63, 39, 62, 41, 58, 58, 62, 54, 54, 78, 78, 78 ; fates_leaf_vcmaxha = 65330, 65330, 65330, 65330, 65330, 65330, 65330, 65330, diff --git a/parameter_files/fates_params_default.cdl b/parameter_files/fates_params_default.cdl index 21383f067f..923774d361 100644 --- a/parameter_files/fates_params_default.cdl +++ b/parameter_files/fates_params_default.cdl @@ -6,7 +6,7 @@ dimensions: fates_history_size_bins = 13 ; fates_hydr_organs = 4 ; fates_prt_organs = 6 ; - fates_leafage_class = 4; + fates_leafage_class = 2; fates_litterclass = 6 ; fates_pft = 2 ; fates_scalar = 1 ; @@ -880,10 +880,9 @@ data: fates_leaf_jmaxse = 495, 495 ; - fates_leaf_long = 0.375, 0.375, - 0.375, 0.375, - 0.375, 0.375, - 0.375, 0.375 ; + fates_leaf_long = 0.75, 0.75, + 0.75, 0.75 ; + fates_leaf_slamax = 0.0954, 0.0954 ; @@ -898,8 +897,6 @@ data: fates_leaf_tpuse = 490, 490 ; fates_leaf_vcmax25top = 50, 50, - 50, 50, - 50, 50, 50, 50; fates_leaf_vcmaxha = 65330, 65330 ; From 22392399e29904232aabfd85b41e5719cac6e77f Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 31 Jan 2019 12:16:17 -0800 Subject: [PATCH 29/34] Updated the hydro default file to have leaf age bins --- parameter_files/fates_params_hydro_default.cdl | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/parameter_files/fates_params_hydro_default.cdl b/parameter_files/fates_params_hydro_default.cdl index e6e3adfbec..59106459e0 100644 --- a/parameter_files/fates_params_hydro_default.cdl +++ b/parameter_files/fates_params_hydro_default.cdl @@ -6,6 +6,7 @@ dimensions: fates_history_size_bins = 13 ; fates_hydr_organs = 4 ; fates_prt_organs = 6 ; + fates_leafage_class = 2; fates_litterclass = 6 ; fates_pft = 2 ; fates_scalar = 1 ; @@ -284,7 +285,7 @@ variables: float fates_leaf_jmaxse(fates_pft) ; fates_leaf_jmaxse:units = "J/mol/K" ; fates_leaf_jmaxse:long_name = "entropy term for jmax" ; - float fates_leaf_long(fates_pft) ; + float fates_leaf_long(fates_leafage_class, fates_pft) ; fates_leaf_long:units = "yr" ; fates_leaf_long:long_name = "Leaf longevity (ie turnover timescale)" ; float fates_leaf_slamax(fates_pft) ; @@ -305,7 +306,7 @@ variables: float fates_leaf_tpuse(fates_pft) ; fates_leaf_tpuse:units = "J/mol/K" ; fates_leaf_tpuse:long_name = "entropy term for tpu" ; - float fates_leaf_vcmax25top(fates_pft) ; + float fates_leaf_vcmax25top(fates_leafage_class, fates_pft) ; fates_leaf_vcmax25top:units = "umol CO2/m^2/s" ; fates_leaf_vcmax25top:long_name = "maximum carboxylation rate of Rub. at 25C, canopy top" ; float fates_leaf_vcmaxha(fates_pft) ; @@ -873,7 +874,8 @@ data: fates_leaf_jmaxse = 495, 495 ; - fates_leaf_long = 1.5, 1.5 ; + fates_leaf_long = 0.75, 0.75, + 0.75, 0.75 ; fates_leaf_slamax = 0.0954, 0.0954 ; @@ -887,7 +889,8 @@ data: fates_leaf_tpuse = 490, 490 ; - fates_leaf_vcmax25top = 50, 50 ; + fates_leaf_vcmax25top = 50, 50, + 50, 50; fates_leaf_vcmaxha = 65330, 65330 ; From ad4712a99802ff1c128368100327cec61c50650d Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 31 Jan 2019 12:21:07 -0800 Subject: [PATCH 30/34] Added parameters to hydro defaults that got lost between commmits --- parameter_files/fates_params_hydro_default.cdl | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/parameter_files/fates_params_hydro_default.cdl b/parameter_files/fates_params_hydro_default.cdl index 59106459e0..d421573e4c 100644 --- a/parameter_files/fates_params_hydro_default.cdl +++ b/parameter_files/fates_params_hydro_default.cdl @@ -172,6 +172,9 @@ variables: float fates_allom_stmode(fates_pft) ; fates_allom_stmode:units = "index" ; fates_allom_stmode:long_name = "storage allometry function index" ; + float fates_allom_frbstor_repro(fates_pft) ; + fates_allom_frbstor_repro:units = "fraction" ; + fates_allom_frbstor_repro:long_name = "fraction of bstore goes to reproduction after plant dies" ; float fates_branch_turnover(fates_pft) ; fates_branch_turnover:units = "yr-1" ; fates_branch_turnover:long_name = "turnover time of branches" ; @@ -588,6 +591,10 @@ variables: float fates_part_dens ; fates_part_dens:units = "kg/m2" ; fates_part_dens:long_name = "spitfire parameter, oven dry particle density, Table A1 Thonicke et al 2010" ; + float fates_soil_salinity(fates_scalar) ; + fates_soil_salinity:units = "ppt" ; + fates_soil_salinity:long_name = "soil salinity used for model when not coupled to dynamic soil salinity" ; + // global attributes: :history = "This file was made from FatesPFTIndexSwapper.py \n", @@ -756,6 +763,8 @@ data: fates_allom_stmode = 1, 1 ; + fates_allom_frbstor_repro = 0.0, 0.0; + fates_branch_turnover = 50, 50 ; fates_prt_nitr_stoich_p1 = @@ -1077,4 +1086,6 @@ data: fates_miner_total = 0.055 ; fates_part_dens = 513 ; + + fates_soil_salinity = 0.4 ; } From 2d4a2fc97a78ddf1e7e08597fd5461b7d6cf9319 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 31 Jan 2019 15:46:29 -0700 Subject: [PATCH 31/34] Reduced maximum number of leaf layers to more reasonable number (LAI=20 in crown), and increased max cohorts per patch to 100 for safety accomodating large numbers of PFTs --- main/EDTypesMod.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 7d8266051c..328f7d73e5 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -16,13 +16,13 @@ module EDTypesMod save integer, parameter :: maxPatchesPerSite = 10 ! maximum number of patches to live on a site - integer, parameter :: maxCohortsPerPatch = 80 ! maximum number of cohorts per patch + integer, parameter :: maxCohortsPerPatch = 100 ! maximum number of cohorts per patch integer, parameter :: nclmax = 3 ! Maximum number of canopy layers integer, parameter :: ican_upper = 1 ! Nominal index for the upper canopy integer, parameter :: ican_ustory = 2 ! Nominal index for understory in two-canopy system - integer, parameter :: nlevleaf = 40 ! number of leaf layers in canopy layer + integer, parameter :: nlevleaf = 20 ! number of leaf layers in canopy layer integer, parameter :: maxpft = 15 ! maximum number of PFTs allowed ! the parameter file may determine that fewer ! are used, but this helps allocate scratch From 699b1812f07c433725f6eb266c5d3c5439cdb093 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 31 Jan 2019 21:48:16 -0700 Subject: [PATCH 32/34] Minor fix on canopy layer type (integer) in restart file --- biogeochem/EDCohortDynamicsMod.F90 | 2 +- main/FatesRestartInterfaceMod.F90 | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 50e7f7fe8b..b8cddf4a3b 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -603,7 +603,7 @@ subroutine terminate_cohorts( currentSite, currentPatch, level ) endif endif - ! In the third canopy layer + ! Outside the maximum canopy layer if (currentCohort%canopy_layer > nclmax ) then terminate = 1 if ( debug ) then diff --git a/main/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index fb9b316f45..8beb537ba4 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -677,8 +677,8 @@ subroutine define_restart_vars(this, initialize_variables) ! 1D cohort Variables ! ----------------------------------------------------------------------------------- - call this%set_restart_var(vname='fates_canopy_layer', vtype=cohort_r8, & - long_name='ed cohort - canopy_layer', units='unitless', flushval = flushzero, & + call this%set_restart_var(vname='fates_canopy_layer', vtype=cohort_int, & + long_name='ed cohort - canopy_layer', units='unitless', flushval = flushinvalid, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_canopy_layer_co ) call this%set_restart_var(vname='fates_canopy_layer_yesterday', vtype=cohort_r8, & @@ -1448,7 +1448,7 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_ncohort_pa => this%rvars(ir_ncohort_pa)%int1d, & rio_solar_zenith_flag_pa => this%rvars(ir_solar_zenith_flag_pa)%int1d, & rio_solar_zenith_angle_pa => this%rvars(ir_solar_zenith_angle_pa)%r81d, & - rio_canopy_layer_co => this%rvars(ir_canopy_layer_co)%r81d, & + rio_canopy_layer_co => this%rvars(ir_canopy_layer_co)%int1d, & rio_canopy_layer_yesterday_co => this%rvars(ir_canopy_layer_yesterday_co)%r81d, & rio_canopy_trim_co => this%rvars(ir_canopy_trim_co)%r81d, & rio_size_class_lasttimestep => this%rvars(ir_size_class_lasttimestep_co)%int1d, & @@ -2144,7 +2144,7 @@ subroutine get_restart_vectors(this, nc, nsites, sites) rio_ncohort_pa => this%rvars(ir_ncohort_pa)%int1d, & rio_solar_zenith_flag_pa => this%rvars(ir_solar_zenith_flag_pa)%int1d, & rio_solar_zenith_angle_pa => this%rvars(ir_solar_zenith_angle_pa)%r81d, & - rio_canopy_layer_co => this%rvars(ir_canopy_layer_co)%r81d, & + rio_canopy_layer_co => this%rvars(ir_canopy_layer_co)%int1d, & rio_canopy_layer_yesterday_co => this%rvars(ir_canopy_layer_yesterday_co)%r81d, & rio_canopy_trim_co => this%rvars(ir_canopy_trim_co)%r81d, & rio_size_class_lasttimestep => this%rvars(ir_size_class_lasttimestep_co)%int1d, & From 1191c88f1c0a7c8ad05ea8ab5316d8072c2b1c31 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Mon, 4 Feb 2019 16:33:20 -0700 Subject: [PATCH 33/34] Moved call to update cohort level biophysical parameters to after the call to set PFT index --- main/FatesRestartInterfaceMod.F90 | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/main/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index 06dbb5ea1f..1c1aafec5b 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -2294,10 +2294,6 @@ subroutine get_restart_vectors(this, nc, nsites, sites) end do end do - - - call UpdateCohortBioPhysRates(ccohort) - !ccohort%vcmax25top !ccohort%jmax25top !ccohort%tpu25top @@ -2335,6 +2331,9 @@ subroutine get_restart_vectors(this, nc, nsites, sites) ccohort%status_coh = rio_status_co(io_idx_co) ccohort%isnew = ( rio_isnew_co(io_idx_co) .eq. new_cohort ) + call UpdateCohortBioPhysRates(ccohort) + + ! Initialize Plant Hydraulics if(hlm_use_planthydro==itrue)then From 1ca837fb4492e211a86284a1b4e16b9471028381 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 5 Feb 2019 09:25:35 -0800 Subject: [PATCH 34/34] Update README.md --- README.md | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/README.md b/README.md index 37d67c497d..423fc32321 100644 --- a/README.md +++ b/README.md @@ -28,3 +28,14 @@ https://github.com/E3SM-Project/E3SM https://github.com/ESCOMP/cesm https://github.com/ESCOMP/ctsm + +## Important Note About Host-Models and Compatible Branches: +------------------------------------------------------------ + +The FATES and E3SM teams maintain compatability of the NGEET/FATES master branch with the **E3SM master** branch. When changes to the FATES API force compatability updates with E3SM, there may be some modest lag time. + +The FATES team maintains compatability of the NGEET/FATES master branch with the **CTSM fates_next_api** branch. Since the FATES team uses this branch for its internal testing, this compatability is tightly (immediately) maintained and these two should always be in sync. However, CTSM master may become out of sync with FATES master for large periods (months) of time. + + + +