handle case where LAI+SAI exceeds array by clipping instead of crashing

biogeochem/EDCanopyStructureMod.F90

@@ -18,8 +18,7 @@ module EDCanopyStructureMod
   use FatesAllometryMod     , only : carea_allom
   use EDCohortDynamicsMod   , only : terminate_cohorts, terminate_cohort, fuse_cohorts
   use EDCohortDynamicsMod   , only : InitPRTObject
-  use FatesAllometryMod     , only : tree_lai
-  use FatesAllometryMod     , only : tree_sai
+  use FatesAllometryMod     , only : tree_lai_sai
   use EDTypesMod            , only : ed_site_type
   use FatesAllometryMod     , only : VegAreaLayer
   use FatesAllometryMod     , only : CrownDepth
@@ -2245,7 +2244,7 @@ subroutine UpdateCohortLAI(currentCohort, canopy_layer_tlai, total_canopy_area)
    ! Update LAI and related variables for a given cohort
 
    ! Uses
-   use EDParamsMod, only : dlower_vai
+   use EDParamsMod, only : dlower_vai, dinc_vai
 
    ! Arguments
    type(fates_cohort_type),intent(inout), target   :: currentCohort
@@ -2254,24 +2253,21 @@ subroutine UpdateCohortLAI(currentCohort, canopy_layer_tlai, total_canopy_area)
 
    ! Local variables
    real(r8) :: leaf_c                              ! leaf carbon [kg]
-
+   real(r8) :: treesai                             ! stem area index within crown m2/m2
+
    ! Obtain the leaf carbon
    leaf_c = currentCohort%prt%GetState(leaf_organ,carbon12_element)
 
-
-   ! Note that tree_lai has an internal check on the canopy locatoin
-   currentCohort%treelai = tree_lai(leaf_c, currentCohort%pft, currentCohort%c_area, &
-        currentCohort%n, currentCohort%canopy_layer,               &
-        canopy_layer_tlai,currentCohort%vcmax25top )
-
+   ! Note that tree_lai has an internal check on the canopy location
+   call  tree_lai_sai(leaf_c, currentCohort%pft, currentCohort%c_area, currentCohort%n,           &
+          currentCohort%canopy_layer, canopy_layer_tlai, currentCohort%vcmax25top, currentCohort%dbh, currentCohort%crowndamage,          &
+          currentCohort%canopy_trim, currentCohort%efstem_coh, 4, currentCohort%treelai, treesai )
+
+   ! Do not update stem area index of SP vegetation
    if (hlm_use_sp .eq. ifalse) then
-      currentCohort%treesai = tree_sai(currentCohort%pft, currentCohort%dbh, currentCohort%crowndamage, &
-                                       currentCohort%canopy_trim, currentCohort%efstem_coh, &
-                                       currentCohort%c_area, currentCohort%n, currentCohort%canopy_layer, &
-                                       canopy_layer_tlai, currentCohort%treelai , &
-                                       currentCohort%vcmax25top,4)
+      currentCohort%treesai = treesai
    end if
-
+   
    ! Number of actual vegetation layers in this cohort's crown
    currentCohort%nv =  count((currentCohort%treelai+currentCohort%treesai) .gt. dlower_vai(:)) + 1
 

biogeochem/EDCohortDynamicsMod.F90

@@ -78,7 +78,6 @@ Module EDCohortDynamicsMod
   use FatesAllometryMod  , only : carea_allom
   use FatesAllometryMod  , only : bstore_allom
   use FatesAllometryMod  , only : ForceDBH
-  use FatesAllometryMod  , only : tree_lai, tree_sai
   use FatesAllometryMod    , only : set_root_fraction
   use PRTGenericMod,          only : prt_carbon_allom_hyp
   use PRTGenericMod,          only : prt_cnp_flex_allom_hyp

biogeochem/EDPhysiologyMod.F90

@@ -45,8 +45,7 @@ module EDPhysiologyMod
   use FatesInterfaceTypesMod, only    : bc_out_type
   use EDCohortDynamicsMod , only : create_cohort, sort_cohorts
   use EDCohortDynamicsMod , only : InitPRTObject
-  use FatesAllometryMod   , only : tree_lai
-  use FatesAllometryMod   , only : tree_sai
+  use FatesAllometryMod   , only : tree_lai_sai
   use FatesAllometryMod   , only : leafc_from_treelai
   use FatesAllometryMod   , only : decay_coeff_vcmax
   use FatesLitterMod      , only : litter_type
@@ -707,18 +706,10 @@ subroutine trim_canopy( currentSite )
 
           leaf_c   = currentCohort%prt%GetState(leaf_organ, carbon12_element)
 
-          currentCohort%treelai = tree_lai(leaf_c, currentCohort%pft, currentCohort%c_area, &
-               currentCohort%n, currentCohort%canopy_layer,               &
-               currentPatch%canopy_layer_tlai,currentCohort%vcmax25top )
-
-          ! We don't need to check on sp mode here since we don't trim_canopy with sp mode
-          currentCohort%treesai = tree_sai(currentCohort%pft, &
-               currentCohort%dbh, currentCohort%crowndamage,  &
-               currentCohort%canopy_trim, &
-               currentCohort%efstem_coh, &
-               currentCohort%c_area, currentCohort%n,currentCohort%canopy_layer,& 
-               currentPatch%canopy_layer_tlai, currentCohort%treelai, &
-               currentCohort%vcmax25top,0 )  
+          call  tree_lai_sai(leaf_c, currentCohort%pft, currentCohort%c_area, currentCohort%n,           &
+               currentCohort%canopy_layer, currentPatch%canopy_layer_tlai, currentCohort%vcmax25top,   &
+               currentCohort%dbh, currentCohort%crowndamage, currentCohort%canopy_trim, &
+               currentCohort%efstem_coh, 0, currentCohort%treelai, currentCohort%treesai )
 
           currentCohort%nv      = count((currentCohort%treelai+currentCohort%treesai) .gt. dlower_vai(:)) + 1
 
@@ -1931,6 +1922,7 @@ subroutine calculate_SP_properties(htop, tlai, tsai, parea, pft, crown_damage,
 
     ! LOCAL VARIABLES:
     real(r8) :: check_treelai       ! check tree LAI against input tlai [m2/m2]
+    real(r8) :: dummy_treesai       ! dummy
     real(r8) :: canopylai(1:nclmax) ! canopy LAI [m2/m2]
     real(r8) :: oldcarea            ! save value of crown area [m2]
 
@@ -1948,17 +1940,17 @@ subroutine calculate_SP_properties(htop, tlai, tsai, parea, pft, crown_damage,
 
     ! calculate leaf carbon from target treelai
     canopylai(:) = 0._r8
-    leaf_c = leafc_from_treelai(tlai, pft, c_area, cohort_n, canopy_layer, vcmax25top)
-
+    leaf_c = leafc_from_treelai(tlai, tsai, pft, c_area, cohort_n, canopy_layer, vcmax25top)
+    
     ! check that the inverse calculation of leafc from treelai is the same as the
     ! standard calculation of treelai from leafc. Maybe can delete eventually?
-    check_treelai = tree_lai(leaf_c, pft, c_area, cohort_n, canopy_layer,                &
-         canopylai, vcmax25top)
 
+    call tree_lai_sai(leaf_c, pft, c_area, cohort_n, canopy_layer, canopylai, vcmax25top, &
+                          dbh, crown_damage, 1.0_r8, 1.0_r8, 11, check_treelai, dummy_treesai)
+
     if (abs(tlai - check_treelai) > area_error_2) then !this is not as precise as nearzero
       write(fates_log(),*) 'error in validate treelai', tlai, check_treelai, tlai - check_treelai
-      write(fates_log(),*) 'tree_lai inputs: ', pft, c_area, cohort_n,                   &
-        canopy_layer, vcmax25top
+      write(fates_log(),*) 'tree_lai inputs: ', pft, c_area, cohort_n, canopy_layer, vcmax25top
        call endrun(msg=errMsg(sourcefile, __LINE__))
     end if
 

biogeochem/FatesAllometryMod.F90

@@ -110,8 +110,6 @@ module FatesAllometryMod
   public :: blmax_allom   ! Generic maximum leaf biomass wrapper
   public :: bleaf         ! Generic actual leaf biomass wrapper
   public :: storage_fraction_of_target ! storage as fraction of leaf biomass
-  public :: tree_lai      ! Calculate tree-level LAI from actual leaf biomass
-  public :: tree_sai      ! Calculate tree-level SAI from tree-level LAI
   public :: bsap_allom    ! Generic sapwood wrapper
   public :: bbgw_allom    ! Generic coarse root wrapper
   public :: bfineroot     ! Generic actual fine root biomass wrapper
@@ -127,7 +125,10 @@ module FatesAllometryMod
                                ! root profiles
   public :: leafc_from_treelai ! Calculate target leaf carbon for a given treelai for SP mode
   public :: VegAreaLayer
-
+
+  public :: tree_lai_sai       ! LAI and SAI calculations must work together, thus they
+                               ! should never be called separately
+
   logical         , parameter :: verbose_logging = .false.
   character(len=*), parameter :: sourcefile = __FILE__
 
@@ -828,39 +829,58 @@ real(r8) function tree_sai(pft, dbh, crowndamage, canopy_trim, elongf_stem, c_ar
 
     tree_sai   =  elongf_stem * prt_params%allom_sai_scaler(pft) * target_lai
 
-    if( (treelai + tree_sai) > (sum(dinc_vai)) )then
-
-       call h_allom(dbh,pft,h)
-
-       write(fates_log(),*) 'The leaf and stem are predicted for a cohort, maxed out the array size'
-       write(fates_log(),*) 'lai: ',treelai
-       write(fates_log(),*) 'sai: ',tree_sai
-       write(fates_log(),*) 'lai+sai: ',treelai+tree_sai
-       write(fates_log(),*) 'target_bleaf: ', target_bleaf
-       write(fates_log(),*) 'area: ', c_area
-       write(fates_log(),*) 'target_lai: ',target_lai
-       write(fates_log(),*) 'dinc_vai:',dinc_vai
-       write(fates_log(),*) 'nlevleaf,sum(dinc_vai):',nlevleaf,sum(dinc_vai)
-       write(fates_log(),*) 'pft: ',pft
-       write(fates_log(),*) 'call id: ',call_id
-       write(fates_log(),*) 'n: ',nplant
-       write(fates_log(),*) 'dbh: ',dbh,' dbh_max: ',prt_params%allom_dbh_maxheight(pft)
-       write(fates_log(),*) 'h: ',h
-       write(fates_log(),*) 'canopy_trim: ',canopy_trim
-       write(fates_log(),*) 'canopy layer: ',cl
-       write(fates_log(),*) 'canopy_tlai: ',canopy_lai(:)
-       write(fates_log(),*) 'vcmax25top: ',vcmax25top
-       call endrun(msg=errMsg(sourcefile, __LINE__))
-    end if
-
+    return
+  end function tree_sai
 
+  ! ============================================================================
+
+  subroutine tree_lai_sai(leaf_c, pft, c_area, nplant, cl, canopy_lai, vcmax25top, &
+                          dbh, crowndamage, canopy_trim, elongf_stem, call_id, & 
+                          treelai, treesai)
 
+    ! This is the public wrapper for calling plant lai and sai
+    ! The purpose of this wrapper is to ensure that they are called together,
+    ! and in the correct order, and that the capping is applied after
+    ! the base calculations are performed.
+
+    real(r8), intent(in) :: leaf_c                    ! plant leaf carbon [kg]
+    integer, intent(in)  :: pft                       ! Plant Functional Type index
+    real(r8), intent(in) :: c_area                    ! areal extent of canopy (m2)
+    real(r8), intent(in) :: nplant                    ! number of individuals in cohort per ha
+    integer, intent(in)  :: cl                        ! canopy layer index
+    real(r8), intent(in) :: canopy_lai(:)             ! total leaf area index of 
+                                                      ! each canopy layer
+    real(r8), intent(in) :: vcmax25top                ! maximum carboxylation rate at canopy
+    real(r8), intent(in) :: dbh
+    integer, intent(in)  :: crowndamage
+    real(r8), intent(in) :: canopy_trim        ! trimming function (0-1)
+    real(r8), intent(in) :: elongf_stem        ! Elongation factor for stems.
+    integer,intent(in)   :: call_id            ! flag specifying where this is called from
+    ! from
+
+    real(r8), intent(out) :: treelai         ! plant LAI [m2 leaf area/m2 crown area]
+    real(r8), intent(out) :: treesai         ! plant SAI [m2 stem area/m2 crown area]
+
+
+    treelai = tree_lai( leaf_c, pft, c_area, nplant, cl, canopy_lai, vcmax25top)
+
+    treesai = tree_sai( pft, dbh, crowndamage, canopy_trim, elongf_stem, c_area, nplant, &
+                             cl, canopy_lai, treelai, vcmax25top, call_id )
+
+    ! Don't allow lai+sai to exceed the vertical discretization bounds
+    if( (treelai + treesai) > (sum(dinc_vai)) )then
+       treelai = sum(dinc_vai) * (1._r8 - prt_params%allom_sai_scaler(pft)) - nearzero
+       treesai = sum(dinc_vai) * prt_params%allom_sai_scaler(pft) - nearzero
+    end if
+
+
     return
-  end function tree_sai
+  end subroutine tree_lai_sai
+
 
 ! =====================================================================================
 
-  real(r8) function leafc_from_treelai( treelai, pft, c_area, nplant, cl, vcmax25top)
+  real(r8) function leafc_from_treelai( treelai, treesai, pft, c_area, nplant, cl, vcmax25top)
 
     ! -----------------------------------------------------------------------------------
     ! Calculates the amount of leaf carbon which is needed to generate a given treelai. 
@@ -869,6 +889,7 @@ real(r8) function leafc_from_treelai( treelai, pft, c_area, nplant, cl, vcmax25t
 
     ! !ARGUMENTS
     real(r8), intent(in) :: treelai                    ! desired tree lai m2/m2
+    real(r8), intent(in) :: treesai                    ! desired tree lai m2/m2
     integer, intent(in)  :: pft                       ! Plant Functional Type index
     real(r8), intent(in) :: c_area                    ! areal extent of canopy (m2)
     real(r8), intent(in) :: nplant                    ! number of individuals in cohort per ha
@@ -909,6 +930,11 @@ real(r8) function leafc_from_treelai( treelai, pft, c_area, nplant, cl, vcmax25t
       call endrun(msg=errMsg(sourcefile, __LINE__))
     endif
 
+    if( (treelai + treesai) > (sum(dinc_vai)) )then
+       write(fates_log(),*) 'SP tree lai cannot exceed sum of dinc_vai',treelai,treesai,sum(dinc_vai)
+       call endrun(msg=errMsg(sourcefile, __LINE__))
+    end if
+
     ! convert PFT-level canopy top and maximum SLA values and convert from m2/gC to m2/kgC
     slat = g_per_kg * prt_params%slatop(pft)
     sla_max = g_per_kg * prt_params%slamax(pft)

biogeochem/FatesCohortMod.F90

@@ -25,7 +25,7 @@ module FatesCohortMod
   use EDPftvarcon,                only : EDPftvarcon_inst
   use FatesSizeAgeTypeIndicesMod, only : sizetype_class_index
   use FatesSizeAgeTypeIndicesMod, only : coagetype_class_index
-  use FatesAllometryMod,          only : carea_allom, tree_lai, tree_sai
+  use FatesAllometryMod,          only : carea_allom, tree_lai_sai
   use PRTAllometricCarbonMod,     only : ac_bc_inout_id_dbh, ac_bc_inout_id_netdc
   use PRTAllometricCarbonMod,     only : ac_bc_in_id_cdamage, ac_bc_in_id_pft
   use PRTAllometricCarbonMod,     only : ac_bc_in_id_ctrim, ac_bc_in_id_lstat
@@ -570,7 +570,8 @@ subroutine Create(this, prt, pft, nn, height, coage, dbh, status,            &
       ! LOCAL VARIABLES:
       integer  :: iage        ! loop counter for leaf age classes
       real(r8) :: leaf_c      ! total leaf carbon [kgC]
-
+      real(r8) :: treesai     ! stem area index within crown [m2/m2]
+
       ! initialize cohort
       call this%Init(prt)
 
@@ -640,15 +641,14 @@ subroutine Create(this, prt, pft, nn, height, coage, dbh, status,            &
       ! Query PARTEH for the leaf carbon [kg]
       leaf_c = this%prt%GetState(leaf_organ, carbon12_element)
 
-      ! calculate tree lai
-      this%treelai = tree_lai(leaf_c, this%pft, this%c_area, this%n,           &
-        this%canopy_layer, can_tlai, this%vcmax25top)
+      call tree_lai_sai(leaf_c, this%pft, this%c_area, this%n,           &
+           this%canopy_layer, can_tlai, this%vcmax25top, this%dbh, this%crowndamage,          &
+           this%canopy_trim, this%efstem_coh, 2, this%treelai, treesai)
 
       if (hlm_use_sp .eq. ifalse) then
-        this%treesai = tree_sai(this%pft, this%dbh, this%crowndamage,          &
-          this%canopy_trim, this%efstem_coh, this%c_area, this%n,              &
-          this%canopy_layer, can_tlai, this%treelai,this%vcmax25top, 2)
+         this%treesai = treesai
       end if
+
 
       call this%InitPRTBoundaryConditions()
 

main/EDMainMod.F90

@@ -85,7 +85,7 @@ module EDMainMod
   use FatesPlantHydraulicsMod  , only : InitPlantHydStates
   use FatesPlantHydraulicsMod  , only : UpdateSizeDepRhizHydProps
   use FatesPlantHydraulicsMod  , only : AccumulateMortalityWaterStorage
-  use FatesAllometryMod        , only : h_allom,tree_sai,tree_lai
+  use FatesAllometryMod        , only : h_allom
   use EDLoggingMortalityMod    , only : IsItLoggingTime
   use EDLoggingMortalityMod    , only : get_harvestable_carbon
   use DamageMainMod            , only : IsItDamageTime

main/FatesConstantsMod.F90

@@ -38,7 +38,9 @@ module FatesConstantsMod
   real(fates_r8), parameter, public :: patchfusion_dbhbin_loweredges(N_DBH_BINS) = &
   (/0._fates_r8, 5._fates_r8, 20._fates_r8, 50._fates_r8, 100._fates_r8, 150._fates_r8/) ! array of bin lower edges for comparing patches
 
-
+
+  real(fates_r8), parameter, public :: min_vai_bin_sum = 5.0_fates_r8   ! The sum of vai increments used to discretize the canopy vertically                                                                                                                           ! must be larger than this number. 
+
   integer , parameter, public :: N_DIST_TYPES = 4          ! Disturbance Modes 1) tree-fall, 2) fire, 3) logging, 4) land-use change
   integer , parameter, public :: dtype_ifall  = 1          ! index for naturally occuring tree-fall generated event
   integer , parameter, public :: dtype_ifire  = 2          ! index for fire generated disturbance event

main/FatesInterfaceMod.F90

@@ -839,7 +839,7 @@ subroutine SetFatesGlobalElements2(use_fates)
       !
       ! --------------------------------------------------------------------------------
 
-      use FatesConstantsMod,      only : fates_check_param_set
+      use FatesConstantsMod,      only : fates_check_param_set, min_vai_bin_sum
 
       logical,intent(in) :: use_fates    ! Is fates turned on?
       integer :: i
@@ -931,6 +931,20 @@ subroutine SetFatesGlobalElements2(use_fates)
             dinc_vai(i) = ED_val_vai_top_bin_width * ED_val_vai_width_increase_factor ** (i-1)
          end do
 
+         if (sum(dinc_vai) < min_vai_bin_sum ) then
+            write(fates_log(), *) 'You specified LAI+SAI bins that add up to a number'
+            write(fates_log(), *) 'that is not reasonably large enough to encapsulate in-canopy'
+            write(fates_log(), *) 'total area indices for large mature trees'
+            write(fates_log(), *) 'sum of vai increments sum(dinc_vai) = ',sum(dinc_vai)
+            write(fates_log(), *) 'minimum allowable user set vai, min_vai_bin_sum = ',min_vai_bin_sum
+            write(fates_log(), *) 'Increase either the top bin width fates_vai_top_bin_width'
+            write(fates_log(), *) 'or the increase factor fates_vai_width_increase_factor'
+            write(fates_log(), *) 'as found in the parameter file.'
+            write(fates_log(), *) 'Or, you can decrease the minimum if you believe its too large'
+            write(fates_log(), *) 'but that is not recommended in most cases.'
+            call endrun(msg=errMsg(sourcefile, __LINE__))
+         end if
+
          ! lower edges of VAI bins       
          do i = 1,nlevleaf
             dlower_vai(i) = sum(dinc_vai(1:i))

main/FatesParameterDerivedMod.F90

@@ -93,7 +93,6 @@ subroutine Init(this,numpft)
     ! local variables
     integer  :: ft                 ! pft index
     integer  :: iage               ! leaf age class index
-    integer  :: c                  ! cwd index
 
     associate( vcmax25top => EDPftvarcon_inst%vcmax25top )