Minor conflict resolution with master, both modified definition of c3…

…psn and clumping index, (master had re-named its entry in the parameter file.

biogeochem/EDCanopyStructureMod.F90

@@ -148,20 +148,6 @@ subroutine canopy_structure( currentSite , bc_in )
          enddo
 
 
-         ! ------------------------------------------------------------------------------
-         ! Check patch area to prevent numerical weirdness
-         ! ------------------------------------------------------------------------------
-
-         if (currentPatch%area .lt. min_patch_area) then  
-
-            write(fates_log(),*) 'An incredibly small patch exists that should'
-            write(fates_log(),*) 'had been fused or culled already'
-            write(fates_log(),*) 'currentPatch%area: ',currentPatch%area
-            write(fates_log(),*) 'min_patch_area: ',min_patch_area
-            call endrun(msg=errMsg(sourcefile, __LINE__))
-
-         end if
-
          ! Does any layer have excess area in it? Keep going until it does not...
          patch_area_counter = 0
          area_not_balanced = .true.
@@ -1099,6 +1085,13 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si)
        currentPatch%canopy_area_profile(:,:,:)  = 0._r8       
        currentPatch%canopy_mask(:,:)            = 0
 
+       ! ------------------------------------------------------------------------------
+       ! It is remotely possible that in deserts we will not have any canopy
+       ! area, ie not plants at all...
+       ! ------------------------------------------------------------------------------
+
+       if (currentPatch%total_canopy_area > tiny(currentPatch%total_canopy_area)) then
+
        currentCohort => currentPatch%shortest
        do while(associated(currentCohort)) 
 
@@ -1230,12 +1223,6 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si)
           ! and canopy area to the accumulators. 
           ! -----------------------------------------------------------------------------
 
-          ! ------------------------------------------------------------------------------
-          ! It is remotely possible that in deserts we will not have any canopy
-          ! area, ie not plants at all...
-          ! ------------------------------------------------------------------------------
-
-          if (currentPatch%total_canopy_area > tiny(currentPatch%total_canopy_area)) then
 
              currentCohort => currentPatch%shortest
              do while(associated(currentCohort))   
@@ -1393,7 +1380,7 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si)
              do cl = 1,currentPatch%NCL_p
                 do iv = 1,currentPatch%ncan(cl,ft)
 
-                   if( sum(currentPatch%canopy_area_profile(cl,:,iv)) > 1.0001_r8 ) then
+                   if( DEBUG .and. sum(currentPatch%canopy_area_profile(cl,:,iv)) > 1.0001_r8 ) then
 
                       write(fates_log(), *) 'FATES: A canopy_area_profile exceeded 1.0'
                       write(fates_log(), *) 'cl: ',cl

biogeochem/EDLoggingMortalityMod.F90

@@ -31,6 +31,7 @@ module EDLoggingMortalityMod
    use EDParamsMod       , only : logging_direct_frac
    use EDParamsMod       , only : logging_mechanical_frac 
    use EDParamsMod       , only : logging_coll_under_frac 
+   use EDParamsMod       , only : logging_dbhmax_infra
    use FatesInterfaceMod , only : hlm_current_year
    use FatesInterfaceMod , only : hlm_current_month
    use FatesInterfaceMod , only : hlm_current_day
@@ -154,7 +155,6 @@ subroutine LoggingMortality_frac( pft_i, dbh, lmort_direct,lmort_collateral,lmor
 
       ! Parameters
       real(r8), parameter   :: adjustment = 1.0 ! adjustment for mortality rates
-      real(r8), parameter   :: logging_dbhmax_infra = 35 !(cm), based on Feldpaush et al. (2005) and Ferry et al. (2010)
 
       if (logging_time) then 
          if(EDPftvarcon_inst%woody(pft_i) == 1)then ! only set logging rates for trees

biogeochem/EDMortalityFunctionsMod.F90

@@ -12,7 +12,6 @@ module EDMortalityFunctionsMod
    use EDTypesMod            , only : ed_patch_type
    use FatesConstantsMod     , only : itrue,ifalse
    use FatesAllometryMod     , only : bleaf
-   use EDParamsMod           , only : ED_val_stress_mort
    use FatesInterfaceMod     , only : bc_in_type
    use FatesInterfaceMod     , only : hlm_use_ed_prescribed_phys
    use FatesInterfaceMod     , only : hlm_freq_day
@@ -62,7 +61,6 @@ subroutine mortality_rates( cohort_in,bc_in,cmort,hmort,bmort,frmort )
     real(r8) :: hf_sm_threshold    ! hydraulic failure soil moisture threshold 
     real(r8) :: temp_dep           ! Temp. function (freezing mortality)
     real(r8) :: temp_in_C          ! Daily averaged temperature in Celcius
-    real(r8),parameter :: frost_mort_scaler = 3.0_r8  ! Scaling factor for freezing mortality
     real(r8),parameter :: frost_mort_buffer = 5.0_r8  ! 5deg buffer for freezing mortality
 
     logical, parameter :: test_zero_mortality = .false. ! Developer test which
@@ -79,7 +77,7 @@ subroutine mortality_rates( cohort_in,bc_in,cmort,hmort,bmort,frmort )
     hf_sm_threshold = EDPftvarcon_inst%hf_sm_threshold(cohort_in%pft)
 
     if(cohort_in%patchptr%btran_ft(cohort_in%pft) <= hf_sm_threshold)then 
-       hmort = ED_val_stress_mort
+       hmort = EDPftvarcon_inst%mort_scalar_hydrfailure(cohort_in%pft)
      else
        hmort = 0.0_r8
      endif 
@@ -89,7 +87,8 @@ subroutine mortality_rates( cohort_in,bc_in,cmort,hmort,bmort,frmort )
        call bleaf(cohort_in%dbh,cohort_in%pft,cohort_in%canopy_trim,b_leaf)
        if( b_leaf > 0._r8 .and. cohort_in%bstore <= b_leaf )then
           frac = cohort_in%bstore/ b_leaf
-          cmort = max(0.0_r8,ED_val_stress_mort*(1.0_r8 - frac))
+          cmort = max(0.0_r8,EDPftvarcon_inst%mort_scalar_cstarvation(cohort_in%pft) * &
+               (1.0_r8 - frac))
         else
           cmort = 0.0_r8
        endif
@@ -109,7 +108,7 @@ subroutine mortality_rates( cohort_in,bc_in,cmort,hmort,bmort,frmort )
     temp_in_C = bc_in%t_veg24_si - tfrz
     temp_dep  = max(0.0,min(1.0,1.0 - (temp_in_C - &
                 EDPftvarcon_inst%freezetol(cohort_in%pft))/frost_mort_buffer) )
-    frmort    = frost_mort_scaler * temp_dep
+    frmort    = EDPftvarcon_inst%mort_scalar_coldstress(cohort_in%pft) * temp_dep
 
 
     !mortality_rates = bmort + hmort + cmort

biogeochem/EDPhysiologyMod.F90

@@ -196,9 +196,12 @@ subroutine trim_canopy( currentSite )
           endif
 
           call bleaf(currentcohort%dbh,ipft,currentcohort%canopy_trim,tar_bl)
-          call bfineroot(currentcohort%dbh,ipft,currentcohort%canopy_trim,tar_bfr)
 
-          bfr_per_bleaf = tar_bfr/tar_bl
+          if ( int(EDPftvarcon_inst%allom_fmode(ipft)) .eq. 1 ) then
+             ! only query fine root biomass if using a fine root allometric model that takes leaf trim into account
+             call bfineroot(currentcohort%dbh,ipft,currentcohort%canopy_trim,tar_bfr)
+             bfr_per_bleaf = tar_bfr/tar_bl
+          endif
 
           !Leaf cost vs netuptake for each leaf layer. 
           do z = 1,nlevleaf
@@ -210,18 +213,27 @@ subroutine trim_canopy( currentSite )
 
 
                    currentCohort%leaf_cost =  1._r8/(EDPftvarcon_inst%slatop(ipft)*1000.0_r8)
-                   currentCohort%leaf_cost = currentCohort%leaf_cost + &
-                        1.0_r8/(EDPftvarcon_inst%slatop(ipft)*1000.0_r8) * &
-                        bfr_per_bleaf / EDPftvarcon_inst%root_long(ipft)
+
+                   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.
+                      currentCohort%leaf_cost = currentCohort%leaf_cost + &
+                           1.0_r8/(EDPftvarcon_inst%slatop(ipft)*1000.0_r8) * &
+                           bfr_per_bleaf / EDPftvarcon_inst%root_long(ipft)
+                   endif
 
                    currentCohort%leaf_cost = currentCohort%leaf_cost * &
                          (EDPftvarcon_inst%grperc(ipft) + 1._r8)
                 else !evergreen costs
                    currentCohort%leaf_cost = 1.0_r8/(EDPftvarcon_inst%slatop(ipft)* &
                         EDPftvarcon_inst%leaf_long(ipft)*1000.0_r8) !convert from sla in m2g-1 to m2kg-1
-                   currentCohort%leaf_cost = currentCohort%leaf_cost + &
-                        1.0_r8/(EDPftvarcon_inst%slatop(ipft)*1000.0_r8) * &
-                        bfr_per_bleaf / EDPftvarcon_inst%root_long(ipft)
+                   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.
+                      currentCohort%leaf_cost = currentCohort%leaf_cost + &
+                           1.0_r8/(EDPftvarcon_inst%slatop(ipft)*1000.0_r8) * &
+                           bfr_per_bleaf / EDPftvarcon_inst%root_long(ipft)
+                   endif
                    currentCohort%leaf_cost = currentCohort%leaf_cost * &
                          (EDPftvarcon_inst%grperc(ipft) + 1._r8)
                 endif
@@ -875,15 +887,9 @@ subroutine PlantGrowth( currentSite, currentCohort, bc_in )
     integer , parameter :: max_substeps = 300
     real(r8), parameter :: max_trunc_error = 1.0_r8
     integer,  parameter :: ODESolve = 2    ! 1=RKF45,  2=Euler
-    real(r8), parameter :: global_branch_turnover = 0.0_r8 ! Temporary branch turnover setting
-                                                           ! Branch-turnover control will be 
-                                                           ! introduced in a later PR
-
 
     ipft = currentCohort%pft
 
-    EDPftvarcon_inst%branch_turnover(ipft) = global_branch_turnover
-
     ! Initialize seed production
     currentCohort%seed_prod  = 0.0_r8
 
@@ -976,35 +982,6 @@ subroutine PlantGrowth( currentSite, currentCohort, bc_in )
              currentCohort%canopy_trim, currentCohort%dbh, currentCohort%hite )
     end if
 
-    ! -----------------------------------------------------------------------------------
-    ! III(a). Calculate the maintenance turnover demands 
-    !       Pre-check, make sure phenology is mutually exclusive and at least one chosen
-    !       (MOVE THIS TO THE PARAMETER READ-IN SECTION)
-    ! -----------------------------------------------------------------------------------
-
-    if (EDPftvarcon_inst%evergreen(ipft) == 1) then
-       if (EDPftvarcon_inst%season_decid(ipft) == 1)then 
-          write(fates_log(),*) 'PFT # ',ipft,' was specified as being both evergreen'
-          write(fates_log(),*) '       and seasonally deciduous, impossible, aborting'
-          call endrun(msg=errMsg(sourcefile, __LINE__))
-       end if
-       if (EDPftvarcon_inst%stress_decid(ipft) == 1)then 
-          write(fates_log(),*) 'PFT # ',ipft,' was specified as being both evergreen'
-          write(fates_log(),*) '       and stress deciduous, impossible, aborting'
-          call endrun(msg=errMsg(sourcefile, __LINE__))
-       end if
-    end if
-    if (EDPftvarcon_inst%stress_decid(ipft) /= 1 .and. &
-        EDPftvarcon_inst%season_decid(ipft) /= 1 .and. &
-        EDPftvarcon_inst%evergreen(ipft)    /= 1) then
-       write(fates_log(),*) 'PFT # ',ipft,' must be defined as having one of three'
-       write(fates_log(),*) 'phenology habits, ie == 1'
-       write(fates_log(),*) 'stress_decid: ',EDPftvarcon_inst%stress_decid(ipft)
-       write(fates_log(),*) 'season_decid: ',EDPftvarcon_inst%season_decid(ipft)
-       write(fates_log(),*) 'evergreen: ',EDPftvarcon_inst%evergreen(ipft)
-       call endrun(msg=errMsg(sourcefile, __LINE__))
-    endif
-
 
     ! -----------------------------------------------------------------------------------
     ! III(b). Calculate the maintenance turnover demands 
@@ -1463,7 +1440,7 @@ function AllomCGrowthDeriv(c_pools,c_mask,cbalance,currentCohort) result(dCdx)
         if (dbh <= EDPftvarcon_inst%dbh_repro_threshold(ipft)) then ! cap on leaf biomass
            repro_fraction = EDPftvarcon_inst%seed_alloc(ipft)
         else
-           repro_fraction = EDPftvarcon_inst%seed_alloc(ipft) + EDPftvarcon_inst%clone_alloc(ipft)
+           repro_fraction = EDPftvarcon_inst%seed_alloc(ipft) + EDPftvarcon_inst%seed_alloc_mature(ipft)
         end if
 
         dCdx = 0.0_r8

biogeochem/FatesAllometryMod.F90

@@ -113,11 +113,9 @@ module FatesAllometryMod
   public :: StructureResetOfDH ! Method to set DBH to sync with structure biomass
   public :: CheckIntegratedAllometries
 
-
   logical         , parameter :: verbose_logging = .false.
   character(len=*), parameter :: sourcefile = __FILE__
 
-
   ! If testing b4b with older versions, do not remove sapwood
   ! Our old methods with saldarriaga did not remove sapwood from the
   ! bdead pool.  But newer allometries are providing total agb
@@ -708,14 +706,23 @@ subroutine bfineroot(d,ipft,canopy_trim,bfr,dbfrdd)
     real(r8) :: slascaler
 
     select case(int(EDPftvarcon_inst%allom_fmode(ipft)))
-    case(1) ! "constant proportionality with bleaf"
+    case(1) ! "constant proportionality with TRIMMED target bleaf"
 
        call blmax_allom(d,ipft,blmax,dblmaxdd)
        call bfrmax_const(d,blmax,dblmaxdd,ipft,bfrmax,dbfrmaxdd)
        bfr    = bfrmax * canopy_trim
        if(present(dbfrdd))then
           dbfrdd = dbfrmaxdd * canopy_trim
        end if
+    case(2) ! "constant proportionality with UNTRIMMED target bleaf"
+
+       call blmax_allom(d,ipft,blmax,dblmaxdd)
+       call bfrmax_const(d,blmax,dblmaxdd,ipft,bfrmax,dbfrmaxdd)
+       bfr    = bfrmax
+       if(present(dbfrdd))then
+          dbfrdd = dbfrmaxdd
+       end if
+
     case DEFAULT 
        write(fates_log(),*) 'An undefined fine root allometry was specified: ', &
             EDPftvarcon_inst%allom_fmode(ipft)

biogeophys/FatesPlantRespPhotosynthMod.F90

@@ -191,7 +191,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime)
     ! Ball-Berry minimum leaf conductance, unstressed (umol H2O/m**2/s)
     ! For C3 and C4 plants
     ! -----------------------------------------------------------------------------------
-    real(r8), dimension(2) :: bbbopt 
+    real(r8), dimension(0:1) :: bbbopt 
 
     associate(  &
          c3psn     => EDPftvarcon_inst%c3psn  , &
@@ -202,9 +202,9 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime)
          frootcn   => EDPftvarcon_inst%frootcn, & ! froot C:N (gc/gN)   ! slope of BB relationship
          q10       => FatesSynchronizedParamsInst%Q10 )
 
+      bbbopt(0) = ED_val_bbopt_c4
       bbbopt(1) = ED_val_bbopt_c3
-      bbbopt(2) = ED_val_bbopt_c4
-
+
       do s = 1,nsites
 
          ! Multi-layer parameters scaled by leaf nitrogen profile.
@@ -324,10 +324,12 @@ 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 
 
-                        if(cl==NCL_p)then !are we in the top canopy layer or a shaded layer?
+                        if(cl==1)then !are we in the top canopy layer or a shaded layer?
                            laican = 0._r8
                         else
-                           laican = sum(currentPatch%canopy_layer_tai(cl+1:NCL_p)) 
+
+                           laican = sum(currentPatch%canopy_layer_tai(1:cl-1)) 
+
                         end if
 
                         ! Loop over leaf-layers
@@ -366,7 +368,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime)
                                  laican = laican + 0.5_r8 * vai
                               else
                                  laican = laican + 0.5_r8 * (currentPatch%elai_profile(cl,ft,iv-1)+ &
-                                       currentPatch%esai_profile(cl,ft,iv-1))+vai
+                                       currentPatch%esai_profile(cl,ft,iv-1)+vai)
                               end if
 
                               ! Scale for leaf nitrogen profile
@@ -719,7 +721,8 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl,         &  ! in
 
    ! Locals
    ! ------------------------------------------------------------------------
-   integer :: pp_type            ! Index for the different photosynthetic pathways C3,C4
+   integer :: c3c4_path_index    ! Index for which photosynthetic pathway 
+                                 ! is active.  C4 = 0,  C3 = 1
    integer :: sunsha             ! Index for differentiating sun and shade
    real(r8) :: gstoma            ! Stomatal Conductance of this leaf layer (m/s)
    real(r8) :: agross            ! co-limited gross leaf photosynthesis (umol CO2/m**2/s)
@@ -756,22 +759,23 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl,         &  ! in
    real(r8),parameter :: init_a2l_co2_c3 = 0.7_r8
    real(r8),parameter :: init_a2l_co2_c4 = 0.4_r8
 
-   ! quantum efficiency, used only for C4 (mol CO2 / mol photons)
-   real(r8),parameter,dimension(2) :: quant_eff = [0.0_r8,0.05_r8]
+   ! quantum efficiency, used only for C4 (mol CO2 / mol photons) (index 0)
+   real(r8),parameter,dimension(0:1) :: quant_eff = [0.05_r8,0.0_r8]
 
    ! empirical curvature parameter for ac, aj photosynthesis co-limitation
-   real(r8),parameter,dimension(2) :: theta_cj  = [0.98_r8,0.80_r8]
+   real(r8),parameter,dimension(0:1) :: theta_cj  = [0.80_r8,0.98_r8]
 
    ! empirical curvature parameter for ap photosynthesis co-limitation
    real(r8),parameter :: theta_ip = 0.95_r8
 
-   associate( bb_slope  => EDPftvarcon_inst%BB_slope ) ! slope of BB relationship
+   associate( bb_slope  => EDPftvarcon_inst%BB_slope)    ! slope of BB relationship
+
+     ! photosynthetic pathway: 0. = c4, 1. = c3
+     c3c4_path_index = nint(EDPftvarcon_inst%c3psn(ft))
 
-     if (nint(EDPftvarcon_inst%c3psn(ft)) == 1) then! photosynthetic pathway: 0. = c4, 1. = c3
-        pp_type = 1
+     if (c3c4_path_index == 1) then
         init_co2_intra_c = init_a2l_co2_c3 * can_co2_ppress
      else
-        pp_type = 2
         init_co2_intra_c = init_a2l_co2_c4 * can_co2_ppress
      end if
 
@@ -843,7 +847,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl,         &  ! in
                  co2_intra_c_old = co2_intra_c
 
                  ! Photosynthesis limitation rate calculations 
-                 if (pp_type == 1)then    
+                 if (c3c4_path_index == 1)then    
 
                     ! C3: Rubisco-limited photosynthesis
                     ac = vcmax * max(co2_intra_c-co2_cpoint, 0._r8) / &
@@ -865,14 +869,14 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl,         &  ! in
                     if(sunsha == 1)then !sunlit
                        !guard against /0's in the night.
                        if((laisun_lsl * canopy_area_lsl) > 0.0000000001_r8) then   
-                          aj = quant_eff(pp_type) * parsun_lsl * 4.6_r8
+                          aj = quant_eff(c3c4_path_index) * parsun_lsl * 4.6_r8
                           !convert from per cohort to per m2 of leaf)
                           aj = aj / (laisun_lsl * canopy_area_lsl)
                        else
                           aj = 0._r8
                        end if
                     else
-                       aj = quant_eff(pp_type) * parsha_lsl * 4.6_r8
+                       aj = quant_eff(c3c4_path_index) * parsha_lsl * 4.6_r8
                        aj = aj / (laisha_lsl * canopy_area_lsl)
                     end if
 
@@ -882,7 +886,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl,         &  ! in
                  end if
 
                  ! Gross photosynthesis smoothing calculations. First co-limit ac and aj. Then co-limit ap
-                 aquad = theta_cj(pp_type)
+                 aquad = theta_cj(c3c4_path_index)
                  bquad = -(ac + aj)
                  cquad = ac * aj
                  call quadratic_f (aquad, bquad, cquad, r1, r2)