Skip to content

Commit

Permalink
Merge pull request ESCOMP#5 from ekluzek/master
Browse files Browse the repository at this point in the history 
Bring sci.1.4.0_api.3.0. to NCAR fates-release
  • Loading branch information
@ekluzek
ekluzek authored Dec 21, 2017
2 parents 8832ce7 + a7e5afe commit dbd4285
Show file tree
Hide file tree
Showing 20 changed files with 1,689 additions and 479 deletions.
48 changes: 20 additions & 28 deletions biogeochem/EDCanopyStructureMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -802,59 +802,51 @@ subroutine canopy_spread( currentSite )
! Calculates the spatial spread of tree canopies based on canopy closure.
!
! !USES:
use EDParamsMod , only : ED_val_maxspread, ED_val_minspread
use EDTypesMod , only : AREA
use EDParamsMod, only : ED_val_canopy_closure_thresh
!
! !ARGUMENTS
type (ed_site_type), intent(inout), target :: currentSite
!
! !LOCAL VARIABLES:
type (ed_cohort_type), pointer :: currentCohort
type (ed_patch_type) , pointer :: currentPatch
real(r8) :: arealayer(nlevleaf) ! Amount of canopy in each layer.
real(r8) :: sitelevel_canopyarea ! Amount of canopy in top layer at the site level
real(r8) :: inc ! Arbitrary daily incremental change in canopy area
integer :: z
!----------------------------------------------------------------------

inc = 0.005_r8
inc = 0.05_r8

currentPatch => currentSite%oldest_patch

sitelevel_canopyarea = 0.0_r8
do while (associated(currentPatch))

!calculate canopy area in each canopy storey...
arealayer = 0.0_r8
!calculate canopy area in each patch...
currentCohort => currentPatch%tallest
do while (associated(currentCohort))
currentCohort%c_area = c_area(currentCohort)
if(EDPftvarcon_inst%woody(currentCohort%pft) == 1)then
arealayer(currentCohort%canopy_layer) = arealayer(currentCohort%canopy_layer) + currentCohort%c_area
if(EDPftvarcon_inst%woody(currentCohort%pft) .eq. 1 .and. currentCohort%canopy_layer .eq. 1 ) then
sitelevel_canopyarea = sitelevel_canopyarea + currentCohort%c_area
endif
currentCohort => currentCohort%shorter
enddo

!If the canopy area is approaching closure, squash the tree canopies and make them taller and thinner
do z = 1,nclmax

if(arealayer(z)/currentPatch%area > 0.9_r8)then
currentPatch%spread(z) = currentPatch%spread(z) - inc
else
currentPatch%spread(z) = currentPatch%spread(z) + inc
endif
if(currentPatch%spread(z) >= ED_val_maxspread)then
currentPatch%spread(z) = ED_val_maxspread
endif
if(currentPatch%spread(z) <= ED_val_minspread)then
currentPatch%spread(z) = ED_val_minspread
endif
enddo !z
!write(fates_log(),*) 'spread',currentPatch%spread(1:2)
!currentPatch%spread(:) = ED_val_maxspread
!FIX(RF,033114) spread is off
!write(fates_log(),*) 'canopy_spread',currentPatch%area,currentPatch%spread(1:2)
currentPatch => currentPatch%younger
currentPatch => currentPatch%younger

enddo !currentPatch

!If the canopy area is approaching closure, squash the tree canopies and make them taller and thinner
if( sitelevel_canopyarea/AREA .gt. ED_val_canopy_closure_thresh ) then
currentSite%spread = currentSite%spread - inc
else
currentSite%spread = currentSite%spread + inc
endif

! put within bounds to make sure it stays between 0 and 1
currentSite%spread = max(min(currentSite%spread, 1._r8), 0._r8)

end subroutine canopy_spread


Expand All @@ -869,7 +861,7 @@ subroutine canopy_summarization( nsites, sites, bc_in )
use FatesInterfaceMod , only : bc_in_type
use EDPatchDynamicsMod , only : set_patchno
use EDPatchDynamicsMod , only : set_root_fraction
use EDTypesMod , only : sizetype_class_index
use FatesSizeAgeTypeIndicesMod, only : sizetype_class_index
use EDGrowthFunctionsMod , only : tree_lai, c_area
use EDtypesMod , only : area
use EDPftvarcon , only : EDPftvarcon_inst
Expand Down
42 changes: 37 additions & 5 deletions biogeochem/EDCohortDynamicsMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@ module EDCohortDynamicsMod
use EDTypesMod , only : nclmax
use EDTypesMod , only : ncwd
use EDTypesMod , only : maxCohortsPerPatch
use EDTypesMod , only : sclass_ed,nlevsclass_ed,AREA
use EDTypesMod , only : AREA
use EDTypesMod , only : min_npm2, min_nppatch
use EDTypesMod , only : min_n_safemath
use FatesInterfaceMod , only : hlm_use_planthydro
Expand All @@ -28,7 +28,8 @@ module EDCohortDynamicsMod
use FatesPlantHydraulicsMod, only : initTreeHydStates
use FatesPlantHydraulicsMod, only : InitHydrCohort
use FatesPlantHydraulicsMod, only : DeallocateHydrCohort
use EDTypesMod , only : sizetype_class_index
use FatesSizeAgeTypeIndicesMod, only : sizetype_class_index


! CIME globals
use shr_log_mod , only : errMsg => shr_log_errMsg
Expand Down Expand Up @@ -415,6 +416,11 @@ subroutine nan_cohort(cc_p)
currentCohort%root_md = nan ! root maintenance demand: kgC/indiv/year
currentCohort%carbon_balance = nan ! carbon remaining for growth and storage: kg/indiv/year
currentCohort%dmort = nan ! proportional mortality rate. (year-1)
currentCohort%lmort_logging = nan
currentCohort%lmort_infra = nan
currentCohort%lmort_collateral = nan


currentCohort%seed_prod = nan ! reproduction seed and clonal: KgC/indiv/year
currentCohort%c_area = nan ! areal extent of canopy (m2)
currentCohort%treelai = nan ! lai of tree (total leaf area (m2) / canopy area (m2)
Expand Down Expand Up @@ -488,7 +494,9 @@ subroutine zero_cohort(cc_p)
currentcohort%dmort = 0._r8
currentcohort%gscan = 0._r8
currentcohort%treesai = 0._r8

currentCohort%lmort_logging = 0._r8
currentCohort%lmort_infra = 0._r8
currentCohort%lmort_collateral = 0._r8
! currentCohort%npp_leaf = 0._r8
! currentCohort%npp_froot = 0._r8
! currentCohort%npp_bsw = 0._r8
Expand Down Expand Up @@ -839,6 +847,13 @@ subroutine fuse_cohorts(patchptr, bc_in)

currentCohort%dmort = (currentCohort%n*currentCohort%dmort + &
nextc%n*nextc%dmort)/newn
currentCohort%lmort_logging = (currentCohort%n*currentCohort%lmort_logging + &
nextc%n*nextc%lmort_logging)/newn
currentCohort%lmort_infra = (currentCohort%n*currentCohort%lmort_infra + &
nextc%n*nextc%lmort_infra)/newn
currentCohort%lmort_collateral = (currentCohort%n*currentCohort%lmort_collateral + &
nextc%n*nextc%lmort_collateral)/newn

currentCohort%fire_mort = (currentCohort%n*currentCohort%fire_mort + &
nextc%n*nextc%fire_mort)/newn
currentCohort%leaf_litter = (currentCohort%n*currentCohort%leaf_litter + &
Expand All @@ -851,6 +866,14 @@ subroutine fuse_cohorts(patchptr, bc_in)
currentCohort%imort = (currentCohort%n*currentCohort%imort + nextc%n*nextc%imort)/newn
currentCohort%fmort = (currentCohort%n*currentCohort%fmort + nextc%n*nextc%fmort)/newn

! logging mortality, Yi Xu
currentCohort%lmort_logging = (currentCohort%n*currentCohort%lmort_logging + &
nextc%n*nextc%lmort_logging)/newn
currentCohort%lmort_collateral = (currentCohort%n*currentCohort%lmort_collateral + &
nextc%n*nextc%lmort_collateral)/newn
currentCohort%lmort_infra = (currentCohort%n*currentCohort%lmort_infra + &
nextc%n*nextc%lmort_infra)/newn

! npp diagnostics
currentCohort%npp_leaf = (currentCohort%n*currentCohort%npp_leaf + nextc%n*nextc%npp_leaf) &
/newn
Expand Down Expand Up @@ -940,10 +963,11 @@ subroutine fuse_cohorts(patchptr, bc_in)
write(fates_log(),*) 'maxcohorts exceeded',dynamic_fusion_tolerance

else

iterate = 0
endif
endif

if ( dynamic_fusion_tolerance .gt. 100._r8) then
if ( dynamic_fusion_tolerance .gt. 100._r8) then
! something has gone terribly wrong and we need to report what
write(fates_log(),*) 'exceeded reasonable expectation of cohort fusion.'
currentCohort => currentPatch%tallest
Expand Down Expand Up @@ -1217,6 +1241,9 @@ subroutine copy_cohort( currentCohort,copyc )
n%root_md = o%root_md
n%carbon_balance = o%carbon_balance
n%dmort = o%dmort
n%lmort_logging = o%lmort_logging
n%lmort_infra = o%lmort_infra
n%lmort_collateral= o%lmort_collateral
n%seed_prod = o%seed_prod
n%treelai = o%treelai
n%treesai = o%treesai
Expand All @@ -1231,6 +1258,11 @@ subroutine copy_cohort( currentCohort,copyc )
n%fmort = o%fmort
n%hmort = o%hmort

! logging mortalities, Yi Xu
n%lmort_logging=o%lmort_logging
n%lmort_collateral =o%lmort_collateral
n%lmort_infra =o%lmort_infra

! Flags
n%isnew = o%isnew

Expand Down
56 changes: 31 additions & 25 deletions biogeochem/EDGrowthFunctionsMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@ module EDGrowthFunctionsMod
use FatesGlobals , only : fates_log
use EDPftvarcon , only : EDPftvarcon_inst
use EDTypesMod , only : ed_cohort_type, nlevleaf, dinc_ed
use FatesConstantsMod , only : itrue,ifalse

implicit none
private
Expand Down Expand Up @@ -82,10 +83,10 @@ real(r8) function Hite( cohort_in )
! if the hite is larger than the maximum allowable height (set by dbhmax) then
! set the height to the maximum value.
! this could do with at least re-factoring and probably re-thinking. RF
if(cohort_in%dbh <= EDPftvarcon_inst%max_dbh(cohort_in%pft)) then
if(cohort_in%dbh <= EDPftvarcon_inst%allom_dbh_maxheight(cohort_in%pft)) then
h = (10.0_r8**(log10(cohort_in%dbh) * m + c))
else
h = (10.0_r8**(log10(EDPftvarcon_inst%max_dbh(cohort_in%pft))*m + c))
h = (10.0_r8**(log10(EDPftvarcon_inst%allom_dbh_maxheight(cohort_in%pft))*m + c))
endif
Hite = h

Expand All @@ -106,26 +107,22 @@ real(r8) function Bleaf( cohort_in )
real(r8) :: dbh2bl_a
real(r8) :: dbh2bl_b
real(r8) :: dbh2bl_c
real(r8) :: slascaler ! changes the target biomass according to the SLA

dbh2bl_a = EDPftvarcon_inst%allom_d2bl1(cohort_in%pft)
dbh2bl_b = EDPftvarcon_inst%allom_d2bl2(cohort_in%pft)
dbh2bl_c = EDPftvarcon_inst%allom_d2bl3(cohort_in%pft)
slascaler = EDPftvarcon_inst%allom_d2bl_slascaler(cohort_in%pft)/EDPftvarcon_inst%slatop(cohort_in%pft)

if(cohort_in%dbh < 0._r8.or.cohort_in%pft == 0.or.cohort_in%dbh > 1000.0_r8)then
write(fates_log(),*) 'problems in bleaf',cohort_in%dbh,cohort_in%pft
endif

if(cohort_in%dbh <= EDPftvarcon_inst%max_dbh(cohort_in%pft))then
if(cohort_in%dbh <= EDPftvarcon_inst%allom_dbh_maxheight(cohort_in%pft))then
bleaf = dbh2bl_a * (cohort_in%dbh**dbh2bl_b) * EDPftvarcon_inst%wood_density(cohort_in%pft)**dbh2bl_c
else
bleaf = dbh2bl_a * (EDPftvarcon_inst%max_dbh(cohort_in%pft)**dbh2bl_b) * &
bleaf = dbh2bl_a * (EDPftvarcon_inst%allom_dbh_maxheight(cohort_in%pft)**dbh2bl_b) * &
EDPftvarcon_inst%wood_density(cohort_in%pft)**dbh2bl_c
endif

bleaf = bleaf * slascaler

!write(fates_log(),*) 'bleaf',bleaf, slascaler,cohort_in%pft

!Adjust for canopies that have become so deep that their bottom layer is not producing any carbon...
Expand Down Expand Up @@ -224,14 +221,13 @@ real(r8) function c_area( cohort_in )
! Function of DBH (cm) canopy spread (m/cm) and number of individuals.
! ============================================================================

use EDParamsMod , only : ED_val_grass_spread
use EDTypesMod , only : nclmax

type(ed_cohort_type), intent(in) :: cohort_in

real(r8) :: dbh ! Tree diameter at breat height. cm.
real(r8) :: crown_area_to_dbh_exponent
integer :: can_layer_index
real(r8) :: spreadterm

! default is to use the same exponent as the dbh to bleaf exponent so that per-plant canopy depth remains invariant during growth,
! but allowed to vary via the allom_blca_expnt_diff term (which has default value of zero)
Expand All @@ -240,15 +236,14 @@ real(r8) function c_area( cohort_in )

if (DEBUG_growth) then
write(fates_log(),*) 'z_area 1',cohort_in%dbh,cohort_in%pft
write(fates_log(),*) 'z_area 2',EDPftvarcon_inst%max_dbh
write(fates_log(),*) 'z_area 2',EDPftvarcon_inst%allom_dbh_maxheight
write(fates_log(),*) 'z_area 3',EDPftvarcon_inst%woody
write(fates_log(),*) 'z_area 4',cohort_in%n
write(fates_log(),*) 'z_area 5',cohort_in%patchptr%spread
write(fates_log(),*) 'z_area 5',cohort_in%siteptr%spread
write(fates_log(),*) 'z_area 6',cohort_in%canopy_layer
write(fates_log(),*) 'z_area 7',ED_val_grass_spread
end if

dbh = min(cohort_in%dbh,EDPftvarcon_inst%max_dbh(cohort_in%pft))
dbh = min(cohort_in%dbh,EDPftvarcon_inst%allom_dbh_maxheight(cohort_in%pft))

! ----------------------------------------------------------------------------------
! The function c_area is called during the process of canopy position demotion
Expand All @@ -259,15 +254,12 @@ real(r8) function c_area( cohort_in )
! So we allow layer index exceedence here and force it down to max.
! (rgk/cdk 05/2017)
! ----------------------------------------------------------------------------------

can_layer_index = min(cohort_in%canopy_layer,nclmax)

if(EDPftvarcon_inst%woody(cohort_in%pft) == 1)then
c_area = 3.142_r8 * cohort_in%n * &
(cohort_in%patchptr%spread(can_layer_index)*dbh)**crown_area_to_dbh_exponent
else
c_area = 3.142_r8 * cohort_in%n * (ED_val_grass_spread*dbh)**crown_area_to_dbh_exponent
end if

! apply site-level spread elasticity to the cohort crown allometry term
spreadterm = cohort_in%siteptr%spread * EDPftvarcon_inst%allom_d2ca_coefficient_max(cohort_in%pft) + &
(1._r8 - cohort_in%siteptr%spread) * EDPftvarcon_inst%allom_d2ca_coefficient_min(cohort_in%pft)
!
c_area = cohort_in%n * spreadterm * dbh ** crown_area_to_dbh_exponent

end function c_area

Expand Down Expand Up @@ -364,7 +356,7 @@ real(r8) function dDbhdBd( cohort_in )
dBD_dDBH = dbh2bd_c*dbh2bd_a*(cohort_in%hite**dbh2bd_b)*(cohort_in%dbh**(dbh2bd_c-1.0_r8))* &
(EDPftvarcon_inst%wood_density(cohort_in%pft)**dbh2bd_d)

if(cohort_in%dbh < EDPftvarcon_inst%max_dbh(cohort_in%pft))then
if(cohort_in%dbh < EDPftvarcon_inst%allom_dbh_maxheight(cohort_in%pft))then
dH_dDBH = (10.0_r8**c)*m*(cohort_in%dbh**(m-1.0_r8))

dBD_dDBH = dBD_dDBH + dbh2bd_b*dbh2bd_a*(cohort_in%hite**(dbh2bd_b - 1.0_r8))* &
Expand Down Expand Up @@ -400,7 +392,7 @@ real(r8) function dDbhdBl( cohort_in )
dblddbh = dbh2bl_b*dbh2bl_a*(cohort_in%dbh**dbh2bl_b)*(EDPftvarcon_inst%wood_density(cohort_in%pft)**dbh2bl_c)
dblddbh = dblddbh*cohort_in%canopy_trim

if( cohort_in%dbh<EDPftvarcon_inst%max_dbh(cohort_in%pft) ) then
if( cohort_in%dbh<EDPftvarcon_inst%allom_dbh_maxheight(cohort_in%pft) ) then
dDbhdBl = 1.0_r8/dblddbh
else
dDbhdBl = 1.0d15 ! At maximum size, the leaf biomass is saturated, dbl=0
Expand All @@ -419,6 +411,7 @@ subroutine mortality_rates( cohort_in,cmort,hmort,bmort )
! ============================================================================

use EDParamsMod, only : ED_val_stress_mort
use FatesInterfaceMod, only : hlm_use_ed_prescribed_phys

type (ed_cohort_type), intent(in) :: cohort_in
real(r8),intent(out) :: bmort ! background mortality : Fraction per year
Expand All @@ -429,6 +422,9 @@ subroutine mortality_rates( cohort_in,cmort,hmort,bmort )

real(r8) :: hf_sm_threshold ! hydraulic failure soil moisture threshold


if (hlm_use_ed_prescribed_phys .eq. ifalse) then

! 'Background' mortality (can vary as a function of density as in ED1.0 and ED2.0, but doesn't here for tractability)
bmort = EDPftvarcon_inst%bmort(cohort_in%pft)

Expand Down Expand Up @@ -457,6 +453,16 @@ subroutine mortality_rates( cohort_in,cmort,hmort,bmort )

!mortality_rates = bmort + hmort + cmort

else ! i.e. hlm_use_ed_prescribed_phys is true
if ( cohort_in%canopy_layer .eq. 1) then
bmort = EDPftvarcon_inst%prescribed_mortality_canopy(cohort_in%pft)
else
bmort = EDPftvarcon_inst%prescribed_mortality_understory(cohort_in%pft)
endif
cmort = 0._r8
hmort = 0._r8
endif

end subroutine mortality_rates

! ============================================================================
Expand Down
Loading

0 comments on commit dbd4285

Please sign in to comment.