merged main and resolved conflicts

biogeochem/EDCanopyStructureMod.F90

@@ -10,6 +10,7 @@ module EDCanopyStructureMod
   use FatesConstantsMod     , only : tinyr8
   use FatesConstantsMod     , only : nearzero
   use FatesConstantsMod     , only : rsnbl_math_prec
+  use FatesConstantsMod     , only : nocomp_bareground
   use FatesGlobals          , only : fates_log
   use EDPftvarcon           , only : EDPftvarcon_inst
   use PRTParametersMod      , only : prt_params
@@ -1369,7 +1370,7 @@ subroutine canopy_summarization( nsites, sites, bc_in )
              endif
 
              ! adding checks for SP and NOCOMP modes.
-             if(currentPatch%nocomp_pft_label.eq.0)then
+             if(currentPatch%nocomp_pft_label.eq.nocomp_bareground)then
                 write(fates_log(),*) 'cohorts in barepatch',currentPatch%total_canopy_area,currentPatch%nocomp_pft_label
                 call endrun(msg=errMsg(sourcefile, __LINE__))
              end if
@@ -1828,7 +1829,7 @@ subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out)
        c = fcolumn(s)
        do while(associated(currentPatch))
 
-          if(currentPatch%nocomp_pft_label.ne.0)then  ! ignore the bare-ground-PFT patch entirely for these BC outs
+          if(currentPatch%nocomp_pft_label.ne.nocomp_bareground)then  ! ignore the bare-ground-PFT patch entirely for these BC outs
 
              ifp = ifp+1
 
@@ -1970,7 +1971,7 @@ subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out)
           currentPatch => sites(s)%oldest_patch
           ifp = 0
           do while(associated(currentPatch))
-             if(currentPatch%nocomp_pft_label.ne.0)then ! for vegetated patches only
+             if(currentPatch%nocomp_pft_label.ne.nocomp_bareground)then ! for vegetated patches only
                 ifp = ifp+1
                 bc_out(s)%canopy_fraction_pa(ifp) = bc_out(s)%canopy_fraction_pa(ifp)/total_patch_area
              endif ! veg patch
@@ -1994,7 +1995,7 @@ subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out)
 
        ! Pass FATES Harvested C to bc_out.
        call UpdateHarvestC(sites(s),bc_out(s))
-       
+
     end do
 
     ! This call to RecruitWaterStorage() makes an accounting of
@@ -2008,7 +2009,6 @@ subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out)
        call RecruitWaterStorage(nsites,sites,bc_out)
     end if
 
-
   end subroutine update_hlm_dynamics
 
   ! =====================================================================================

biogeochem/EDMortalityFunctionsMod.F90

@@ -225,7 +225,8 @@ end subroutine mortality_rates
 
  ! ============================================================================
 
- subroutine Mortality_Derivative( currentSite, currentCohort, bc_in, frac_site_primary)
+ subroutine Mortality_Derivative( currentSite, currentCohort, bc_in, frac_site_primary, &
+         harvestable_forest_c, harvest_tag)
 
     !
     ! !DESCRIPTION:
@@ -234,14 +235,21 @@ subroutine Mortality_Derivative( currentSite, currentCohort, bc_in, frac_site_pr
     ! elsewhere).
     !
     ! !USES:
-
     use FatesInterfaceTypesMod, only : hlm_freq_day
     !
     ! !ARGUMENTS    
     type(ed_site_type), intent(inout), target  :: currentSite
     type(ed_cohort_type),intent(inout), target :: currentCohort
     type(bc_in_type), intent(in)               :: bc_in
     real(r8), intent(in)                       :: frac_site_primary
+
+    real(r8), intent(in) :: harvestable_forest_c(:)   ! total carbon available for logging, kgC site-1
+    integer, intent(out) :: harvest_tag(:)    ! tag to record the harvest status
+                                              ! for the calculation of harvest debt in C-based
+                                              ! harvest mode
+                                              ! 0 - successful;
+                                              ! 1 - unsuccessful since not enough carbon
+                                              ! 2 - not applicable
     !
     ! !LOCAL VARIABLES:
     real(r8) :: cmort    ! starvation mortality rate (fraction per year)
@@ -271,10 +279,7 @@ subroutine Mortality_Derivative( currentSite, currentCohort, bc_in, frac_site_pr
                                bc_in%hlm_harvest_units, &
                                currentCohort%patchptr%anthro_disturbance_label, &
                                currentCohort%patchptr%age_since_anthro_disturbance, &
-                               frac_site_primary)
-
-
-
+                               frac_site_primary, harvestable_forest_c, harvest_tag)
 
     if (currentCohort%canopy_layer > 1)then 
        ! Include understory logging mortality rates not associated with disturbance

biogeochem/EDPatchDynamicsMod.F90

@@ -39,6 +39,7 @@ module EDPatchDynamicsMod
   use EDTypesMod           , only : AREA_INV
   use FatesConstantsMod    , only : rsnbl_math_prec
   use FatesConstantsMod    , only : fates_tiny
+  use FatesConstantsMod    , only : nocomp_bareground
   use FatesInterfaceTypesMod    , only : hlm_use_planthydro
   use FatesInterfaceTypesMod    , only : hlm_numSWb
   use FatesInterfaceTypesMod    , only : bc_in_type
@@ -47,6 +48,7 @@ module EDPatchDynamicsMod
   use FatesInterfaceTypesMod    , only : hlm_use_sp
   use FatesInterfaceTypesMod    , only : hlm_use_nocomp
   use FatesInterfaceTypesMod    , only : hlm_use_fixed_biogeog
+  use FatesInterfaceTypesMod    , only : hlm_num_lu_harvest_cats
   use FatesGlobals         , only : endrun => fates_endrun
   use FatesConstantsMod    , only : r8 => fates_r8
   use FatesConstantsMod    , only : itrue, ifalse
@@ -57,6 +59,9 @@ module EDPatchDynamicsMod
   use EDLoggingMortalityMod, only : logging_litter_fluxes 
   use EDLoggingMortalityMod, only : logging_time
   use EDLoggingMortalityMod, only : get_harvest_rate_area
+  use EDLoggingMortalityMod, only : get_harvest_rate_carbon
+  use EDLoggingMortalityMod, only : get_harvestable_carbon
+  use EDLoggingMortalityMod, only : get_harvest_debt
   use EDParamsMod          , only : fates_mortality_disturbance_fraction
   use FatesAllometryMod    , only : carea_allom
   use FatesAllometryMod    , only : set_root_fraction
@@ -69,6 +74,7 @@ module EDPatchDynamicsMod
   use FatesConstantsMod    , only : n_anthro_disturbance_categories
   use FatesConstantsMod    , only : fates_unset_r8
   use FatesConstantsMod    , only : fates_unset_int
+  use FatesConstantsMod    , only : hlm_harvest_carbon
   use EDCohortDynamicsMod  , only : InitPRTObject
   use EDCohortDynamicsMod  , only : InitPRTBoundaryConditions
   use ChecksBalancesMod,      only : SiteMassStock
@@ -186,9 +192,12 @@ subroutine disturbance_rates( site_in, bc_in)
     real(r8) :: dist_rate_ldist_notharvested
     integer  :: threshold_sizeclass
     integer  :: i_dist
+    integer  :: h_index
     real(r8) :: frac_site_primary
     real(r8) :: harvest_rate
     real(r8) :: tempsum
+    real(r8) :: harvestable_forest_c(hlm_num_lu_harvest_cats)
+    integer  :: harvest_tag(hlm_num_lu_harvest_cats)
 
     !----------------------------------------------------------------------------------------------
     ! Calculate Mortality Rates (these were previously calculated during growth derivatives)
@@ -197,6 +206,9 @@ subroutine disturbance_rates( site_in, bc_in)
 
     ! first calculate the fractino of the site that is primary land
     call get_frac_site_primary(site_in, frac_site_primary)
+
+    ! get available biomass for harvest for all patches
+    call get_harvestable_carbon(site_in, bc_in%site_area, bc_in%hlm_harvest_catnames, harvestable_forest_c)
 
     currentPatch => site_in%oldest_patch
     do while (associated(currentPatch))   
@@ -227,7 +239,9 @@ subroutine disturbance_rates( site_in, bc_in)
                 bc_in%hlm_harvest_units, &
                 currentPatch%anthro_disturbance_label, &
                 currentPatch%age_since_anthro_disturbance, &
-                frac_site_primary)
+                frac_site_primary, &
+                harvestable_forest_c, &
+                harvest_tag)
 
           currentCohort%lmort_direct     = lmort_direct
           currentCohort%lmort_collateral = lmort_collateral
@@ -236,9 +250,12 @@ subroutine disturbance_rates( site_in, bc_in)
 
           currentCohort => currentCohort%taller
        end do
+
        currentPatch => currentPatch%younger
     end do
 
+    call get_harvest_debt(site_in, bc_in, harvest_tag)
+
     ! ---------------------------------------------------------------------------------------------
     ! Calculate Disturbance Rates based on the mortality rates just calculated
     ! ---------------------------------------------------------------------------------------------
@@ -288,6 +305,11 @@ subroutine disturbance_rates( site_in, bc_in)
                                currentCohort%lmort_infra +                           &
                                currentCohort%l_degrad ) *                            &
                                currentCohort%c_area/currentPatch%area
+
+             if(currentPatch%disturbance_rates(dtype_ilog)>1.0) then
+                 write(fates_log(),*) 'See luc mortalities:', currentCohort%lmort_direct, &
+                     currentCohort%lmort_collateral, currentCohort%lmort_infra, currentCohort%l_degrad
+             end if
 
              ! Non-harvested part of the logging disturbance rate
              dist_rate_ldist_notharvested = dist_rate_ldist_notharvested + currentCohort%l_degrad * &
@@ -298,13 +320,19 @@ subroutine disturbance_rates( site_in, bc_in)
        enddo !currentCohort
 
        ! for non-closed-canopy areas subject to logging, add an additional increment of area disturbed
-       ! equivalent to the fradction logged to account for transfer of interstitial ground area to new secondary lands
+       ! equivalent to the fraction logged to account for transfer of interstitial ground area to new secondary lands
        if ( logging_time .and. &
             (currentPatch%area - currentPatch%total_canopy_area) .gt. fates_tiny ) then
           ! The canopy is NOT closed. 
 
-          call get_harvest_rate_area (currentPatch%anthro_disturbance_label, bc_in%hlm_harvest_catnames, &
-               bc_in%hlm_harvest_rates, frac_site_primary, currentPatch%age_since_anthro_disturbance, harvest_rate)
+          if(bc_in%hlm_harvest_units == hlm_harvest_carbon) then
+             call get_harvest_rate_carbon (currentPatch%anthro_disturbance_label, bc_in%hlm_harvest_catnames, &
+                   bc_in%hlm_harvest_rates, currentPatch%age_since_anthro_disturbance, harvestable_forest_c, &
+                   harvest_rate, harvest_tag)
+          else
+             call get_harvest_rate_area (currentPatch%anthro_disturbance_label, bc_in%hlm_harvest_catnames, &
+                  bc_in%hlm_harvest_rates, frac_site_primary, currentPatch%age_since_anthro_disturbance, harvest_rate)
+          end if
 
           currentPatch%disturbance_rates(dtype_ilog) = currentPatch%disturbance_rates(dtype_ilog) + &
                (currentPatch%area - currentPatch%total_canopy_area) * harvest_rate / currentPatch%area
@@ -314,6 +342,12 @@ subroutine disturbance_rates( site_in, bc_in)
                (currentPatch%area - currentPatch%total_canopy_area) * harvest_rate / currentPatch%area
        endif
 
+       ! For nocomp mode, we need to prevent producing too small patches, which may produce small patches
+       if ((hlm_use_nocomp .eq. itrue) .and. &
+           (currentPatch%disturbance_rates(dtype_ilog)*currentPatch%area .lt. min_patch_area_forced)) then
+          currentPatch%disturbance_rates(dtype_ilog) = 0._r8
+       end if
+
        ! fraction of the logging disturbance rate that is non-harvested
        if (currentPatch%disturbance_rates(dtype_ilog) .gt. nearzero) then
           currentPatch%fract_ldist_not_harvested = dist_rate_ldist_notharvested / &
@@ -1282,11 +1316,11 @@ subroutine set_patchno( currentSite )
        currentPatch => currentPatch%younger
     enddo
 
-    if(hlm_use_sp.eq.itrue)then
+    if(hlm_use_fixed_biogeog.eq.itrue .and. hlm_use_nocomp.eq.itrue)then
       patchno = 1
       currentPatch => currentSite%oldest_patch
       do while(associated(currentPatch))
-        if(currentPatch%nocomp_pft_label.eq.0)then
+        if(currentPatch%nocomp_pft_label.eq.nocomp_bareground)then
          ! for bareground patch, we make the patch number 0
          ! we also do not count this in the veg. patch numbering scheme.
           currentPatch%patchno = 0
@@ -2027,7 +2061,7 @@ subroutine create_patch(currentSite, new_patch, age, areap, label,nocomp_pft)
     real(r8), intent(in) :: age                  ! notional age of this patch in years
     real(r8), intent(in) :: areap                ! initial area of this patch in m2. 
     integer, intent(in)  :: label                ! anthropogenic disturbance label
-    integer, intent(in)  :: nocomp_pft
+    integer, intent(in)  :: nocomp_pft           ! no competition mode pft label
 
 
     ! Until bc's are pointed to by sites give veg a default temp [K]

biogeochem/EDPhysiologyMod.F90

@@ -23,6 +23,7 @@ module EDPhysiologyMod
   use FatesInterfaceTypesMod, only    : hlm_use_tree_damage
   use FatesConstantsMod, only    : r8 => fates_r8
   use FatesConstantsMod, only    : nearzero
+  use FatesConstantsMod, only    : nocomp_bareground
   use EDPftvarcon      , only    : EDPftvarcon_inst
   use PRTParametersMod , only    : prt_params
   use EDPftvarcon      , only    : GetDecompyFrac
@@ -1752,7 +1753,6 @@ subroutine SeedIn( currentSite, bc_in )
     ! !USES:
     use EDTypesMod, only : area
     use EDTypesMod, only : homogenize_seed_pfts
-
     !
     ! !ARGUMENTS
     type(ed_site_type), intent(inout), target  :: currentSite
@@ -2668,7 +2668,7 @@ subroutine fragmentation_scaler( currentPatch, bc_in)
     catanf(t1) = 11.75_r8 +(29.7_r8 / pi) * atan( pi * 0.031_r8  * ( t1 - 15.4_r8 ))
     catanf_30 = catanf(30._r8)
 
-    if(currentPatch%nocomp_pft_label.ne.0)then
+    if(currentPatch%nocomp_pft_label.ne.nocomp_bareground)then
 
        ! Use the hlm temp and moisture decomp fractions by default
        if ( use_hlm_soil_scalar ) then

biogeochem/FatesAllometryMod.F90

@@ -1085,10 +1085,10 @@ subroutine bstore_allom(d,ipft,crowndamage, canopy_trim,bstore,dbstoredd)
 
      real(r8) :: bl          ! Allometric target leaf biomass
      real(r8) :: dbldd       ! Allometric target change in leaf biomass per cm
+     real(r8) :: blmax       ! Allometric target leaf biomass (UNTRIMMED)
+     real(r8) :: dblmaxdd    ! Allometric target change in leaf biomass per cm (UNTRIMMED)
 
 
-     ! TODO: allom_stmode needs to be added to the parameter file
-
      associate( allom_stmode => prt_params%allom_stmode(ipft), &
                 cushion      => prt_params%cushion(ipft) )
 
@@ -1098,7 +1098,12 @@ subroutine bstore_allom(d,ipft,crowndamage, canopy_trim,bstore,dbstoredd)
 
           call bleaf(d,ipft, crowndamage, canopy_trim, bl, dbldd)
           call bstore_blcushion(d,bl,dbldd,cushion,ipft,bstore,dbstoredd)
-
+
+       case(2) ! Storage is constant proportionality of untrimmed maximum leaf
+          ! biomass (ie cushion * bleaf_max)
+          call blmax_allom(d,ipft,blmax,dblmaxdd)
+          call bstore_blcushion(d,blmax,dblmaxdd,cushion,ipft,bstore,dbstoredd)
+
        case DEFAULT 
           write(fates_log(),*) 'An undefined fine storage allometry was specified: ', &
                 allom_stmode
@@ -2515,7 +2520,6 @@ subroutine ForceDBH( ipft, crowndamage, canopy_trim, d, h, bdead, bl )
      end if
 
      call h_allom(d,ipft,h)
-
      if(counter>20)then
         write(fates_log(),*) 'dbh counter: ',counter,' is woody: ',&
              (prt_params%woody(ipft) == itrue)

biogeophys/EDAccumulateFluxesMod.F90

@@ -13,6 +13,7 @@ module EDAccumulateFluxesMod
   use FatesGlobals, only      : fates_log
   use shr_log_mod , only      : errMsg => shr_log_errMsg
   use FatesConstantsMod , only : r8 => fates_r8
+  use FatesConstantsMod , only : nocomp_bareground
 
 
   implicit none
@@ -64,7 +65,7 @@ subroutine AccumulateFluxes_ED(nsites, sites, bc_in, bc_out, dt_time)
 
        cpatch => sites(s)%oldest_patch
        do while (associated(cpatch))                 
-          if(cpatch%nocomp_pft_label.ne.0)then
+          if(cpatch%nocomp_pft_label.ne.nocomp_bareground)then
              ifp = ifp+1
 
              if( bc_in(s)%filter_photo_pa(ifp) == 3 ) then

biogeophys/EDBtranMod.F90

@@ -8,6 +8,7 @@ module EDBtranMod
   use EDPftvarcon       , only : EDPftvarcon_inst
   use FatesConstantsMod , only : tfrz => t_water_freeze_k_1atm 
   use FatesConstantsMod , only : itrue,ifalse,nearzero
+  use FatesConstantsMod , only : nocomp_bareground
   use EDTypesMod        , only : ed_site_type,       &
        ed_patch_type,      &
        ed_cohort_type,     &
@@ -138,7 +139,7 @@ subroutine btran_ed( nsites, sites, bc_in, bc_out)
        ifp = 0
        cpatch => sites(s)%oldest_patch
        do while (associated(cpatch))                 
-          if(cpatch%nocomp_pft_label.ne.0)then ! only for veg patches
+          if(cpatch%nocomp_pft_label.ne.nocomp_bareground)then ! only for veg patches
              ifp=ifp+1
 
              ! THIS SHOULD REALLY BE A COHORT LOOP ONCE WE HAVE rootfr_ft FOR COHORTS (RGK)

biogeophys/EDSurfaceAlbedoMod.F90

@@ -15,6 +15,7 @@ module EDSurfaceRadiationMod
   use FatesConstantsMod , only : r8 => fates_r8
   use FatesConstantsMod , only : itrue
   use FatesConstantsMod , only : pi_const
+  use FatesConstantsMod , only : nocomp_bareground
   use FatesInterfaceTypesMod , only : bc_in_type
   use FatesInterfaceTypesMod , only : bc_out_type
   use FatesInterfaceTypesMod , only : hlm_numSWb
@@ -100,7 +101,7 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out )
        ifp = 0
        currentpatch => sites(s)%oldest_patch
        do while (associated(currentpatch))
-          if(currentpatch%nocomp_pft_label.ne.0)then
+          if(currentpatch%nocomp_pft_label.ne.nocomp_bareground)then
              ! do not do albedo calculations for bare ground patch in SP mode
              ! and (more impotantly) do not iterate ifp or it will mess up the indexing wherein
              ! ifp=1 is the first vegetated patch.
@@ -1147,7 +1148,7 @@ subroutine ED_SunShadeFracs(nsites, sites,bc_in,bc_out)
      cpatch => sites(s)%oldest_patch
 
      do while (associated(cpatch))
-        if(cpatch%nocomp_pft_label.ne.0)then !only for veg patches
+        if(cpatch%nocomp_pft_label.ne.nocomp_bareground)then !only for veg patches
            ! do not do albedo calculations for bare ground patch in SP mode
            ! and (more impotantly) do not iterate ifp or it will mess up the indexing wherein
            ! ifp=1 is the first vegetated patch.

biogeophys/FatesPlantHydraulicsMod.F90

@@ -42,6 +42,7 @@ module FatesPlantHydraulicsMod
   use FatesConstantsMod, only : cm3_per_m3
   use FatesConstantsMod, only : kg_per_g
   use FatesConstantsMod, only : fates_unset_r8
+  use FatesConstantsMod, only : nocomp_bareground
 
   use EDParamsMod       , only : hydr_kmax_rsurf1
   use EDParamsMod       , only : hydr_kmax_rsurf2
@@ -2546,7 +2547,7 @@ subroutine hydraulics_bc ( nsites, sites, bc_in, bc_out, dtime)
      ifp = 0
      cpatch => sites(s)%oldest_patch
      do while (associated(cpatch))
-        if(cpatch%nocomp_pft_label.ne.0)then
+        if(cpatch%nocomp_pft_label.ne.nocomp_bareground)then
            ifp = ifp + 1
 
            ! ----------------------------------------------------------------------------