Resolved conflict in inventory init between phen refactor and allowing pft=0

Author: rgknox

biogeochem/EDCohortDynamicsMod.F90

@@ -14,6 +14,7 @@ module EDCohortDynamicsMod
   use FatesConstantsMod     , only : itrue,ifalse
   use FatesConstantsMod     , only : fates_unset_r8
   use FatesConstantsMod     , only : nearzero
+  use FatesConstantsMod     , only : calloc_abs_error
   use FatesInterfaceMod     , only : hlm_days_per_year
   use FatesInterfaceMod     , only : nleafage
   use EDPftvarcon           , only : EDPftvarcon_inst
@@ -47,6 +48,10 @@ module EDCohortDynamicsMod
   use FatesSizeAgeTypeIndicesMod, only : sizetype_class_index
   use FatesAllometryMod  , only : bleaf
   use FatesAllometryMod  , only : bfineroot
+  use FatesAllometryMod  , only : bsap_allom
+  use FatesAllometryMod  , only : bagw_allom
+  use FatesAllometryMod  , only : bbgw_allom
+  use FatesAllometryMod  , only : bdead_allom
   use FatesAllometryMod  , only : h_allom
   use FatesAllometryMod  , only : carea_allom
   use FatesAllometryMod  , only : ForceDBH
@@ -93,6 +98,7 @@ module EDCohortDynamicsMod
   public :: count_cohorts
   public :: InitPRTCohort
   public :: UpdateCohortBioPhysRates
+  public :: EvaluateAndCorrectDBH
 
   logical, parameter :: debug  = .false. ! local debug flag
 
@@ -1706,4 +1712,90 @@ end subroutine UpdateCohortBioPhysRates
 
   ! ============================================================================
 
+
+  subroutine EvaluateAndCorrectDBH(currentCohort,delta_dbh,delta_hite)
+
+    ! -----------------------------------------------------------------------------------
+    ! If the current diameter of a plant is somehow less than what is allometrically 
+    ! consistent with stuctural biomass (or, in the case of grasses, leaf biomass) 
+    ! then correct (increase) the dbh to match that.
+    ! -----------------------------------------------------------------------------------
+
+    ! argument
+    type(ed_cohort_type),intent(inout) :: currentCohort
+    real(r8),intent(out)               :: delta_dbh
+    real(r8),intent(out)               :: delta_hite
+    
+    ! locals
+    real(r8) :: dbh
+    real(r8) :: canopy_trim
+    integer  :: ipft
+    real(r8) :: sapw_area
+    real(r8) :: target_sapw_c
+    real(r8) :: target_agw_c
+    real(r8) :: target_bgw_c
+    real(r8) :: target_struct_c
+    real(r8) :: target_leaf_c
+    real(r8) :: struct_c
+    real(r8) :: hite_out
+    real(r8) :: leaf_c
+    
+    dbh  = currentCohort%dbh
+    ipft = currentCohort%pft
+    canopy_trim = currentCohort%canopy_trim
+
+    delta_dbh   = 0._r8
+    delta_hite  = 0._r8
+    
+    if( EDPftvarcon_inst%woody(ipft) == itrue) then
+
+       struct_c = currentCohort%prt%GetState(struct_organ, all_carbon_elements)
+    
+       ! Target sapwood biomass according to allometry and trimming [kgC]
+       call bsap_allom(dbh,ipft,canopy_trim,sapw_area,target_sapw_c)
+       
+       ! Target total above ground biomass in woody/fibrous tissues  [kgC]
+       call bagw_allom(dbh,ipft,target_agw_c)
+       
+       ! Target total below ground biomass in woody/fibrous tissues [kgC] 
+       call bbgw_allom(dbh,ipft,target_bgw_c)
+       
+       ! Target total dead (structrual) biomass [kgC]
+       call bdead_allom( target_agw_c, target_bgw_c, target_sapw_c, ipft, target_struct_c)
+       
+       ! ------------------------------------------------------------------------------------
+       ! If structure is larger than target, then we need to correct some integration errors
+       ! by slightly increasing dbh to match it.
+       ! For grasses, if leaf biomass is larger than target, then we reset dbh to match
+       ! -----------------------------------------------------------------------------------
+       
+       if( (struct_c - target_struct_c ) > calloc_abs_error ) then
+          call ForceDBH( ipft, canopy_trim, dbh, hite_out, bdead=struct_c )
+          delta_dbh = dbh - currentCohort%dbh 
+          delta_hite = hite_out - currentCohort%hite
+          currentCohort%dbh  = dbh
+          currentCohort%hite = hite_out
+       end if
+       
+    else
+
+       ! This returns the sum of leaf carbon over all (age) bins
+       leaf_c  = currentCohort%prt%GetState(leaf_organ, all_carbon_elements)
+
+       ! Target leaf biomass according to allometry and trimming
+       call bleaf(dbh,ipft,canopy_trim,target_leaf_c)
+
+       if( ( leaf_c - target_leaf_c ) > calloc_abs_error ) then
+          call ForceDBH( ipft, canopy_trim, dbh, hite_out, bl=leaf_c )
+          delta_dbh = dbh - currentCohort%dbh 
+          delta_hite = hite_out - currentCohort%hite
+          currentCohort%dbh = dbh
+          currentCohort%hite = hite_out
+       end if
+       
+    end if
+    return
+  end subroutine EvaluateAndCorrectDBH
+  
+
 end module EDCohortDynamicsMod

biogeochem/EDPhysiologyMod.F90

@@ -987,6 +987,10 @@ subroutine seeds_in( currentSite, cp_pnt )
                  EDPftvarcon_inst%seed_rain(p) !KgC/m2/year
        currentSite%seed_rain_flux(p) = currentSite%seed_rain_flux(p) + &
                  EDPftvarcon_inst%seed_rain(p) * currentPatch%area/AREA !KgC/m2/year
+
+       currentSite%flux_in = currentSite%flux_in + &
+             EDPftvarcon_inst%seed_rain(p) * currentPatch%area * hlm_freq_day
+
     enddo
 
 

biogeochem/FatesAllometryMod.F90

@@ -2200,6 +2200,7 @@ end function decay_coeff_kn
 
   ! =====================================================================================
 
+
   subroutine ForceDBH( ipft, canopy_trim, d, h, bdead, bl )
 
      ! =========================================================================

main/EDMainMod.F90

@@ -24,6 +24,7 @@ module EDMainMod
   use EDCohortDynamicsMod      , only : fuse_cohorts
   use EDCohortDynamicsMod      , only : sort_cohorts
   use EDCohortDynamicsMod      , only : count_cohorts
+  use EDCohortDynamicsMod      , only : EvaluateAndCorrectDBH
   use EDPatchDynamicsMod       , only : disturbance_rates
   use EDPatchDynamicsMod       , only : fuse_patches
   use EDPatchDynamicsMod       , only : spawn_patches
@@ -273,7 +274,9 @@ subroutine ed_integrate_state_variables(currentSite, bc_in )
     real(r8) :: hite_old              ! height of plant before daily PRT [m]
     logical  :: is_drought            ! logical for if the plant (site) is in a drought state
     real(r8) :: leaf_c
-    
+    real(r8) :: delta_dbh             ! correction for dbh
+    real(r8) :: delta_hite            ! correction for hite
+
     !-----------------------------------------------------------------------
 
     small_no = 0.0000000000_r8  ! Obviously, this is arbitrary.  RF - changed to zero
@@ -317,8 +320,7 @@ subroutine ed_integrate_state_variables(currentSite, bc_in )
           ! Apply Plant Allocation and Reactive Transport
           ! -----------------------------------------------------------------------------
 
-          hite_old = currentCohort%hite
-          dbh_old  = currentCohort%dbh
+          
 
           ! -----------------------------------------------------------------------------
           !  Identify the net carbon gain for this dynamics interval
@@ -380,6 +382,15 @@ subroutine ed_integrate_state_variables(currentSite, bc_in )
                currentPatch%canopy_layer_tlai, currentCohort%treelai,currentCohort%vcmax25top,7 )
 
 
+          ! If the current diameter of a plant is somehow less than what is consistent
+          ! with what is allometrically consistent with the stuctural biomass, then
+          ! correct the dbh to match.
+
+          call EvaluateAndCorrectDBH(currentCohort,delta_dbh,delta_hite)
+
+          hite_old = currentCohort%hite
+          dbh_old  = currentCohort%dbh
+
           ! Conduct Growth (parteh)
           call currentCohort%prt%DailyPRT()
           call currentCohort%prt%CheckMassConservation(ft,5)

main/FatesInventoryInitMod.F90

@@ -809,6 +809,7 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, &
       real(r8) :: b_dead
       real(r8) :: b_store 
       real(r8) :: a_sapwood  ! area of sapwood at reference height [m2]
+      integer                                     :: i_pft, ncohorts_to_create
 
 
       character(len=128),parameter    :: wr_fmt = &
@@ -857,9 +858,9 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, &
          call endrun(msg=errMsg(sourcefile, __LINE__))
       end if
 
-      if (c_pft <= 0 ) then
+      if (c_pft < 0 ) then
          write(fates_log(), *) 'inventory pft: ',c_pft
-         write(fates_log(), *) 'The inventory produced a cohort with <=0 pft index'
+         write(fates_log(), *) 'The inventory produced a cohort with <0 pft index'
          call endrun(msg=errMsg(sourcefile, __LINE__))
       end if
 
@@ -886,64 +887,91 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, &
          write(fates_log(), *) 'The inventory produced a cohort with very large density /m2'
          call endrun(msg=errMsg(sourcefile, __LINE__))
       end if
-      
-      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
-      temp_cohort%dbh         = c_dbh
-      call h_allom(c_dbh,c_pft,temp_cohort%hite)
-      temp_cohort%canopy_trim = 1.0_r8
+      if (c_pft .eq. 0 ) then
+         write(fates_log(), *) 'inventory pft: ',c_pft
+         write(fates_log(), *) 'SPECIAL CASE TRIGGERED: PFT == 0 and therefore this subroutine'
+         write(fates_log(), *) 'will assign a cohort with n = n_orig/numpft to every cohort in range 1 to numpft'
+         ncohorts_to_create = numpft
+      else
+         ncohorts_to_create = 1
+      end if
 
-      ! Calculate total above-ground biomass from allometry
+      do i_pft = 1,ncohorts_to_create
+         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)
 
-      call bagw_allom(temp_cohort%dbh,c_pft,b_agw)
-      ! Calculate coarse root biomass from allometry
-      call bbgw_allom(temp_cohort%dbh,c_pft,b_bgw)
-      
-      ! Calculate the leaf biomass (calculates a maximum first, then applies canopy trim
-      ! and sla scaling factors)
-      call bleaf(temp_cohort%dbh,c_pft,temp_cohort%canopy_trim,b_leaf)
-      
-      ! Calculate fine root biomass
-      call bfineroot(temp_cohort%dbh,c_pft,temp_cohort%canopy_trim,b_fineroot)
-      
-      ! Calculate sapwood biomass
-      call bsap_allom(temp_cohort%dbh,c_pft,temp_cohort%canopy_trim, a_sapwood, b_sapwood)
-      
-      call bdead_allom( b_agw, b_bgw, b_sapwood, c_pft, b_dead )
+         if (c_pft .ne. 0 ) then
+            ! normal case: assign each cohort to its specified PFT
+            temp_cohort%pft         = c_pft
+         else
+            ! special case, make an identical cohort for each PFT
+            temp_cohort%pft         = i_pft
+         endif
 
-      call bstore_allom(temp_cohort%dbh, c_pft, temp_cohort%canopy_trim, b_store)
-      
-      temp_cohort%laimemory = 0._r8
-      cstatus = leaves_on
-            
-      if( EDPftvarcon_inst%season_decid(c_pft) == itrue .and. &
-           any(csite%cstatus == [phen_cstat_nevercold,phen_cstat_iscold])) then
-         temp_cohort%laimemory = b_leaf
-         b_leaf  = 0._r8
-         cstatus = leaves_off
-      endif
-      
-      if ( EDPftvarcon_inst%stress_decid(c_pft) == itrue .and. &
-           any(csite%dstatus == [phen_dstat_timeoff,phen_dstat_moistoff])) 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.
-      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, equal_leaf_aclass, bc_in)
+         temp_cohort%n           = c_nplant * cpatch%area / real(ncohorts_to_create,r8)
+         temp_cohort%dbh         = c_dbh
 
-      
-      deallocate(temp_cohort) ! get rid of temporary cohort
+         call h_allom(c_dbh,temp_cohort%pft,temp_cohort%hite)
+         temp_cohort%canopy_trim = 1.0_r8
+
+         ! Calculate total above-ground biomass from allometry
+
+         call bagw_allom(temp_cohort%dbh,temp_cohort%pft,b_agw)
+         ! Calculate coarse root biomass from allometry
+         call bbgw_allom(temp_cohort%dbh,temp_cohort%pft,b_bgw)
+
+         ! Calculate the leaf biomass (calculates a maximum first, then applies canopy trim
+         ! and sla scaling factors)
+         call bleaf(temp_cohort%dbh,temp_cohort%pft,temp_cohort%canopy_trim,b_leaf)
+
+         ! Calculate fine root biomass
+         call bfineroot(temp_cohort%dbh,temp_cohort%pft,temp_cohort%canopy_trim,b_fineroot)
+
+         ! Calculate sapwood biomass
+         call bsap_allom(temp_cohort%dbh,temp_cohort%pft,temp_cohort%canopy_trim, a_sapwood, b_sapwood)
+
+         call bdead_allom( b_agw, b_bgw, b_sapwood, temp_cohort%pft, b_dead )
+
+         call bstore_allom(temp_cohort%dbh, temp_cohort%pft, temp_cohort%canopy_trim, b_store)
+
+         temp_cohort%laimemory = 0._r8
+         cstatus = leaves_on
+
+
+         if( EDPftvarcon_inst%season_decid(temp_cohort%pft) == itrue .and. &
+              any(csite%cstatus == [phen_cstat_nevercold,phen_cstat_iscold])) then
+            temp_cohort%laimemory = b_leaf
+            b_leaf  = 0._r8
+            cstatus = leaves_off
+         endif
+         
+         if ( EDPftvarcon_inst%stress_decid(temp_cohort%pft) == itrue .and. &
+              any(csite%dstatus == [phen_dstat_timeoff,phen_dstat_moistoff])) 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.
+         if( debug_inv) then
+            write(fates_log(),*) 'calling create_cohort: ', 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, cstatus, rstatus, temp_cohort%canopy_trim, &
+                 1, csite%spread
+         endif
+
+         call create_cohort(csite, cpatch, 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, cstatus, rstatus, temp_cohort%canopy_trim, &
+              1, csite%spread, equal_leaf_aclass, bc_in)
+
+         deallocate(temp_cohort) ! get rid of temporary cohort
+      end do
 
       return
-   end subroutine set_inventory_edcohort_type1
+    end subroutine set_inventory_edcohort_type1
 
 end module FatesInventoryInitMod