changed cohort fusion to conserve crown area

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 : fates_unset_r8
   use FatesInterfaceMod     , only : hlm_days_per_year
   use EDPftvarcon           , only : EDPftvarcon_inst
   use EDTypesMod            , only : ed_site_type, ed_patch_type, ed_cohort_type
@@ -729,6 +730,7 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in)
      real(r8) :: newn
      real(r8) :: diff
      real(r8) :: dynamic_fusion_tolerance
+     real(r8) :: dbh
 
      integer  :: largersc, smallersc, sc_i        ! indices for tracking the growth flux caused by fusion
      real(r8) :: larger_n, smaller_n
@@ -794,6 +796,7 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in)
                              if( currentCohort%isnew.eqv.nextc%isnew ) then
 
                                 newn = currentCohort%n + nextc%n
+
                                 fusion_took_place = 1         
 
                                 if ( fuse_debug .and. currentCohort%isnew ) then
@@ -822,26 +825,28 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in)
                                 currentCohort%laimemory   = (currentCohort%n*currentCohort%laimemory   &
                                       + nextc%n*nextc%laimemory)/newn
 
-                                currentCohort%dbh         = (currentCohort%n*currentCohort%dbh         &
-                                      + nextc%n*nextc%dbh)/newn
-
-                                call h_allom(currentCohort%dbh,currentCohort%pft,currentCohort%hite)
-
                                 currentCohort%canopy_trim = (currentCohort%n*currentCohort%canopy_trim &
                                       + nextc%n*nextc%canopy_trim)/newn
 
-                                ! -----------------------------------------------------------------
-                                ! If fusion pushed structural biomass to be larger than
-                                ! the allometric target value derived by diameter, we
-                                ! then increase diameter and height until the allometric 
-                                ! target matches actual bdead. (if it is the other way around
-                                ! we then just let the carbon pools grow to fill-out allometry)
-                                ! -----------------------------------------------------------------
+                                ! conserve total crown area during the fusion step, and then calculate dbh of the
+                                ! fused cohort as that which conserves both crown area and the dbh to crown area allometry.  
+                                ! dbh will be updated in the next growth step in the (likely) event that dbh to structural 
+                                ! biomass allometry is exceeded. if using a capped crown area allometry and above the cap, 
+                                ! then calculate as the weighted average of fusing cohorts' dbh
+
+                                currentCohort%c_area = currentCohort%c_area + nextc%c_area
+
+                                call carea_allom(dbh,newn,currentSite%spread,currentCohort%pft,&
+                                     currentCohort%c_area,inverse=.true.)
 
-                                if( EDPftvarcon_inst%woody(currentCohort%pft) == itrue ) then
-                                   call StructureResetOfDH( currentCohort%prt%GetState(struct_organ,all_carbon_elements), currentCohort%pft, &
-                                         currentCohort%canopy_trim, currentCohort%dbh, currentCohort%hite )
-                                end if
+                                if (dbh .eq. fates_unset_r8) then
+                                   currentCohort%dbh = (currentCohort%n*currentCohort%dbh         &
+                                        + nextc%n*nextc%dbh)/newn
+                                else
+                                   currentCohort%dbh = dbh
+                                endif
+
+                                call h_allom(currentCohort%dbh,currentCohort%pft,currentCohort%hite)
 
                                 call sizetype_class_index(currentCohort%dbh,currentCohort%pft, &
                                       currentCohort%size_class,currentCohort%size_by_pft_class)

biogeochem/FatesAllometryMod.F90

@@ -90,6 +90,7 @@ module FatesAllometryMod
   use FatesConstantsMod, only : cm2_per_m2
   use FatesConstantsMod, only : kg_per_Megag
   use FatesConstantsMod, only : calloc_abs_error
+  use FatesConstantsMod, only : fates_unset_r8
   use shr_log_mod      , only : errMsg => shr_log_errMsg
   use FatesGlobals     , only : fates_log
   use FatesGlobals     , only : endrun => fates_endrun
@@ -448,15 +449,17 @@ end subroutine blmax_allom
   ! Generic crown area allometry wrapper
   ! ============================================================================
 
-  subroutine carea_allom(d,nplant,site_spread,ipft,c_area)
+  subroutine carea_allom(d,nplant,site_spread,ipft,c_area,inverse)
 
-     real(r8),intent(in)    :: d           ! plant diameter [cm]
+     real(r8),intent(inout) :: d           ! plant diameter [cm]
      real(r8),intent(in)    :: site_spread ! site level spread factor (crowdedness)
      real(r8),intent(in)    :: nplant      ! number of plants [1/ha]
      integer(i4),intent(in) :: ipft        ! PFT index
-     real(r8),intent(out)   :: c_area      ! crown area per cohort (m2)
+     real(r8),intent(inout) :: c_area      ! crown area per cohort (m2)
+     logical,optional,intent(in) :: inverse! if true, calculate dbh from crown area instead of crown area from dbh
 
      real(r8)               :: d_eff     ! Effective diameter (cm)
+     logical                :: do_inverse, capped_allom
 
      associate( dbh_maxh    => EDPftvarcon_inst%allom_dbh_maxheight(ipft), &
                 allom_lmode => EDPftvarcon_inst%allom_lmode(ipft),  &
@@ -465,24 +468,45 @@ subroutine carea_allom(d,nplant,site_spread,ipft,c_area)
                 d2ca_min    => EDPftvarcon_inst%allom_d2ca_coefficient_min(ipft), &
                 d2ca_max    => EDPftvarcon_inst%allom_d2ca_coefficient_max(ipft))
 
+       if( .not. present(inverse) ) then 
+          do_inverse = .false.
+       else
+          do_inverse = inverse
+       endif
+
+       if (do_inverse) then
+          c_area = c_area / nplant
+       endif
+
        select case(int(allom_lmode))
        case(1)
           d_eff = min(d,dbh_maxh)
-          call carea_2pwr(d_eff,site_spread,d2bl_p2,d2bl_ediff,d2ca_min,d2ca_max,c_area)
+          call carea_2pwr(d_eff,site_spread,d2bl_p2,d2bl_ediff,d2ca_min,d2ca_max,c_area,do_inverse)
+          capped_allom = .true.
        case(2)   ! "2par_pwr")
-          call carea_2pwr(d,site_spread,d2bl_p2,d2bl_ediff,d2ca_min,d2ca_max,c_area)
+          call carea_2pwr(d,site_spread,d2bl_p2,d2bl_ediff,d2ca_min,d2ca_max,c_area,do_inverse)
+          capped_allom = .false.
        case(3)
           d_eff = min(d,dbh_maxh)
-          call carea_2pwr(d_eff,site_spread,d2bl_p2,d2bl_ediff,d2ca_min,d2ca_max,c_area)
+          call carea_2pwr(d_eff,site_spread,d2bl_p2,d2bl_ediff,d2ca_min,d2ca_max,c_area,do_inverse)
+          capped_allom = .true.
        case DEFAULT
-         write(fates_log(),*) 'An undefined leaf allometry was specified: ', &
+          write(fates_log(),*) 'An undefined leaf allometry was specified: ', &
                allom_lmode
-         write(fates_log(),*) 'Aborting'
-         call endrun(msg=errMsg(sourcefile, __LINE__))
-      end select
+          write(fates_log(),*) 'Aborting'
+          call endrun(msg=errMsg(sourcefile, __LINE__))
+       end select
+
 
-      c_area = c_area * nplant
+       if (capped_allom .and. do_inverse) then
+          if (d_eff .lt. dbh_maxh) then
+             d = d_eff
+          else
+             d = fates_unset_r8
+          endif
+       endif
 
+       c_area = c_area * nplant
 
     end associate
     return
@@ -1848,20 +1872,21 @@ end subroutine h2d_martcano
   ! =============================================================================
 
 
-  subroutine carea_2pwr(d,spread,d2bl_p2,d2bl_ediff,d2ca_min,d2ca_max,c_area)
+  subroutine carea_2pwr(d,spread,d2bl_p2,d2bl_ediff,d2ca_min,d2ca_max,c_area,inverse)
 
      ! ============================================================================
      ! Calculate area of ground covered by entire cohort. (m2)
      ! Function of DBH (cm) canopy spread (m/cm) and number of individuals. 
      ! ============================================================================
 
-     real(r8),intent(in) :: d           ! diameter [cm]
+     real(r8),intent(inout) :: d           ! diameter [cm]
      real(r8),intent(in) :: spread      ! site level relative spread score [0-1]
      real(r8),intent(in) :: d2bl_p2     ! parameter 2 in the diameter->bleaf allometry (exponent)
      real(r8),intent(in) :: d2bl_ediff  ! area difference factor in the diameter-bleaf allometry (exponent)
      real(r8),intent(in) :: d2ca_min    ! minimum diameter to crown area scaling factor
      real(r8),intent(in) :: d2ca_max    ! maximum diameter to crown area scaling factor
-     real(r8),intent(out) :: c_area     ! crown area for one plant [m2]
+     real(r8),intent(inout) :: c_area     ! crown area for one plant [m2]
+     logical,intent(in)  :: inverse     ! if true, calculate dbh from crown area rather than its reverse
 
      real(r8)            :: crown_area_to_dbh_exponent
      real(r8)            :: spreadterm  ! Effective 2bh to crown area scaling factor
@@ -1885,7 +1910,11 @@ subroutine carea_2pwr(d,spread,d2bl_p2,d2bl_ediff,d2ca_min,d2ca_max,c_area)
 
      spreadterm = spread * d2ca_max + (1._r8 - spread) * d2ca_min
 
-     c_area = spreadterm * d ** crown_area_to_dbh_exponent
+     if ( .not. inverse) then
+        c_area = spreadterm * d ** crown_area_to_dbh_exponent
+     else
+        d = (c_area / spreadterm) ** (1./crown_area_to_dbh_exponent)
+     endif
 
   end subroutine carea_2pwr
 

main/FatesConstantsMod.F90

@@ -17,7 +17,7 @@ module FatesConstantsMod
 
   ! Unset and various other 'special' values
   integer, parameter :: fates_unset_int = -9999
-
+  real(fates_r8), parameter :: fates_unset_r8 = -9999._fates_r8
 
   ! Integer equivalent of true  (in case some compilers dont auto convert)
   integer, parameter :: itrue = 1