Skip to content

Commit

Permalink
Minor conflict resolution with master, both modified definition of c3…
Browse files Browse the repository at this point in the history 
…psn and clumping index, (master had re-named its entry in the parameter file.
  • Loading branch information
@rgknox
rgknox committed Apr 17, 2018
2 parents eb2d66e + a29f894 commit d7520d5
Show file tree
Hide file tree
Showing 13 changed files with 2,485 additions and 203 deletions.
29 changes: 8 additions & 21 deletions biogeochem/EDCanopyStructureMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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))

Expand Down Expand Up @@ -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))
Expand Down Expand Up @@ -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
Expand Down
2 changes: 1 addition & 1 deletion biogeochem/EDLoggingMortalityMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
9 changes: 4 additions & 5 deletions biogeochem/EDMortalityFunctionsMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
Expand Down
65 changes: 21 additions & 44 deletions biogeochem/EDPhysiologyMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
13 changes: 10 additions & 3 deletions biogeochem/FatesAllometryMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down
40 changes: 22 additions & 18 deletions biogeophys/FatesPlantRespPhotosynthMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 , &
Expand All @@ -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.
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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) / &
Expand All @@ -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

Expand All @@ -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)
Expand Down
Loading

0 comments on commit d7520d5

Please sign in to comment.