diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index ed56b1fdbc..9c769255a4 100755 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -802,7 +802,8 @@ 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 @@ -810,51 +811,42 @@ subroutine canopy_spread( 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 @@ -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 diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 80e83cc3ca..6a6be90fac 100755 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -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 @@ -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 @@ -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) @@ -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 @@ -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 + & @@ -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 @@ -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 @@ -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 @@ -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 diff --git a/biogeochem/EDGrowthFunctionsMod.F90 b/biogeochem/EDGrowthFunctionsMod.F90 index 47bbb9591c..075c216a33 100755 --- a/biogeochem/EDGrowthFunctionsMod.F90 +++ b/biogeochem/EDGrowthFunctionsMod.F90 @@ -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 @@ -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 @@ -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... @@ -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) @@ -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 @@ -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 @@ -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))* & @@ -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 fates_r8 + use EDTypesMod , only : ed_cohort_type + use EDTypesMod , only : ed_patch_type + use EDTypesMod , only : ncwd + use EDTypesMod , only : ed_site_type + use EDTypesMod , only : ed_resources_management_type + use EDTypesMod , only : dtype_ilog + use EDTypesMod , only : dtype_ifall + use EDTypesMod , only : dtype_ifire + use EDPftvarcon , only : EDPftvarcon_inst + use EDParamsMod , only : logging_event_code + use EDParamsMod , only : logging_dbhmin + use EDParamsMod , only : logging_collateral_frac + use EDParamsMod , only : logging_direct_frac + use EDParamsMod , only : logging_mechanical_frac + use EDParamsMod , only : ED_val_understorey_death + use FatesInterfaceMod , only : hlm_current_year + use FatesInterfaceMod , only : hlm_current_month + use FatesInterfaceMod , only : hlm_current_day + use FatesInterfaceMod , only : hlm_model_day + use FatesInterfaceMod , only : hlm_day_of_year + use FatesInterfaceMod , only : hlm_days_per_year + use FatesInterfaceMod , only : hlm_use_logging + use FatesConstantsMod , only : itrue,ifalse + use FatesGlobals , only : endrun => fates_endrun + use FatesGlobals , only : fates_log + use shr_log_mod , only : errMsg => shr_log_errMsg + + implicit none + private + + logical, protected :: logging_time ! If true, logging should be + ! performed during the current time-step + + character(len=*), parameter, private :: sourcefile = & + __FILE__ + + public :: LoggingMortality_frac + public :: logging_litter_fluxes + public :: logging_time + public :: IsItLoggingTime + +contains + + subroutine IsItLoggingTime(is_master,currentSite) + + ! ------------------------------------------------------------------------------- + ! This subroutine determines if the current dynamics step should enact + ! the logging module. + ! This is done by comparing the current model time to the logging event + ! ids. If there is a match, it is logging time. + ! ------------------------------------------------------------------------------- + + integer, intent(in) :: is_master + type(ed_site_type), intent(inout), target :: currentSite ! site structure + + integer :: icode ! Integer equivalent of the event code (parameter file only allows reals) + integer :: log_date ! Day of month for logging exctracted from event code + integer :: log_month ! Month of year for logging extraced from event code + integer :: log_year ! Year for logging extracted from event code + character(len=64) :: fmt = '(a,i2.2,a,i2.2,a,i4.4)' + + logging_time = .false. + icode = int(logging_event_code) + + if(hlm_use_logging.eq.ifalse) return + + if(icode .eq. 1) then + ! Logging is turned off + logging_time = .false. + + else if(icode .eq. 2) then + ! Logging event on the first step + if( hlm_model_day.eq.1 ) then + logging_time = .true. + end if + + else if(icode .eq. 3) then + ! Logging event every day + logging_time = .true. + + else if(icode .eq. 4) then + ! logging event once a month + if(hlm_current_day.eq.1 ) then + logging_time = .true. + end if + + else if(icode < 0 .and. icode > -366) then + ! Logging event every year on specific day of year + if(hlm_day_of_year .eq. icode ) then + logging_time = .true. + end if + + else if(icode > 10000 ) then + ! Specific Event: YYYYMMDD + log_date = icode - int(100* floor(real(icode)/100)) + log_year = floor(real(icode)/10000) + log_month = floor(real(icode)/100) - log_year*100 + + if( hlm_current_day.eq.log_date .and. & + hlm_current_month.eq.log_month .and. & + hlm_current_year.eq.log_year ) then + logging_time = .true. + end if + else + ! Bad logging event flag + write(fates_log(),*) 'An invalid logging code was specified in fates_params' + write(fates_log(),*) 'Check EDLoggingMortalityMod.F90:IsItLoggingTime()' + write(fates_log(),*) 'for a breakdown of the valid codes and change' + write(fates_log(),*) 'fates_logging_event_code in the file accordingly.' + write(fates_log(),*) 'exiting' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + ! Initialize some site level diagnostics that are calculated for each event + currentSite%resources_management%delta_litter_stock = 0.0_r8 + currentSite%resources_management%delta_biomass_stock = 0.0_r8 + currentSite%resources_management%delta_individual = 0.0_r8 + + if(logging_time .and. (is_master.eq.itrue) ) then + write(fates_log(),fmt) 'Logging Event Enacted on date: ', & + hlm_current_month,'-',hlm_current_day,'-',hlm_current_year + end if + return + end subroutine IsItLoggingTime + + ! ====================================================================================== + + subroutine LoggingMortality_frac( pft_i, dbh, lmort_logging,lmort_collateral,lmort_infra ) + + ! Arguments + integer, intent(in) :: pft_i ! pft index + real(r8), intent(in) :: dbh ! diameter at breast height (cm) + real(r8), intent(out) :: lmort_logging ! direct (harvestable) mortality fraction + real(r8), intent(out) :: lmort_collateral ! collateral damage mortality fraction + real(r8), intent(out) :: lmort_infra ! infrastructure mortality fraction + + ! Parameters + real(r8), parameter :: adjustment = 1.0 ! adjustment for mortality rates + + if (logging_time) then + if(EDPftvarcon_inst%woody(pft_i) == 1)then ! only set logging rates for trees + + ! Pass logging rates to cohort level + + if (dbh >= logging_dbhmin ) then + lmort_logging = logging_direct_frac * adjustment + lmort_collateral = logging_collateral_frac * adjustment + else + lmort_logging = 0.0_r8 + lmort_collateral = 0.0_r8 + end if + + lmort_infra = logging_mechanical_frac * adjustment + !damage rates for size class < & > threshold_size need to be specified seperately + + ! Collateral damage to smaller plants below the direct logging size threshold + ! will be applied via "understory_death" via the disturbance algorithm + + else + lmort_logging = 0.0_r8 + lmort_collateral = 0.0_r8 + lmort_infra = 0.0_r8 + end if + else + lmort_logging = 0.0_r8 + lmort_collateral = 0.0_r8 + lmort_infra = 0.0_r8 + end if + + end subroutine LoggingMortality_frac + + ! ============================================================================ + + subroutine logging_litter_fluxes(currentSite, currentPatch, newPatch, patch_site_areadis) + + ! ------------------------------------------------------------------------------------------- + ! + ! DESCRIPTION: + ! Carbon going from ongoing mortality into CWD pools. + ! This module includes only those fluxes associated with a disturbance generated by logging. + ! Purpose: + ! 1) move logging-associated carbon to CWD and litter pool + ! 2) move the logging trunk from live into product pool + ! 3) generate fluxes used in carbon balance checking + ! E.g,: + ! Remove trunk of logged trees from litter/CWD + ! Add other parts of logged trees and all parts of collaterally and mechanically + ! damaged trees into CWD/litter + ! + ! This routine is only called if logging disturbance is the dominant disturbance. + ! + ! + ! Note: The litter losses due to disturbance in the logging case is almost + ! exactly like the natural tree-fall case. The big differences are that + ! the mortality rates governing the fluxes, follow a different rule set. + ! We also compute an export flux (product) that does not go to litter. + ! + ! Trunk Product Flux: Only usable wood is exported from a site. This is the above-ground + ! portion of the bole, and only boles associated with direct-logging, + ! not inftrastructure or collateral damage mortality. + ! + ! ------------------------------------------------------------------------------------------- + + + !USES: + use SFParamsMod, only : SF_val_cwd_frac + use EDtypesMod, only : area + use EDtypesMod, only : ed_site_type + use EDtypesMod, only : ed_patch_type + use EDtypesMod, only : ed_cohort_type + use EDGrowthFunctionsMod, only : c_area + + ! !ARGUMENTS: + type(ed_site_type) , intent(inout), target :: currentSite + type(ed_patch_type) , intent(inout), target :: currentPatch + type(ed_patch_type) , intent(inout), target :: newPatch + real(r8) , intent(in) :: patch_site_areadis + + !LOCAL VARIABLES: + type(ed_cohort_type), pointer :: currentCohort + real(r8) :: litter_area ! area over which to distribute this litter (m2/site). + real(r8) :: np_mult ! Fraction of the new patch which came from the current patch + real(r8) :: direct_dead ! Mortality count through direct logging + real(r8) :: indirect_dead ! Mortality count through: impacts, infrastructure and collateral damage + real(r8) :: trunk_product_site ! flux of carbon in trunk products exported off site [ kgC/site ] + ! (note we are accumulating over the patch, but scale is site level) + real(r8) :: delta_litter_stock ! flux of carbon in total litter flux [ kgC/site ] + real(r8) :: delta_biomass_stock ! total flux of carbon through mortality (litter+product) [ kgC/site ] + real(r8) :: delta_individual ! change in plant number through mortality [ plants/site ] + real(r8) :: cwd_litter_density ! Component woody biomass transferred through mortality [kgC/m2] + ! (works with canopy_mortality_woody_litter, breaks into CWD partition + ! and converts units to /m2) + real(r8) :: woody_litter ! Woody biomass transferred through mortality [kgC/site] + real(r8) :: leaf_litter ! Leafy biomass transferred through mortality [kgC/site] + real(r8) :: root_litter ! Rooty + storage biomass transferred through mort [kgC/site] + real(r8) :: agb_frac ! local copy of the above ground biomass fraction [fraction] + integer :: p ! pft index + integer :: c ! cwd index + + + ! Zero some site level accumulator diagnsotics + trunk_product_site = 0.0_r8 + delta_litter_stock = 0.0_r8 + delta_biomass_stock = 0.0_r8 + delta_individual = 0.0_r8 + + + currentCohort => currentPatch%shortest + do while(associated(currentCohort)) + p = currentCohort%pft + + + if(currentCohort%canopy_layer == 1)then + direct_dead = currentCohort%n * currentCohort%lmort_logging + indirect_dead = currentCohort%n * & + (currentCohort%lmort_collateral + currentCohort%lmort_infra) + + else + if(EDPftvarcon_inst%woody(currentCohort%pft) == 1)then + direct_dead = 0.0_r8 + indirect_dead = ED_val_understorey_death * currentCohort%n * & + (patch_site_areadis/currentPatch%area) !kgC/site/day + else + ! If the cohort of interest is grass, it will not experience + ! any mortality associated with the logging disturbance + direct_dead = 0.0_r8 + indirect_dead = 0.0_r8 + end if + end if + + agb_frac = EDPftvarcon_inst%allom_agb_frac(currentCohort%pft) + litter_area = currentPatch%area + np_mult = patch_site_areadis/newPatch%area + + ! ---------------------------------------------------------------------------------------- + ! Handle woody litter flux for non-bole components of biomass + ! This litter is distributed between the current and new patches, & + ! not to any other patches. This is really the eventually area of the current patch & + ! (currentPatch%area-patch_site_areadis) +patch_site_areadis... + ! For the new patch, only some fraction of its land area (patch_areadis/np%area) is + ! derived from the current patch, so we need to multiply by patch_areadis/np%area + ! ---------------------------------------------------------------------------------------- + + do c = 1,ncwd-1 + woody_litter = (direct_dead+indirect_dead) * & + (currentCohort%bdead+currentCohort%bsw) + cwd_litter_density = SF_val_CWD_frac(c) * woody_litter / litter_area + + newPatch%cwd_ag(c) = newPatch%cwd_ag(c) + agb_frac * cwd_litter_density * np_mult + currentPatch%cwd_ag(c) = currentPatch%cwd_ag(c) + agb_frac * cwd_litter_density + newPatch%cwd_bg(c) = newPatch%cwd_bg(c) + (1._r8-agb_frac) * cwd_litter_density * np_mult + currentPatch%cwd_bg(c) = currentPatch%cwd_bg(c) + (1._r8-agb_frac) * cwd_litter_density + + ! Diagnostics on fluxes into the AG and BG CWD pools + currentSite%CWD_AG_diagnostic_input_carbonflux(c) = & + currentSite%CWD_AG_diagnostic_input_carbonflux(c) + & + SF_val_CWD_frac(c) * woody_litter * hlm_days_per_year * agb_frac/ AREA + + currentSite%CWD_BG_diagnostic_input_carbonflux(c) = & + currentSite%CWD_BG_diagnostic_input_carbonflux(c) + & + SF_val_CWD_frac(c) * woody_litter * hlm_days_per_year * (1.0_r8 - agb_frac) / AREA + + ! Diagnostic specific to resource management code + delta_litter_stock = delta_litter_stock + woody_litter * SF_val_CWD_frac(c) + + enddo + + ! ---------------------------------------------------------------------------------------- + ! Handle litter flux for the boles of infrastucture and collateral damage mort + ! In this case the boles from direct logging are exported off-site and are not added + ! to the litter pools. That is why we handle this outside the loop above. Only the + ! collateral damange and infrastructure logging is applied to bole litter + ! ---------------------------------------------------------------------------------------- + + woody_litter = indirect_dead * (currentCohort%bdead+currentCohort%bsw) + + cwd_litter_density = SF_val_CWD_frac(ncwd) * woody_litter / litter_area + + newPatch%cwd_ag(ncwd) = newPatch%cwd_ag(ncwd) + agb_frac * cwd_litter_density * np_mult + currentPatch%cwd_ag(ncwd) = currentPatch%cwd_ag(ncwd) + agb_frac * cwd_litter_density + + newPatch%cwd_bg(ncwd) = newPatch%cwd_bg(ncwd) + (1._r8-agb_frac) * cwd_litter_density * np_mult + currentPatch%cwd_bg(ncwd) = currentPatch%cwd_bg(ncwd) + (1._r8-agb_frac) * cwd_litter_density + + currentSite%CWD_AG_diagnostic_input_carbonflux(ncwd) = & + currentSite%CWD_AG_diagnostic_input_carbonflux(ncwd) + & + SF_val_CWD_frac(ncwd) * woody_litter * hlm_days_per_year * agb_frac/ AREA + + currentSite%CWD_BG_diagnostic_input_carbonflux(ncwd) = & + currentSite%CWD_BG_diagnostic_input_carbonflux(ncwd) + & + SF_val_CWD_frac(ncwd) * woody_litter * hlm_days_per_year * (1.0_r8 - agb_frac) / AREA + + delta_litter_stock = delta_litter_stock + woody_litter * SF_val_CWD_frac(ncwd) + + ! ---------------------------------------------------------------------------------------- + ! Handle litter flux for the belowground portion of directly logged boles + ! ---------------------------------------------------------------------------------------- + + woody_litter = direct_dead * (currentCohort%bdead+currentCohort%bsw) + cwd_litter_density = SF_val_CWD_frac(ncwd) * woody_litter / litter_area + + newPatch%cwd_bg(ncwd) = newPatch%cwd_bg(ncwd) + & + (1._r8-agb_frac) * cwd_litter_density * np_mult + + currentPatch%cwd_bg(ncwd) = currentPatch%cwd_bg(ncwd) + & + (1._r8-agb_frac) * cwd_litter_density + + currentSite%CWD_BG_diagnostic_input_carbonflux(ncwd) = & + currentSite%CWD_BG_diagnostic_input_carbonflux(ncwd) + & + SF_val_CWD_frac(ncwd) * woody_litter * hlm_days_per_year * (1.0_r8 - agb_frac) / AREA + + + ! ---------------------------------------------------------------------------------------- + ! Handle harvest (export, flux-out) flux for the above ground boles + ! In this case the boles from direct logging are exported off-site and are not added + ! to the litter pools. That is why we handle this outside the loop above. Only the + ! collateral damange and infrastructure logging is applied to litter + ! + ! Losses to the system as a whole, for C-balancing (kGC/site/day) + ! Site level product, (kgC/site, accumulated over simulation) + ! ---------------------------------------------------------------------------------------- + + trunk_product_site = trunk_product_site + & + SF_val_CWD_frac(ncwd) * agb_frac * direct_dead * (currentCohort%bdead+currentCohort%bsw) + + + ! ---------------------------------------------------------------------------------------- + ! Handle fluxes of leaf, root and storage carbon into litter pools. + ! (none of these are exported) + ! ---------------------------------------------------------------------------------------- + + leaf_litter = (direct_dead+indirect_dead)*currentCohort%bl + root_litter = (direct_dead+indirect_dead)*(currentCohort%br+currentCohort%bstore) + + newPatch%leaf_litter(p) = newPatch%leaf_litter(p) + leaf_litter / litter_area * np_mult + newPatch%root_litter(p) = newPatch%root_litter(p) + root_litter / litter_area * np_mult + + currentPatch%leaf_litter(p) = currentPatch%leaf_litter(p) + leaf_litter / litter_area + currentPatch%root_litter(p) = currentPatch%root_litter(p) + root_litter / litter_area + + ! track as diagnostic fluxes + currentSite%leaf_litter_diagnostic_input_carbonflux(p) = & + currentSite%leaf_litter_diagnostic_input_carbonflux(p) + & + leaf_litter * hlm_days_per_year / AREA + + currentSite%root_litter_diagnostic_input_carbonflux(p) = & + currentSite%root_litter_diagnostic_input_carbonflux(p) + & + root_litter * hlm_days_per_year / AREA + + + ! Logging specific diagnostics + ! ---------------------------------------------------------------------------------------- + + ! Note that litter stock also has terms above in the CWD loop + delta_litter_stock = delta_litter_stock + & + leaf_litter + & + root_litter + + delta_biomass_stock = delta_biomass_stock + & + leaf_litter + & + root_litter + & + (direct_dead+indirect_dead) * (currentCohort%bdead+currentCohort%bsw) + + delta_individual = delta_individual + & + direct_dead + & + indirect_dead + + currentCohort => currentCohort%taller + end do + + ! Update the amount of carbon exported from the site through logging + ! operations. Currently we assume only above-ground portion + ! of the tree bole that experienced "direct" logging is exported + ! This portion is known as "trunk_product_site + + + currentSite%flux_out = currentSite%flux_out + trunk_product_site + + currentSite%resources_management%trunk_product_site = & + currentSite%resources_management%trunk_product_site + & + trunk_product_site + + currentSite%resources_management%delta_litter_stock = & + currentSite%resources_management%delta_litter_stock + & + delta_litter_stock + + currentSite%resources_management%delta_biomass_stock = & + currentSite%resources_management%delta_biomass_stock + & + delta_biomass_stock + + currentSite%resources_management%delta_individual = & + currentSite%resources_management%delta_individual + & + delta_individual + + currentCohort => newPatch%shortest + do while(associated(currentCohort)) + currentCohort%c_area = c_area(currentCohort) + currentCohort => currentCohort%taller + enddo + + end subroutine logging_litter_fluxes + +end module EDLoggingMortalityMod diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 old mode 100755 new mode 100644 index c84bdadf73..8376d1b226 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -13,6 +13,9 @@ module EDPatchDynamicsMod use EDTypesMod , only : min_patch_area use EDTypesMod , only : nclmax use EDTypesMod , only : maxpft + use EDTypesMod , only : dtype_ifall + use EDTypesMod , only : dtype_ilog + use EDTypesMod , only : dtype_ifire use FatesInterfaceMod , only : hlm_use_planthydro use FatesInterfaceMod , only : hlm_numlevgrnd use FatesInterfaceMod , only : hlm_numlevsoil @@ -25,11 +28,15 @@ module EDPatchDynamicsMod use FatesConstantsMod , only : itrue use FatesPlantHydraulicsMod, only : InitHydrCohort use FatesPlantHydraulicsMod, only : DeallocateHydrCohort + use EDLoggingMortalityMod, only : logging_litter_fluxes + use EDLoggingMortalityMod, only : logging_time use EDParamsMod , only : fates_mortality_disturbance_fraction + ! CIME globals use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=) - + use shr_log_mod , only : errMsg => shr_log_errMsg + ! implicit none private @@ -46,6 +53,8 @@ module EDPatchDynamicsMod public :: set_root_fraction private:: fuse_2_patches + character(len=*), parameter, private :: sourcefile = & + __FILE__ ! 10/30/09: Created by Rosie Fisher ! ============================================================================ @@ -60,30 +69,42 @@ subroutine disturbance_rates( site_in) ! and then determines which is the larger at the patch scale (for now, there an only ! be one disturbance type for each timestep. ! all disturbance rates here are per daily timestep. - ! + + ! 2016-2017 + ! Modify to add logging disturbance + ! !USES: use EDGrowthFunctionsMod , only : c_area, mortality_rates - ! + ! loging flux + use EDLoggingMortalityMod , only : LoggingMortality_frac + + ! !ARGUMENTS: type(ed_site_type) , intent(inout), target :: site_in ! ! !LOCAL VARIABLES: type (ed_patch_type) , pointer :: currentPatch type (ed_cohort_type), pointer :: currentCohort + real(r8) :: cmort real(r8) :: bmort real(r8) :: hmort - !--------------------------------------------------------------------- - !MORTALITY - site_in%disturbance_mortality = 0.0_r8 + real(r8) :: lmort_logging + real(r8) :: lmort_collateral + real(r8) :: lmort_infra - currentPatch => site_in%oldest_patch + integer :: threshold_sizeclass + !---------------------------------------------------------------------------------------------- + ! Calculate Mortality Rates (these were previously calculated during growth derivatives) + ! And the same rates in understory plants have already been applied to %dndt + !---------------------------------------------------------------------------------------------- + + currentPatch => site_in%oldest_patch do while (associated(currentPatch)) currentCohort => currentPatch%shortest - do while(associated(currentCohort)) ! Mortality for trees in the understorey. currentCohort%patchptr => currentPatch @@ -99,42 +120,109 @@ subroutine disturbance_rates( site_in) currentCohort%imort = 0.0_r8 ! Impact mortality is always zero except in new patches currentCohort%fmort = 0.0_r8 ! Fire mortality is initialized as zero, but may be changed + call LoggingMortality_frac(currentCohort%pft, currentCohort%dbh, & + lmort_logging,lmort_collateral,lmort_infra ) + + currentCohort%lmort_logging = lmort_logging + currentCohort%lmort_collateral = lmort_collateral + currentCohort%lmort_infra = lmort_infra + + + currentCohort => currentCohort%taller + end do + currentPatch => currentPatch%younger + end do + + ! --------------------------------------------------------------------------------------------- + ! Calculate Disturbance Rates based on the mortality rates just calculated + ! --------------------------------------------------------------------------------------------- + + currentPatch => site_in%oldest_patch + do while (associated(currentPatch)) + + currentPatch%disturbance_rates(dtype_ifall) = 0.0_r8 + currentPatch%disturbance_rates(dtype_ilog) = 0.0_r8 + + currentCohort => currentPatch%shortest + do while(associated(currentCohort)) + if(currentCohort%canopy_layer == 1)then - currentPatch%disturbance_rates(1) = currentPatch%disturbance_rates(1) + & + ! Treefall Disturbance Rate + currentPatch%disturbance_rates(dtype_ifall) = currentPatch%disturbance_rates(dtype_ifall) + & fates_mortality_disturbance_fraction * & min(1.0_r8,currentCohort%dmort)*hlm_freq_day*currentCohort%c_area/currentPatch%area + ! Logging Disturbance Rate + currentPatch%disturbance_rates(dtype_ilog) = currentPatch%disturbance_rates(dtype_ilog) + & + min(1.0_r8, currentCohort%lmort_logging + & + currentCohort%lmort_collateral + & + currentCohort%lmort_infra ) * & + currentCohort%c_area/currentPatch%area + endif - currentCohort => currentCohort%taller - enddo !currentCohort - ! if fires occur at site + ! Fire Disturbance Rate ! Fudge - fires can't burn the whole patch, as this causes /0 errors. ! This is accumulating the daily fires over the whole 30 day patch generation phase. - currentPatch%disturbance_rates(2) = min(0.99_r8,currentPatch%disturbance_rates(2) + currentPatch%frac_burnt) + currentPatch%disturbance_rates(dtype_ifire) = & + min(0.99_r8,currentPatch%disturbance_rates(dtype_ifire) + currentPatch%frac_burnt) - if (currentPatch%disturbance_rates(2) > 0.98_r8)then - write(fates_log(),*) 'very high fire areas',currentPatch%disturbance_rates(2),currentPatch%frac_burnt + if (currentPatch%disturbance_rates(dtype_ifire) > 0.98_r8)then + write(fates_log(),*) 'very high fire areas', & + currentPatch%disturbance_rates(dtype_ifire),currentPatch%frac_burnt endif - !Only use larger of two natural disturbance modes WHY? - if(currentPatch%disturbance_rates(2) > currentPatch%disturbance_rates(1))then ! DISTURBANCE IS FIRE - currentPatch%disturbance_rate = currentPatch%disturbance_rates(2) - ! RGK 02-18-2014 - ! Since treefall mortality is not actually being applied - ! Go through and zero the diagnostic rates + + ! ------------------------------------------------------------------------------------------ + ! Determine which disturbance is dominant, and force mortality diagnostics in the upper + ! canopy to be zero for the non-dominant mode. Note: upper-canopy tree-fall mortality is + ! not always disturbance generating, so when tree-fall mort is non-dominant, make sure + ! to still diagnose and track the non-disturbance rate + ! ------------------------------------------------------------------------------------------ + + + if (currentPatch%disturbance_rates(dtype_ilog) > currentPatch%disturbance_rates(dtype_ifall) .and. & + currentPatch%disturbance_rates(dtype_ilog) > currentPatch%disturbance_rates(dtype_ifire) ) then + + currentPatch%disturbance_rate = currentPatch%disturbance_rates(dtype_ilog) + + ! Update diagnostics currentCohort => currentPatch%shortest do while(associated(currentCohort)) if(currentCohort%canopy_layer == 1)then - currentCohort%cmort=0.0_r8 - currentCohort%hmort=0.0_r8 - currentCohort%bmort=0.0_r8 + currentCohort%fmort = 0.0_r8 + currentCohort%cmort = currentCohort%cmort*(1.0_r8 - fates_mortality_disturbance_fraction) + currentCohort%hmort = currentCohort%hmort*(1.0_r8 - fates_mortality_disturbance_fraction) + currentCohort%bmort = currentCohort%bmort*(1.0_r8 - fates_mortality_disturbance_fraction) + currentCohort%dmort = currentCohort%dmort*(1.0_r8 - fates_mortality_disturbance_fraction) + ! currentCohort%imort will likely exist with logging end if + currentCohort => currentCohort%taller + enddo !currentCohort + + + elseif (currentPatch%disturbance_rates(dtype_ifire) > currentPatch%disturbance_rates(dtype_ifall) .and. & + currentPatch%disturbance_rates(dtype_ifire) > currentPatch%disturbance_rates(dtype_ilog) ) then ! DISTURBANCE IS FIRE + + currentPatch%disturbance_rate = currentPatch%disturbance_rates(dtype_ifire) + ! Update diagnostics, zero non-fire mortality rates + currentCohort => currentPatch%shortest + do while(associated(currentCohort)) + if(currentCohort%canopy_layer == 1)then + currentCohort%cmort = currentCohort%cmort*(1.0_r8 - fates_mortality_disturbance_fraction) + currentCohort%hmort = currentCohort%hmort*(1.0_r8 - fates_mortality_disturbance_fraction) + currentCohort%bmort = currentCohort%bmort*(1.0_r8 - fates_mortality_disturbance_fraction) + currentCohort%dmort = currentCohort%dmort*(1.0_r8 - fates_mortality_disturbance_fraction) + currentCohort%lmort_logging = 0.0_r8 + currentCohort%lmort_collateral = 0.0_r8 + currentCohort%lmort_infra = 0.0_r8 + end if + ! This may be counter-intuitive, but the diagnostic fire-mortality rate ! will stay zero in the patch that undergoes fire, this is because ! the actual cohorts who experience the fire are only those in the @@ -144,29 +232,30 @@ subroutine disturbance_rates( site_in) currentCohort => currentCohort%taller enddo !currentCohort - else - currentPatch%disturbance_rate = currentPatch%disturbance_rates(1) ! DISTURBANCE IS MORTALITY + else ! If fire and loggin are not greater than treefall, just set disturbance rate to tree-fall + ! which is most likely a 0.0 + + currentPatch%disturbance_rate = currentPatch%disturbance_rates(dtype_ifall) + + ! Update diagnostics, zero non-treefall mortality rates + currentCohort => currentPatch%shortest + do while(associated(currentCohort)) + if(currentCohort%canopy_layer == 1)then + currentCohort%lmort_logging = 0.0_r8 + currentCohort%lmort_collateral = 0.0_r8 + currentCohort%lmort_infra = 0.0_r8 + currentCohort%fmort = 0.0_r8 + end if + currentCohort => currentCohort%taller + enddo !currentCohort + + endif - site_in%disturbance_mortality = site_in%disturbance_mortality + & - currentPatch%disturbance_rates(1)*currentPatch%area/area currentPatch => currentPatch%younger enddo !patch loop - ! FIRE - site_in%disturbance_fire = site_in%frac_burnt/AREA - - ! Use largest disturbance mode and ignore the other... This is necessary to - ! have a single type of disturbance and to calculate the survival rates etc... - if (site_in%disturbance_fire > site_in%disturbance_mortality) then - site_in%disturbance_rate = site_in%disturbance_fire - site_in%dist_type = 2 - else - site_in%disturbance_rate = site_in%disturbance_mortality - site_in%dist_type = 1 - endif - end subroutine disturbance_rates ! ============================================================================ @@ -187,9 +276,9 @@ subroutine spawn_patches( currentSite, bc_in) ! ! !USES: - use EDParamsMod , only : ED_val_maxspread, ED_val_understorey_death + use EDParamsMod , only : ED_val_understorey_death use EDCohortDynamicsMod , only : zero_cohort, copy_cohort, terminate_cohorts - + ! ! !ARGUMENTS: type (ed_site_type), intent(inout), target :: currentSite @@ -211,7 +300,6 @@ subroutine spawn_patches( currentSite, bc_in) real(r8) :: leaf_litter_local(maxpft) ! initial value of leaf litter. KgC/m2 real(r8) :: cwd_ag_local(ncwd) ! initial value of above ground coarse woody debris. KgC/m2 real(r8) :: cwd_bg_local(ncwd) ! initial value of below ground coarse woody debris. KgC/m2 - real(r8) :: spread_local(nclmax) ! initial value of canopy spread parameter.no units !--------------------------------------------------------------------- storesmallcohort => null() ! storage of the smallest cohort for insertion routine @@ -239,18 +327,11 @@ subroutine spawn_patches( currentSite, bc_in) cwd_bg_local = 0.0_r8 leaf_litter_local = 0.0_r8 root_litter_local = 0.0_r8 - spread_local(1:nclmax) = ED_val_maxspread age = 0.0_r8 allocate(new_patch) - -! This is called inside "create_patch" -! create_patch must first allocate some vector spaces before -! zero'ing can occur (RGK) -! call zero_patch(new_patch) - call create_patch(currentSite, new_patch, age, site_areadis, & - spread_local, cwd_ag_local, cwd_bg_local, leaf_litter_local, & + cwd_ag_local, cwd_bg_local, leaf_litter_local, & root_litter_local) new_patch%tallest => null() @@ -259,13 +340,26 @@ subroutine spawn_patches( currentSite, bc_in) currentPatch => currentSite%oldest_patch ! loop round all the patches that contribute surviving indivduals and litter pools to the new patch. do while(associated(currentPatch)) - patch_site_areadis = currentPatch%area * currentPatch%disturbance_rate ! how much land is disturbed in this donor patch? - call average_patch_properties(currentPatch, new_patch, patch_site_areadis) ! MAY BE REDUNDANT CALL - if (currentSite%disturbance_mortality > currentSite%disturbance_fire) then !mortality is dominant disturbance - call mortality_litter_fluxes(currentSite, currentPatch, new_patch, patch_site_areadis) - else + ! This is the amount of patch area that is disturbed, and donated by the donor + patch_site_areadis = currentPatch%area * currentPatch%disturbance_rate + + call average_patch_properties(currentPatch, new_patch, patch_site_areadis) + + if (currentPatch%disturbance_rates(dtype_ilog) > currentPatch%disturbance_rates(dtype_ifall) .and. & + currentPatch%disturbance_rates(dtype_ilog) > currentPatch%disturbance_rates(dtype_ifire) ) then + + call logging_litter_fluxes(currentSite, currentPatch, new_patch, patch_site_areadis) + + elseif (currentPatch%disturbance_rates(dtype_ifire) > currentPatch%disturbance_rates(dtype_ifall) .and. & + currentPatch%disturbance_rates(dtype_ifire) > currentPatch%disturbance_rates(dtype_ilog) ) then + call fire_litter_fluxes(currentSite, currentPatch, new_patch, patch_site_areadis) + + else + + call mortality_litter_fluxes(currentSite, currentPatch, new_patch, patch_site_areadis) + endif !INSERT SURVIVORS FROM DISTURBANCE INTO NEW PATCH @@ -285,9 +379,12 @@ subroutine spawn_patches( currentSite, bc_in) nc%canopy_layer = 1 nc%canopy_layer_yesterday = 1._r8 - !mortality is dominant disturbance - if(currentPatch%disturbance_rates(1) > currentPatch%disturbance_rates(2))then + ! treefall mortality is the dominant disturbance + if(currentPatch%disturbance_rates(dtype_ifall) > currentPatch%disturbance_rates(dtype_ifire) .and. & + currentPatch%disturbance_rates(dtype_ifall) > currentPatch%disturbance_rates(dtype_ilog))then + if(currentCohort%canopy_layer == 1)then + ! In the donor patch we are left with fewer trees because the area has decreased ! the plant density for large trees does not actually decrease in the donor patch ! because this is the part of the original patch where no trees have actually fallen @@ -296,22 +393,32 @@ subroutine spawn_patches( currentSite, bc_in) currentCohort%n = currentCohort%n * (1.0_r8 - fates_mortality_disturbance_fraction * & min(1.0_r8,currentCohort%dmort * hlm_freq_day)) - nc%n = 0.0_r8 ! kill all of the trees who caused the disturbance. + nc%n = 0.0_r8 ! kill all of the trees who caused the disturbance. + nc%cmort = nan ! The mortality diagnostics are set to nan because the cohort should dissappear nc%hmort = nan nc%bmort = nan nc%fmort = nan nc%imort = nan + nc%lmort_logging = nan + nc%lmort_collateral = nan + nc%lmort_infra = nan + else ! small trees if(EDPftvarcon_inst%woody(currentCohort%pft) == 1)then - ! Number of trees in the understory of new patch, before we impose impact mortality and survivorship - nc%n = currentCohort%n * patch_site_areadis/currentPatch%area + ! Survivorship of undestory woody plants. Two step process. + ! Step 1: Reduce current number of plants to reflect the change in area. + ! The number density per square are doesn't change, but since the patch is smaller + ! and cohort counts are absolute, reduce this number. + nc%n = currentCohort%n * patch_site_areadis/currentPatch%area + + ! Step 2: Apply survivor ship function based on the understory death fraction ! remaining of understory plants of those that are knocked over by the overstorey trees dying... - nc%n = (1.0_r8 - ED_val_understorey_death) * currentCohort%n * patch_site_areadis/currentPatch%area - + nc%n = nc%n * (1.0_r8 - ED_val_understorey_death) + ! since the donor patch split and sent a fraction of its members ! to the new patch and a fraction to be preserved in itself, ! when reporting diagnostic rates, we must carry over the mortality rates from @@ -320,11 +427,15 @@ subroutine spawn_patches( currentSite, bc_in) ! number density in EDCLMLink, and the number density of this new patch is donated ! so with the number density must come the effective mortality rates. - nc%fmort = 0.0_r8 ! Should had also been zero in the donor - nc%imort = ED_val_understorey_death/hlm_freq_day ! This was zero in the donor - nc%cmort = currentCohort%cmort - nc%hmort = currentCohort%hmort - nc%bmort = currentCohort%bmort + nc%fmort = 0.0_r8 ! Should had also been zero in the donor + nc%imort = ED_val_understorey_death/hlm_freq_day ! This was zero in the donor + nc%cmort = currentCohort%cmort + nc%hmort = currentCohort%hmort + nc%bmort = currentCohort%bmort + nc%dmort = currentCohort%dmort + nc%lmort_logging = currentCohort%lmort_logging + nc%lmort_collateral = currentCohort%lmort_collateral + nc%lmort_infra = currentCohort%lmort_infra ! understory trees that might potentially be knocked over in the disturbance. ! The existing (donor) patch should not have any impact mortality, it should @@ -340,15 +451,22 @@ subroutine spawn_patches( currentSite, bc_in) ! Those remaining in the existing currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area) - nc%fmort = nan ! These should not make it to the diagnostics - nc%imort = nan ! If they do.. they should invalidate it - nc%cmort = nan ! - nc%hmort = nan ! - nc%bmort = nan ! - + nc%fmort = 0.0_r8 + nc%imort = 0.0_r8 + nc%cmort = currentCohort%cmort + nc%hmort = currentCohort%hmort + nc%bmort = currentCohort%bmort + nc%dmort = currentCohort%dmort + nc%lmort_logging = currentCohort%lmort_logging + nc%lmort_collateral = currentCohort%lmort_collateral + nc%lmort_infra = currentCohort%lmort_infra + endif endif - else !fire + + ! Fire is the dominant disturbance + elseif (currentPatch%disturbance_rates(dtype_ifire) > currentPatch%disturbance_rates(dtype_ifall) .and. & + currentPatch%disturbance_rates(dtype_ifire) > currentPatch%disturbance_rates(dtype_ilog)) then !fire ! Number of members in the new patch, before we impose fire survivorship nc%n = currentCohort%n * patch_site_areadis/currentPatch%area @@ -359,13 +477,103 @@ subroutine spawn_patches( currentSite, bc_in) ! loss of individual from fire in new patch. nc%n = nc%n * (1.0_r8 - currentCohort%fire_mort) - nc%fmort = currentCohort%fire_mort/hlm_freq_day - nc%imort = 0.0_r8 - nc%cmort = currentCohort%cmort - nc%hmort = currentCohort%hmort - nc%bmort = currentCohort%bmort + nc%fmort = currentCohort%fire_mort/hlm_freq_day + nc%imort = 0.0_r8 + + nc%cmort = currentCohort%cmort + nc%hmort = currentCohort%hmort + nc%bmort = currentCohort%bmort + nc%dmort = currentCohort%dmort + nc%lmort_logging = currentCohort%lmort_logging + nc%lmort_collateral = currentCohort%lmort_collateral + nc%lmort_infra = currentCohort%lmort_infra + + ! Logging is the dominant disturbance + elseif (currentPatch%disturbance_rates(dtype_ilog) > currentPatch%disturbance_rates(dtype_ifall) .and. & + currentPatch%disturbance_rates(dtype_ilog) > currentPatch%disturbance_rates(dtype_ifire)) then ! Logging + + ! If this cohort is in the upper canopy. It generated + if(currentCohort%canopy_layer == 1)then + + ! Trees generating this disturbance are not there by definition + nc%n = 0.0_r8 + + ! Reduce counts in the existing/donor patch according to the logging rate + currentCohort%n = currentCohort%n * (1.0_r8 - min(1.0_r8,(currentCohort%lmort_logging + & + currentCohort%lmort_collateral + & + currentCohort%lmort_infra))) + + ! The mortality diagnostics are set to nan because the cohort should dissappear + nc%cmort = nan + nc%hmort = nan + nc%bmort = nan + nc%fmort = nan + nc%imort = nan + nc%lmort_logging = nan + nc%lmort_collateral = nan + nc%lmort_infra = nan + + else + + ! WHat to do with cohorts in the understory of a logging generated + ! disturbance patch? - endif + if(EDPftvarcon_inst%woody(currentCohort%pft) == 1)then + + + ! Survivorship of undestory woody plants. Two step process. + ! Step 1: Reduce current number of plants to reflect the change in area. + ! The number density per square are doesn't change, but since the patch is smaller + ! and cohort counts are absolute, reduce this number. + nc%n = currentCohort%n * patch_site_areadis/currentPatch%area + + ! Step 2: Apply survivor ship function based on the understory death fraction + + ! remaining of understory plants of those that are knocked over by the overstorey trees dying... + ! CURRENTLY ASSUMING THAT LOGGING SURVIVORSHIP OF UNDERSTORY PLANTS IS SAME AS NATURAL + ! TREEFALL (STILL BEING DISCUSSED) + nc%n = nc%n * (1.0_r8 - ED_val_understorey_death) + + ! Step 3: Reduce the number count of cohorts in the original/donor/non-disturbed patch + ! to reflect the area change + currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area) + + + nc%fmort = 0.0_r8 + nc%imort = ED_val_understorey_death/hlm_freq_day + nc%cmort = currentCohort%cmort + nc%hmort = currentCohort%hmort + nc%bmort = currentCohort%bmort + nc%dmort = currentCohort%dmort + nc%lmort_logging = currentCohort%lmort_logging + nc%lmort_collateral = currentCohort%lmort_collateral + nc%lmort_infra = currentCohort%lmort_infra + + else + + ! grass is not killed by mortality disturbance events. Just move it into the new patch area. + ! Just split the grass into the existing and new patch structures + nc%n = currentCohort%n * patch_site_areadis/currentPatch%area + + ! Those remaining in the existing + currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area) + + ! No grass impact mortality imposed on the newly created patch + nc%fmort = 0.0_r8 + nc%imort = 0.0_r8 + nc%cmort = currentCohort%cmort + nc%hmort = currentCohort%hmort + nc%bmort = currentCohort%bmort + nc%dmort = currentCohort%dmort + nc%lmort_logging = currentCohort%lmort_logging + nc%lmort_collateral = currentCohort%lmort_collateral + nc%lmort_infra = currentCohort%lmort_infra + + endif ! is/is-not woody + + endif ! Select canopy layer + + end if ! Select disturbance mode if (nc%n > 0.0_r8) then storebigcohort => new_patch%tallest @@ -419,6 +627,7 @@ subroutine spawn_patches( currentSite, bc_in) enddo ! currentPatch patch loop. + !*************************/ !** INSERT NEW PATCH INTO LINKED LIST !**********`***************/ @@ -539,8 +748,6 @@ subroutine average_patch_properties( currentPatch, newPatch, patch_site_areadis enddo - newPatch%spread = newPatch%spread + currentPatch%spread * patch_site_areadis/newPatch%area - end subroutine average_patch_properties ! ============================================================================ @@ -780,7 +987,7 @@ subroutine mortality_litter_fluxes(currentSite, cp_target, new_patch_target, pat currentCohort => currentPatch%shortest do while(associated(currentCohort)) p = currentCohort%pft - if(currentPatch%disturbance_rates(1) > currentPatch%disturbance_rates(2))then !mortality is dominant disturbance + if(currentCohort%canopy_layer == 1)then !currentCohort%dmort = mortality_rates(currentCohort) !the disturbance calculations are done with the previous n, c_area and d_mort. So it's probably & @@ -812,7 +1019,7 @@ subroutine mortality_litter_fluxes(currentSite, cp_target, new_patch_target, pat ! no-op endif endif - endif + currentCohort => currentCohort%taller @@ -870,7 +1077,7 @@ subroutine mortality_litter_fluxes(currentSite, cp_target, new_patch_target, pat end subroutine mortality_litter_fluxes ! ============================================================================ - subroutine create_patch(currentSite, new_patch, age, areap, spread_local,cwd_ag_local,cwd_bg_local, & + subroutine create_patch(currentSite, new_patch, age, areap,cwd_ag_local,cwd_bg_local, & leaf_litter_local,root_litter_local) ! ! !DESCRIPTION: @@ -887,7 +1094,6 @@ subroutine create_patch(currentSite, new_patch, age, areap, spread_local,cwd_ag_ real(r8), intent(in) :: cwd_bg_local(:) ! initial value of below ground coarse woody debris. KgC/m2 real(r8), intent(in) :: root_litter_local(:)! initial value of root litter. KgC/m2 real(r8), intent(in) :: leaf_litter_local(:)! initial value of leaf litter. KgC/m2 - real(r8), intent(in) :: spread_local(:) ! initial value of canopy spread parameter.no units ! ! !LOCAL VARIABLES: !--------------------------------------------------------------------- @@ -917,7 +1123,6 @@ subroutine create_patch(currentSite, new_patch, age, areap, spread_local,cwd_ag_ new_patch%age = age new_patch%age_class = 1 new_patch%area = areap - new_patch%spread = spread_local new_patch%cwd_ag = cwd_ag_local new_patch%cwd_bg = cwd_bg_local new_patch%leaf_litter = leaf_litter_local @@ -1012,7 +1217,6 @@ subroutine zero_patch(cp_p) currentPatch%nrad(:,:) = 999 ! number of exposed leaf layers for each canopy layer and pft currentPatch%ncan(:,:) = 999 ! number of total leaf layers for each canopy layer and pft currentPatch%lai = nan ! leaf area index of patch - currentPatch%spread(:) = nan ! dynamic ratio of dbh to canopy area. currentPatch%pft_agb_profile(:,:) = nan ! DISTURBANCE @@ -1074,7 +1278,6 @@ subroutine zero_patch(cp_p) currentPatch%sabs_dif(:) = 0.0_r8 currentPatch%zstar = 0.0_r8 - end subroutine zero_patch ! ============================================================================ @@ -1236,7 +1439,7 @@ subroutine fuse_2_patches(dp, rp) ! associated with the secnd patch ! ! !USES: - use EDTypesMod, only: get_age_class_index + use FatesSizeAgeTypeIndicesMod, only: get_age_class_index ! ! !ARGUMENTS: type (ed_patch_type) , intent(inout), pointer :: dp ! Donor Patch diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index dcdcdf437d..7bc7f1a393 100755 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -31,6 +31,7 @@ module EDPhysiologyMod use FatesGlobals , only : fates_log use FatesGlobals , only : endrun => fates_endrun use EDParamsMod , only : fates_mortality_disturbance_fraction + use FatesConstantsMod , only : itrue,ifalse implicit none private @@ -124,7 +125,7 @@ subroutine non_canopy_derivs( currentSite, currentPatch, bc_in ) call seed_germination(currentSite, currentPatch) ! update fragmenting pool fluxes - call cwd_input(currentPatch) + call cwd_input( currentSite, currentPatch) call cwd_out( currentSite, currentPatch, bc_in) do p = 1,numpft @@ -762,6 +763,8 @@ subroutine Growth_Derivatives( currentSite, currentCohort, bc_in) ! ! !USES: use EDGrowthFunctionsMod , only : Bleaf, dDbhdBd, dhdbd, hite, mortality_rates,dDbhdBl + use FatesInterfaceMod, only : hlm_use_ed_prescribed_phys + use EDLoggingMortalityMod, only : LoggingMortality_frac ! ! !ARGUMENTS @@ -784,14 +787,31 @@ subroutine Growth_Derivatives( currentSite, currentCohort, bc_in) real(r8) :: cmort ! starvation mortality rate (fraction per year) real(r8) :: bmort ! background mortality rate (fraction per year) real(r8) :: hmort ! hydraulic failure mortality rate (fraction per year) + + real(r8) :: lmort_logging ! Mortality fraction associated with direct logging + real(r8) :: lmort_collateral ! Mortality fraction associated with logging collateral damage + real(r8) :: lmort_infra ! Mortality fraction associated with logging infrastructure + real(r8) :: dndt_logging ! Mortality rate (per day) associated with the a logging event + real(r8) :: balive_loss !---------------------------------------------------------------------- ! Mortality for trees in the understorey. !if trees are in the canopy, then their death is 'disturbance'. This probably needs a different terminology call mortality_rates(currentCohort,cmort,hmort,bmort) + call LoggingMortality_frac(currentCohort%pft, currentCohort%dbh, & + currentCohort%lmort_logging, & + currentCohort%lmort_collateral, & + currentCohort%lmort_infra ) + if (currentCohort%canopy_layer > 1)then - currentCohort%dndt = -1.0_r8 * (cmort+hmort+bmort) * currentCohort%n + + ! Include understory logging mortality rates not associated with disturbance + dndt_logging = (currentCohort%lmort_logging + & + currentCohort%lmort_collateral + & + currentCohort%lmort_infra)/hlm_freq_day + + currentCohort%dndt = -1.0_r8 * (cmort+hmort+bmort+dndt_logging) * currentCohort%n else currentCohort%dndt = -(1.0_r8 - fates_mortality_disturbance_fraction) & * (cmort+hmort+bmort) * currentCohort%n @@ -815,12 +835,22 @@ subroutine Growth_Derivatives( currentSite, currentCohort, bc_in) ! NPP if ( DEBUG ) write(fates_log(),*) 'EDphys 716 ',currentCohort%npp_acc - ! convert from kgC/indiv/day into kgC/indiv/year + ! convert from kgC/indiv/day into kgC/indiv/year ! TODO: CONVERT DAYS_PER_YEAR TO DBLE (HOLDING FOR B4B COMPARISONS, RGK-01-2017) currentCohort%npp_acc_hold = currentCohort%npp_acc * hlm_days_per_year currentCohort%gpp_acc_hold = currentCohort%gpp_acc * hlm_days_per_year currentCohort%resp_acc_hold = currentCohort%resp_acc * hlm_days_per_year + if (hlm_use_ed_prescribed_phys .eq. itrue) then + if (currentCohort%canopy_layer .eq. 1) then + currentCohort%npp_acc_hold = EDPftvarcon_inst%prescribed_npp_canopy(currentCohort%pft) & + * currentCohort%c_area / currentCohort%n + else + currentCohort%npp_acc_hold = EDPftvarcon_inst%prescribed_npp_understory(currentCohort%pft) & + * currentCohort%c_area / currentCohort%n + endif + endif + currentSite%flux_in = currentSite%flux_in + currentCohort%npp_acc * currentCohort%n ! Maintenance demands @@ -944,7 +974,7 @@ subroutine Growth_Derivatives( currentSite, currentCohort, bc_in) !only if carbon balance is +ve if ((currentCohort%balive >= target_balive).AND.(currentCohort%carbon_balance > 0._r8))then ! fraction of carbon going into active vs structural carbon - if (currentCohort%dbh <= EDPftvarcon_inst%max_dbh(currentCohort%pft))then ! cap on leaf biomass + if (currentCohort%dbh <= EDPftvarcon_inst%allom_dbh_maxheight(currentCohort%pft))then ! cap on leaf biomass dbldbd = dDbhdBd(currentCohort)/dDbhdBl(currentCohort) dbrdbd = EDPftvarcon_inst%allom_l2fr(currentCohort%pft) * dbldbd dhdbd_fn = dhdbd(currentCohort) @@ -980,6 +1010,7 @@ subroutine Growth_Derivatives( currentSite, currentCohort, bc_in) if ( DEBUG ) write(fates_log(),*) 'EDPhys dbstoredt I ',currentCohort%dbstoredt currentCohort%seed_prod = (1.0_r8 - gr_fract) * currentCohort%carbon_balance + if (abs(currentCohort%npp_acc_hold-(currentCohort%dbalivedt+currentCohort%dbdeaddt+currentCohort%dbstoredt+ & currentCohort%seed_prod+currentCohort%md)) > 0.0000000001_r8)then write(fates_log(),*) 'error in carbon check growth derivs',currentCohort%npp_acc_hold- & @@ -1015,16 +1046,16 @@ subroutine Growth_Derivatives( currentSite, currentCohort, bc_in) end subroutine Growth_Derivatives ! ============================================================================ - subroutine recruitment( t, currentSite, currentPatch, bc_in ) + subroutine recruitment( currentSite, currentPatch, bc_in ) ! ! !DESCRIPTION: ! spawn new cohorts of juveniles of each PFT ! ! !USES: use EDGrowthFunctionsMod, only : bdead,dbh, Bleaf + use FatesInterfaceMod, only : hlm_use_ed_prescribed_phys ! ! !ARGUMENTS - integer, intent(in) :: t type(ed_site_type), intent(inout), target :: currentSite type(ed_patch_type), intent(inout), pointer :: currentPatch type(bc_in_type), intent(in) :: bc_in @@ -1049,14 +1080,13 @@ subroutine recruitment( t, currentSite, currentPatch, bc_in ) + EDpftvarcon_inst%allom_latosa_int(ft)*temp_cohort%hite) temp_cohort%bstore = EDPftvarcon_inst%cushion(ft)*(temp_cohort%balive/ (1.0_r8 + EDPftvarcon_inst%allom_l2fr(ft) & + EDpftvarcon_inst%allom_latosa_int(ft)*temp_cohort%hite)) - temp_cohort%n = currentPatch%area * currentPatch%seed_germination(ft)*hlm_freq_day & - / (temp_cohort%bdead+temp_cohort%balive+temp_cohort%bstore) - - if (t == 1)then - write(fates_log(),*) 'filling in cohorts where there are none left; this will break carbon balance', & - currentPatch%patchno,currentPatch%area - temp_cohort%n = 0.1_r8*currentPatch%area - write(fates_log(),*) 'cohort n',ft,temp_cohort%n + + if (hlm_use_ed_prescribed_phys .eq. ifalse) then + temp_cohort%n = currentPatch%area * currentPatch%seed_germination(ft)*hlm_freq_day & + / (temp_cohort%bdead+temp_cohort%balive+temp_cohort%bstore) + else + ! prescribed recruitment rates. number per sq. meter per year + temp_cohort%n = currentPatch%area * EDPftvarcon_inst%prescribed_recruitment(ft) * hlm_freq_day endif temp_cohort%laimemory = 0.0_r8 @@ -1093,7 +1123,7 @@ subroutine recruitment( t, currentSite, currentPatch, bc_in ) end subroutine recruitment ! ============================================================================ - subroutine CWD_Input( currentPatch) + subroutine CWD_Input( currentSite, currentPatch) ! ! !DESCRIPTION: ! Generate litter fields from turnover. @@ -1103,13 +1133,18 @@ subroutine CWD_Input( currentPatch) ! ! !ARGUMENTS + type(ed_site_type), intent(inout), target :: currentSite type(ed_patch_type),intent(inout), target :: currentPatch ! ! !LOCAL VARIABLES: type(ed_cohort_type), pointer :: currentCohort integer :: c,p - real(r8) :: not_dead_n !projected remaining number of trees in understorey cohort after turnover - real(r8) :: dead_n !understorey dead tree density + real(r8) :: dead_n ! total understorey dead tree density + real(r8) :: dead_n_dlogging ! direct logging understory dead-tree density + real(r8) :: dead_n_ilogging ! indirect understory dead-tree density (logging) + real(r8) :: dead_n_natural ! understory dead density not associated + ! with direct logging + real(r8) :: trunk_product ! carbon flux into trunk products kgC/day/site integer :: pft !---------------------------------------------------------------------- @@ -1146,29 +1181,103 @@ subroutine CWD_Input( currentPatch) ! ================================================ ! Litter fluxes for understorey mortality. KgC/m2/year ! ================================================ + + ! Total number of dead understory (n/m2) dead_n = -1.0_r8 * currentCohort%dndt / currentPatch%area + ! Total number of dead understory from direct logging + ! (it is possible that large harvestable trees are in the understory) + dead_n_dlogging = ( currentCohort%lmort_logging) * & + currentCohort%n/hlm_freq_day/currentPatch%area + + ! Total number of dead understory from indirect logging + dead_n_ilogging = ( currentCohort%lmort_collateral + currentCohort%lmort_infra) * & + currentCohort%n/hlm_freq_day/currentPatch%area + + dead_n_natural = dead_n - dead_n_dlogging - dead_n_ilogging + + currentPatch%leaf_litter_in(pft) = currentPatch%leaf_litter_in(pft) + & - (currentCohort%bl+currentCohort%leaf_litter/hlm_freq_day)* dead_n + currentCohort%bl * dead_n currentPatch%root_litter_in(pft) = currentPatch%root_litter_in(pft) + & (currentCohort%br+currentCohort%bstore) * dead_n + ! Update diagnostics that track resource management + currentSite%resources_management%delta_litter_stock = & + currentSite%resources_management%delta_litter_stock + & + (currentCohort%bl+currentCohort%br+currentCohort%bstore) * & + (dead_n_ilogging+dead_n_dlogging) * & + hlm_freq_day * currentPatch%area + ! Update diagnostics that track resource management + currentSite%resources_management%delta_biomass_stock = & + currentSite%resources_management%delta_biomass_stock + & + (currentCohort%bl+currentCohort%br+currentCohort%bstore) * & + (dead_n_ilogging+dead_n_dlogging) * & + hlm_freq_day * currentPatch%area + do c = 1,ncwd - currentPatch%cwd_AG_in(c) = currentPatch%cwd_AG_in(c) + (currentCohort%bdead+currentCohort%bsw) * & - SF_val_CWD_frac(c) * dead_n * EDPftvarcon_inst%allom_agb_frac(currentCohort%pft) + currentPatch%cwd_BG_in(c) = currentPatch%cwd_BG_in(c) + (currentCohort%bdead+currentCohort%bsw) * & - SF_val_CWD_frac(c) * dead_n * (1.0_r8-EDPftvarcon_inst%allom_agb_frac(currentCohort%pft)) + SF_val_CWD_frac(c) * dead_n * (1.0_r8-EDPftvarcon_inst%allom_agb_frac(currentCohort%pft)) + + ! Send AGB component of boles from non direct-logging activities to AGB litter pool + if (c==ncwd) then + + currentPatch%cwd_AG_in(c) = currentPatch%cwd_AG_in(c) + (currentCohort%bdead+currentCohort%bsw) * & + SF_val_CWD_frac(c) * (dead_n_natural+dead_n_ilogging) * & + EDPftvarcon_inst%allom_agb_frac(currentCohort%pft) + + else + currentPatch%cwd_AG_in(c) = currentPatch%cwd_AG_in(c) + (currentCohort%bdead+currentCohort%bsw) * & + SF_val_CWD_frac(c) * dead_n * & + EDPftvarcon_inst%allom_agb_frac(currentCohort%pft) + + ! Send AGB component of boles from direct-logging activities to export/harvest pool + ! Generate trunk product (kgC/day/site) + trunk_product = (currentCohort%bdead+currentCohort%bsw) * & + SF_val_CWD_frac(c) * dead_n_dlogging * EDPftvarcon_inst%allom_agb_frac(currentCohort%pft) * & + hlm_freq_day * currentPatch%area + + currentSite%flux_out = currentSite%flux_out + trunk_product + + ! Update diagnostics that track resource management + currentSite%resources_management%trunk_product_site = & + currentSite%resources_management%trunk_product_site + & + trunk_product + ! Update diagnostics that track resource management + currentSite%resources_management%trunk_product_site = & + currentSite%resources_management%trunk_product_site + & + trunk_product + end if + + ! Update diagnostics that track resource management + currentSite%resources_management%delta_litter_stock = & + currentSite%resources_management%delta_litter_stock + & + (currentCohort%bdead+currentCohort%bsw) * & + SF_val_CWD_frac(c) * (dead_n_natural+dead_n_ilogging) * & + hlm_freq_day * currentPatch%area + ! Update diagnostics that track resource management + currentSite%resources_management%delta_biomass_stock = & + currentSite%resources_management%delta_biomass_stock + & + (currentCohort%bdead+currentCohort%bsw) * & + SF_val_CWD_frac(c) * dead_n * & + hlm_freq_day * currentPatch%area + if (currentPatch%cwd_AG_in(c) < 0.0_r8)then write(fates_log(),*) 'negative CWD in flux',currentPatch%cwd_AG_in(c), & - (currentCohort%bdead+currentCohort%bsw), dead_n + (currentCohort%bdead+currentCohort%bsw), dead_n endif - enddo + end do + ! Update diagnostics that track resource management + currentSite%resources_management%delta_individual = & + currentSite%resources_management%delta_individual + & + (dead_n_dlogging+dead_n_ilogging) * hlm_freq_day * currentPatch%area + endif !canopy layer - + currentCohort => currentCohort%taller - enddo ! end loop over cohorts do p = 1,numpft diff --git a/main/ChecksBalancesMod.F90 b/main/ChecksBalancesMod.F90 index 7a9e3e7467..49a0351fd0 100644 --- a/main/ChecksBalancesMod.F90 +++ b/main/ChecksBalancesMod.F90 @@ -1,13 +1,16 @@ module ChecksBalancesMod - use shr_kind_mod , only : r8 => shr_kind_r8 + use shr_kind_mod , only : r8 => shr_kind_r8 use shr_const_mod, only: SHR_CONST_CDAY - + use EDtypesMod , only : ed_site_type,ed_patch_type,ed_cohort_type + use EDTypesMod , only : AREA + implicit none private public :: SummarizeNetFluxes public :: FATES_BGC_Carbon_Balancecheck + public :: SiteCarbonStock contains @@ -24,7 +27,7 @@ subroutine SummarizeNetFluxes( nsites, sites, bc_in, is_beg_day ) ! ! !USES: use FatesInterfaceMod , only : bc_in_type - use EDtypesMod , only : ed_site_type,ed_patch_type,ed_cohort_type + use EDtypesMod , only : AREA ! implicit none @@ -234,5 +237,46 @@ subroutine FATES_BGC_Carbon_Balancecheck(nsites, sites, bc_in, is_beg_day, dtime return end subroutine FATES_BGC_Carbon_Balancecheck + ! ============================================================================================== + + subroutine SiteCarbonStock(currentSite,total_stock,biomass_stock,litter_stock,seed_stock) + + type(ed_site_type),intent(inout),target :: currentSite + real(r8),intent(out) :: total_stock + real(r8),intent(out) :: litter_stock + real(r8),intent(out) :: biomass_stock + real(r8),intent(out) :: seed_stock + + type(ed_patch_type), pointer :: currentPatch + type(ed_cohort_type), pointer :: currentCohort + + litter_stock = 0.0_r8 + biomass_stock = 0.0_r8 + seed_stock = sum(currentSite%seed_bank)*AREA + + currentPatch => currentSite%oldest_patch + do while(associated(currentPatch)) + litter_stock = litter_stock + currentPatch%area * & + (sum(currentPatch%cwd_ag) + & + sum(currentPatch%cwd_bg) + & + sum(currentPatch%leaf_litter) + & + sum(currentPatch%root_litter)) + + currentCohort => currentPatch%tallest + do while(associated(currentCohort)) + biomass_stock = biomass_stock + (currentCohort%bdead + currentCohort%balive + & + currentCohort%bstore) * currentCohort%n + currentCohort => currentCohort%shorter + enddo !end cohort loop + currentPatch => currentPatch%younger + enddo !end patch loop + + total_stock = biomass_stock + seed_stock + litter_stock + + return + end subroutine SiteCarbonStock + + + end module ChecksBalancesMod diff --git a/main/EDInitMod.F90 b/main/EDInitMod.F90 index 73bda67f5d..085aee9f10 100755 --- a/main/EDInitMod.F90 +++ b/main/EDInitMod.F90 @@ -15,14 +15,17 @@ module EDInitMod use EDGrowthFunctionsMod , only : bdead, bleaf, dbh use EDCohortDynamicsMod , only : create_cohort, fuse_cohorts, sort_cohorts use EDPatchDynamicsMod , only : create_patch - use EDTypesMod , only : ed_site_type, ed_patch_type, ed_cohort_type, area + use EDTypesMod , only : ed_site_type, ed_patch_type, ed_cohort_type use EDTypesMod , only : ncwd use EDTypesMod , only : nuMWaterMem use EDTypesMod , only : maxpft + use EDTypesMod , only : AREA use FatesInterfaceMod , only : bc_in_type use FatesInterfaceMod , only : hlm_use_planthydro use FatesInterfaceMod , only : hlm_use_inventory_init use FatesInterfaceMod , only : numpft + use ChecksBalancesMod , only : SiteCarbonStock + use FatesInterfaceMod , only : nlevsclass ! CIME GLOBALS use shr_log_mod , only : errMsg => shr_log_errMsg @@ -36,6 +39,7 @@ module EDInitMod __FILE__ public :: zero_site + public :: init_site_vars public :: init_patches public :: set_site_properties private :: init_cohorts @@ -46,6 +50,24 @@ module EDInitMod ! ============================================================================ + subroutine init_site_vars( site_in ) + ! + ! !DESCRIPTION: + ! + ! + ! !ARGUMENTS + type(ed_site_type), intent(inout) :: site_in + ! + ! !LOCAL VARIABLES: + !---------------------------------------------------------------------- + ! + allocate(site_in%terminated_nindivs(1:nlevsclass,1:numpft,2)) + allocate(site_in%demotion_rate(1:nlevsclass)) + allocate(site_in%promotion_rate(1:nlevsclass)) + ! + end subroutine init_site_vars + + ! ============================================================================ subroutine zero_site( site_in ) ! ! !DESCRIPTION: @@ -63,9 +85,7 @@ subroutine zero_site( site_in ) site_in%youngest_patch => null() ! pointer to yngest patch at the site ! DISTURBANCE - site_in%disturbance_rate = 0._r8 ! site level disturbance rates from mortality and fire. - site_in%dist_type = 0 ! disturbance dist_type id. - site_in%total_burn_flux_to_atm = 0._r8 ! + site_in%total_burn_flux_to_atm = 0._r8 ! PHENOLOGY site_in%status = 0 ! are leaves in this pixel on or off? @@ -107,6 +127,12 @@ subroutine zero_site( site_in ) site_in%CWD_BG_diagnostic_input_carbonflux(:) = 0._r8 site_in%leaf_litter_diagnostic_input_carbonflux(:) = 0._r8 site_in%root_litter_diagnostic_input_carbonflux(:) = 0._r8 + + ! Resources management (logging/harvesting, etc) + site_in%resources_management%trunk_product_site = 0.0_r8 + + ! canopy spread + site_in%spread = 0._r8 end subroutine zero_site @@ -184,7 +210,7 @@ subroutine set_site_properties( nsites, sites) sites(s)%frac_burnt = 0.0_r8 sites(s)%old_stock = 0.0_r8 - + sites(s)%spread = 1.0_r8 end do return @@ -201,7 +227,6 @@ subroutine init_patches( nsites, sites, bc_in) ! - use EDParamsMod , only : ED_val_maxspread use FatesPlantHydraulicsMod, only : updateSizeDepRhizHydProps use FatesInventoryInitMod, only : initialize_sites_by_inventory @@ -215,10 +240,15 @@ subroutine init_patches( nsites, sites, bc_in) integer :: s real(r8) :: cwd_ag_local(ncwd) real(r8) :: cwd_bg_local(ncwd) - real(r8) :: spread_local(nclmax) real(r8) :: leaf_litter_local(maxpft) real(r8) :: root_litter_local(maxpft) real(r8) :: age !notional age of this patch + + ! dummy locals + real(r8) :: biomass_stock + real(r8) :: litter_stock + real(r8) :: seed_stock + type(ed_patch_type), pointer :: newp ! List out some nominal patch values that are used for Near Bear Ground initializations @@ -228,7 +258,6 @@ subroutine init_patches( nsites, sites, bc_in) cwd_bg_local(:) = 0.0_r8 !ED_val_init_litter leaf_litter_local(:) = 0.0_r8 root_litter_local(:) = 0.0_r8 - spread_local(:) = ED_val_maxspread age = 0.0_r8 ! --------------------------------------------------------------------------------------------- @@ -244,8 +273,12 @@ subroutine init_patches( nsites, sites, bc_in) if (hlm_use_planthydro.eq.itrue) then call updateSizeDepRhizHydProps(sites(s), bc_in(s)) end if + ! For carbon balance checks, we need to initialize the + ! total carbon stock + call SiteCarbonStock(sites(s),sites(s)%old_stock,biomass_stock,litter_stock,seed_stock) + enddo - + else !FIX(SPM,032414) clean this up...inits out of this loop @@ -263,7 +296,7 @@ subroutine init_patches( nsites, sites, bc_in) ! make new patch... call create_patch(sites(s), newp, age, AREA, & - spread_local, cwd_ag_local, cwd_bg_local, leaf_litter_local, & + cwd_ag_local, cwd_bg_local, leaf_litter_local, & root_litter_local) call init_cohorts(newp, bc_in(s)) @@ -275,6 +308,11 @@ subroutine init_patches( nsites, sites, bc_in) call updateSizeDepRhizHydProps(sites(s), bc_in(s)) end if + + ! For carbon balance checks, we need to initialize the + ! total carbon stock + call SiteCarbonStock(sites(s),sites(s)%old_stock,biomass_stock,litter_stock,seed_stock) + enddo end if @@ -362,4 +400,7 @@ subroutine init_cohorts( patch_in, bc_in) end subroutine init_cohorts + ! =============================================================================================== + + end module EDInitMod diff --git a/main/EDMainMod.F90 b/main/EDMainMod.F90 index a9141461e4..6fa2970d0d 100755 --- a/main/EDMainMod.F90 +++ b/main/EDMainMod.F90 @@ -14,6 +14,7 @@ module EDMainMod use FatesInterfaceMod , only : hlm_current_month use FatesInterfaceMod , only : hlm_current_day use FatesInterfaceMod , only : hlm_use_planthydro + use FatesInterfaceMod , only : hlm_reference_date use FatesInterfaceMod , only : hlm_use_ed_st3 use FatesInterfaceMod , only : bc_in_type use FatesInterfaceMod , only : hlm_masterproc @@ -33,12 +34,13 @@ module EDMainMod use EDPhysiologyMod , only : recruitment use EDPhysiologyMod , only : trim_canopy use SFMainMod , only : fire_model - use EDTypesMod , only : get_age_class_index + use FatesSizeAgeTypeIndicesMod, only : get_age_class_index use EDtypesMod , only : ncwd use EDtypesMod , only : ed_site_type use EDtypesMod , only : ed_patch_type use EDtypesMod , only : ed_cohort_type use EDTypesMod , only : do_ed_phenology + use EDTypesMod , only : AREA use FatesConstantsMod , only : itrue,ifalse use FatesPlantHydraulicsMod , only : do_growthrecruiteffects use FatesPlantHydraulicsMod , only : updateSizeDepTreeHydProps @@ -46,8 +48,13 @@ module EDMainMod use FatesPlantHydraulicsMod , only : initTreeHydStates use FatesPlantHydraulicsMod , only : updateSizeDepRhizHydProps ! use FatesPlantHydraulicsMod , only : updateSizeDepRhizHydStates - + use EDLoggingMortalityMod , only : IsItLoggingTime + use FatesGlobals , only : endrun => fates_endrun + use ChecksBalancesMod , only : SiteCarbonStock + ! CIME Globals + use shr_log_mod , only : errMsg => shr_log_errMsg + use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=) implicit none private @@ -64,6 +71,10 @@ module EDMainMod private :: bypass_dynamics logical :: DEBUG = .false. + + character(len=*), parameter, private :: sourcefile = & + __FILE__ + ! ! 10/30/09: Created by Rosie Fisher !----------------------------------------------------------------------- @@ -87,10 +98,14 @@ subroutine ed_ecosystem_dynamics(currentSite, bc_in) if ( hlm_masterproc==itrue ) write(fates_log(),'(A,I4,A,I2.2,A,I2.2)') 'FATES Dynamics: ',& hlm_current_year,'-',hlm_current_month,'-',hlm_current_day + ! Call a routine that simply identifies if logging should occur + ! This is limited to a global event until more structured event handling is enabled + call IsItLoggingTime(hlm_masterproc,currentSite) + !************************************************************************** ! Fire, growth, biogeochemistry. !************************************************************************** - + !FIX(SPM,032414) take this out. On startup these values are all zero and on restart it !zeros out values read in the restart file @@ -104,6 +119,7 @@ subroutine ed_ecosystem_dynamics(currentSite, bc_in) call fire_model(currentSite, bc_in) ! Calculate disturbance and mortality based on previous timestep vegetation. + ! disturbance_rates calls logging mortality and other mortalities, Yi Xu call disturbance_rates(currentSite) end if @@ -130,7 +146,7 @@ subroutine ed_ecosystem_dynamics(currentSite, bc_in) do while (associated(currentPatch)) ! adds small cohort of each PFT - call recruitment(0, currentSite, currentPatch, bc_in) + call recruitment(currentSite, currentPatch, bc_in) currentPatch => currentPatch%younger enddo @@ -209,7 +225,6 @@ subroutine ed_integrate_state_variables(currentSite, bc_in ) ! FIX(SPM,032414) refactor so everything goes through interface ! ! !USES: - use EDTypesMod, only : ageclass_ed ! ! !ARGUMENTS: type(ed_site_type) , intent(inout) :: currentSite @@ -365,7 +380,6 @@ subroutine ed_integrate_state_variables(currentSite, bc_in ) endif enddo - end subroutine ed_integrate_state_variables !-------------------------------------------------------------------------------! @@ -458,6 +472,7 @@ subroutine ed_total_balance_check (currentSite, call_index ) real(r8) :: total_stock ! total ED carbon in KgC/site real(r8) :: change_in_stock ! Change since last time we set ed_allsites_inst%old_stock in this routine. KgC/site real(r8) :: error ! How much carbon did we gain or lose (should be zero!) + real(r8) :: error_frac ! Error as a fraction of total biomass real(r8) :: net_flux ! Difference between recorded fluxes in and out. KgC/site ! nb. There is no time associated with these variables @@ -471,48 +486,78 @@ subroutine ed_total_balance_check (currentSite, call_index ) !----------------------------------------------------------------------- change_in_stock = 0.0_r8 - biomass_stock = 0.0_r8 - litter_stock = 0.0_r8 - - seed_stock = sum(currentSite%seed_bank) - - currentPatch => currentSite%oldest_patch - do while(associated(currentPatch)) - litter_stock = litter_stock + currentPatch%area * (sum(currentPatch%cwd_ag)+ & - sum(currentPatch%cwd_bg)+sum(currentPatch%leaf_litter)+sum(currentPatch%root_litter)) - currentCohort => currentPatch%tallest; - - do while(associated(currentCohort)) - - biomass_stock = biomass_stock + (currentCohort%bdead + currentCohort%balive + & - currentCohort%bstore) * currentCohort%n - currentCohort => currentCohort%shorter; - - enddo !end cohort loop - currentPatch => currentPatch%younger + call SiteCarbonStock(currentSite,total_stock,biomass_stock,litter_stock,seed_stock) - enddo !end patch loop - total_stock = biomass_stock + seed_stock +litter_stock change_in_stock = total_stock - currentSite%old_stock + net_flux = currentSite%flux_in - currentSite%flux_out error = abs(net_flux - change_in_stock) - ! We are not closing this error within 10e-6 very often - ! but this is filling up the logs too much - ! Encapsulating print statements and making new issue (RGK 09-11-2017) - - if ( abs(error) > 10e-6 .and. DEBUG ) then - write(fates_log(),*) 'total error: call index: ',call_index, & - 'in: ',currentSite%flux_in, & - 'out: ',currentSite%flux_out, & - 'net: ',net_flux, & - 'dstock: ',change_in_stock, & - 'error=net_flux-dstock:', error - write(fates_log(),*) 'biomass,litter,seeds', biomass_stock,litter_stock,seed_stock + if(change_in_stock>0.0)then + error_frac = error/abs(total_stock) + else + error_frac = 0.0_r8 + end if + + ! ----------------------------------------------------------------------------------- + ! Terms: + ! %flux_in: accumulates npp over all cohorts, + ! currentSite%flux_in = currentSite%flux_in + & + ! currentCohort%npp_acc * currentCohort%n + ! %flux_out: coarse woody debris going into fragmentation pools: + ! currentSite%flux_out + sum(currentPatch%leaf_litter_out) * & + ! currentPatch%area *hlm_freq_day!kgC/site/day + ! burn fractions: + ! currentSite%flux_out = currentSite%flux_out + & + ! burned_litter * new_patch%area !kG/site/day + ! ----------------------------------------------------------------------------------- + + if ( error_frac > 10e-6 ) then + write(fates_log(),*) 'carbon balance error detected' + write(fates_log(),*) 'error fraction relative to biomass stock:',error_frac + write(fates_log(),*) 'call index: ',call_index + write(fates_log(),*) 'flux in (npp): ',currentSite%flux_in + write(fates_log(),*) 'flux out (fragmentation/harvest): ',currentSite%flux_out + write(fates_log(),*) 'net: ',net_flux + write(fates_log(),*) 'dstock: ',change_in_stock + write(fates_log(),*) 'error=net_flux-dstock:', error + write(fates_log(),*) 'biomass', biomass_stock + write(fates_log(),*) 'litter',litter_stock + write(fates_log(),*) 'seeds',seed_stock + write(fates_log(),*) 'previous total',currentSite%old_stock + + if(DEBUG)then + change_in_stock = 0.0_r8 + biomass_stock = 0.0_r8 + litter_stock = 0.0_r8 + + seed_stock = sum(currentSite%seed_bank)*AREA + currentPatch => currentSite%oldest_patch + do while(associated(currentPatch)) + write(fates_log(),*) '---------------------------------------' + write(fates_log(),*) currentPatch%area , sum(currentPatch%cwd_ag), sum(currentPatch%cwd_bg) + write(fates_log(),*) sum(currentPatch%leaf_litter),sum(currentPatch%root_litter) + write(fates_log(),*)'---' + currentCohort => currentPatch%tallest + do while(associated(currentCohort)) + write(fates_log(),*) currentCohort%bdead,currentCohort%balive,currentCohort%bstore,currentCohort%n + currentCohort => currentCohort%shorter; + enddo !end cohort loop + currentPatch => currentPatch%younger + enddo !end patch loop + end if + write(fates_log(),*) 'lat lon',currentSite%lat,currentSite%lon + + ! If this is the first day of simulation, carbon balance reports but does not end the run + if( int(hlm_current_year*10000 + hlm_current_month*100 + hlm_current_day).ne.hlm_reference_date ) then + write(fates_log(),*) 'aborting on date:',hlm_current_year,hlm_current_month,hlm_current_day + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + endif currentSite%flux_in = 0.0_r8 diff --git a/main/EDParamsMod.F90 b/main/EDParamsMod.F90 index 27a54564c0..b6da0529e7 100644 --- a/main/EDParamsMod.F90 +++ b/main/EDParamsMod.F90 @@ -23,11 +23,8 @@ module EDParamsMod real(r8),protected :: ED_size_diagnostic_scale ! Flag to switch between a linear and exponential ! scale on the plant size axis in diagnostics (NOT USED YET) real(r8),protected :: fates_mortality_disturbance_fraction ! the fraction of canopy mortality that results in disturbance - real(r8),protected :: ED_val_grass_spread real(r8),protected :: ED_val_comp_excln real(r8),protected :: ED_val_stress_mort - real(r8),protected :: ED_val_maxspread ! maximum ratio of dbh to canopy area (cm/m2) - real(r8),protected :: ED_val_minspread ! minimum ratio of dbh to canopy area (cm/m2) real(r8),protected :: ED_val_init_litter real(r8),protected :: ED_val_nignitions real(r8),protected :: ED_val_understorey_death @@ -47,14 +44,16 @@ module EDParamsMod real(r8),protected :: ED_val_phen_coldtemp real(r8),protected :: ED_val_cohort_fusion_tol real(r8),protected :: ED_val_patch_fusion_tol + real(r8),protected :: ED_val_canopy_closure_thresh ! site-level canopy closure point where trees take on forest (narrow) versus savannah (wide) crown allometry + + ! two special parameters whose size is defined in the parameter file + real(r8),protected,allocatable :: ED_val_history_sizeclass_bin_edges(:) + real(r8),protected,allocatable :: ED_val_history_ageclass_bin_edges(:) character(len=param_string_length),parameter :: ED_name_size_diagnostic_scale = "fates_size_diagnostic_scale" character(len=param_string_length),parameter :: ED_name_mort_disturb_frac = "fates_mort_disturb_frac" - character(len=param_string_length),parameter :: ED_name_grass_spread = "fates_grass_spread" character(len=param_string_length),parameter :: ED_name_comp_excln = "fates_comp_excln" character(len=param_string_length),parameter :: ED_name_stress_mort = "fates_stress_mort" - character(len=param_string_length),parameter :: ED_name_maxspread = "fates_maxspread" - character(len=param_string_length),parameter :: ED_name_minspread = "fates_minspread" character(len=param_string_length),parameter :: ED_name_init_litter = "fates_init_litter" character(len=param_string_length),parameter :: ED_name_nignitions = "fates_nignitions" character(len=param_string_length),parameter :: ED_name_understorey_death = "fates_understorey_death" @@ -73,7 +72,12 @@ module EDParamsMod character(len=param_string_length),parameter :: ED_name_phen_ncolddayslim= "fates_phen_ncolddayslim" character(len=param_string_length),parameter :: ED_name_phen_coldtemp= "fates_phen_coldtemp" character(len=param_string_length),parameter :: ED_name_cohort_fusion_tol= "fates_cohort_fusion_tol" - character(len=param_string_length),parameter :: ED_name_patch_fusion_tol= "fates_patch_fusion_tol" + character(len=param_string_length),parameter :: ED_name_patch_fusion_tol= "fates_patch_fusion_tol" + character(len=param_string_length),parameter :: ED_name_canopy_closure_thresh= "fates_canopy_closure_thresh" + + ! non-scalar parameter names + character(len=param_string_length),parameter :: ED_name_history_sizeclass_bin_edges= "fates_history_sizeclass_bin_edges" + character(len=param_string_length),parameter :: ED_name_history_ageclass_bin_edges= "fates_history_ageclass_bin_edges" ! Hydraulics Control Parameters (ONLY RELEVANT WHEN USE_FATES_HYDR = TRUE) ! ---------------------------------------------------------------------------------------------- @@ -92,6 +96,10 @@ module EDParamsMod real(r8),protected :: logging_collateral_frac ! Ratio of collateral mortality to direct logging mortality character(len=param_string_length),parameter :: logging_name_collateral_frac = "fates_logging_collateral_frac" + + real(r8),protected :: logging_coll_under_frac ! Fraction of understory plants that die when logging disturbance + ! is generated + character(len=param_string_length),parameter :: logging_name_coll_under_frac = "fates_logging_coll_under_frac" real(r8),protected :: logging_direct_frac ! Fraction of stems logged per event character(len=param_string_length),parameter :: logging_name_direct_frac = "fates_logging_direct_frac" @@ -121,11 +129,8 @@ subroutine FatesParamsInit() ED_size_diagnostic_scale = nan fates_mortality_disturbance_fraction = nan - ED_val_grass_spread = nan ED_val_comp_excln = nan ED_val_stress_mort = nan - ED_val_maxspread = nan - ED_val_minspread = nan ED_val_init_litter = nan ED_val_nignitions = nan ED_val_understorey_death = nan @@ -145,6 +150,7 @@ subroutine FatesParamsInit() ED_val_phen_coldtemp = nan ED_val_cohort_fusion_tol = nan ED_val_patch_fusion_tol = nan + ED_val_canopy_closure_thresh = nan hydr_psi0 = nan hydr_psicap = nan @@ -164,12 +170,15 @@ subroutine FatesRegisterParams(fates_params) ! that need to be synced with host values. use FatesParametersInterface, only : fates_parameters_type, dimension_name_scalar1d, dimension_shape_1d + use FatesParametersInterface, only : dimension_name_history_size_bins, dimension_name_history_age_bins implicit none class(fates_parameters_type), intent(inout) :: fates_params character(len=param_string_length), parameter :: dim_names(1) = (/dimension_name_scalar1d/) + character(len=param_string_length), parameter :: dim_names_sizeclass(1) = (/dimension_name_history_size_bins/) + character(len=param_string_length), parameter :: dim_names_ageclass(1) = (/dimension_name_history_age_bins/) call FatesParamsInit() @@ -179,21 +188,12 @@ subroutine FatesRegisterParams(fates_params) call fates_params%RegisterParameter(name=ED_name_mort_disturb_frac, dimension_shape=dimension_shape_1d, & dimension_names=dim_names) - call fates_params%RegisterParameter(name=ED_name_grass_spread, dimension_shape=dimension_shape_1d, & - dimension_names=dim_names) - call fates_params%RegisterParameter(name=ED_name_comp_excln, dimension_shape=dimension_shape_1d, & dimension_names=dim_names) call fates_params%RegisterParameter(name=ED_name_stress_mort, dimension_shape=dimension_shape_1d, & dimension_names=dim_names) - call fates_params%RegisterParameter(name=ED_name_maxspread, dimension_shape=dimension_shape_1d, & - dimension_names=dim_names) - - call fates_params%RegisterParameter(name=ED_name_minspread, dimension_shape=dimension_shape_1d, & - dimension_names=dim_names) - call fates_params%RegisterParameter(name=ED_name_init_litter, dimension_shape=dimension_shape_1d, & dimension_names=dim_names) @@ -251,6 +251,9 @@ subroutine FatesRegisterParams(fates_params) call fates_params%RegisterParameter(name=ED_name_patch_fusion_tol, dimension_shape=dimension_shape_1d, & dimension_names=dim_names) + call fates_params%RegisterParameter(name=ED_name_canopy_closure_thresh, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names) + call fates_params%RegisterParameter(name=hydr_name_psi0, dimension_shape=dimension_shape_1d, & dimension_names=dim_names) @@ -263,6 +266,9 @@ subroutine FatesRegisterParams(fates_params) call fates_params%RegisterParameter(name=logging_name_collateral_frac, dimension_shape=dimension_shape_1d, & dimension_names=dim_names) + call fates_params%RegisterParameter(name=logging_name_coll_under_frac, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names) + call fates_params%RegisterParameter(name=logging_name_direct_frac, dimension_shape=dimension_shape_1d, & dimension_names=dim_names) @@ -273,6 +279,13 @@ subroutine FatesRegisterParams(fates_params) dimension_names=dim_names) + ! non-scalar parameters + call fates_params%RegisterParameter(name=ED_name_history_sizeclass_bin_edges, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names_sizeclass) + + call fates_params%RegisterParameter(name=ED_name_history_ageclass_bin_edges, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names_ageclass) + end subroutine FatesRegisterParams @@ -291,21 +304,12 @@ subroutine FatesReceiveParams(fates_params) call fates_params%RetreiveParameter(name=ED_name_mort_disturb_frac, & data=fates_mortality_disturbance_fraction) - call fates_params%RetreiveParameter(name=ED_name_grass_spread, & - data=ED_val_grass_spread) - call fates_params%RetreiveParameter(name=ED_name_comp_excln, & data=ED_val_comp_excln) call fates_params%RetreiveParameter(name=ED_name_stress_mort, & data=ED_val_stress_mort) - call fates_params%RetreiveParameter(name=ED_name_maxspread, & - data=ED_val_maxspread) - - call fates_params%RetreiveParameter(name=ED_name_minspread, & - data=ED_val_minspread) - call fates_params%RetreiveParameter(name=ED_name_init_litter, & data=ED_val_init_litter) @@ -363,6 +367,9 @@ subroutine FatesReceiveParams(fates_params) call fates_params%RetreiveParameter(name=ED_name_patch_fusion_tol, & data=ED_val_patch_fusion_tol) + call fates_params%RetreiveParameter(name=ED_name_canopy_closure_thresh, & + data=ED_val_canopy_closure_thresh) + call fates_params%RetreiveParameter(name=hydr_name_psi0, & data=hydr_psi0) @@ -374,6 +381,9 @@ subroutine FatesReceiveParams(fates_params) call fates_params%RetreiveParameter(name=logging_name_collateral_frac, & data=logging_collateral_frac) + + call fates_params%RetreiveParameter(name=logging_name_coll_under_frac, & + data=logging_coll_under_frac) call fates_params%RetreiveParameter(name=logging_name_direct_frac, & data=logging_direct_frac) @@ -384,6 +394,14 @@ subroutine FatesReceiveParams(fates_params) call fates_params%RetreiveParameter(name=logging_name_event_code, & data=logging_event_code) + ! parameters that are arrays of size defined within the params file and thus need allocating as well + call fates_params%RetreiveParameterAllocate(name=ED_name_history_sizeclass_bin_edges, & + data=ED_val_history_sizeclass_bin_edges) + + call fates_params%RetreiveParameterAllocate(name=ED_name_history_ageclass_bin_edges, & + data=ED_val_history_ageclass_bin_edges) + + end subroutine FatesReceiveParams ! ===================================================================================== @@ -400,11 +418,8 @@ subroutine FatesReportParams(is_master) write(fates_log(),*) '----------- FATES Scalar Parameters -----------------' write(fates_log(),fmt0) 'ED_size_diagnostic_scale = ',ED_size_diagnostic_scale write(fates_log(),fmt0) 'fates_mortality_disturbance_fraction = ',fates_mortality_disturbance_fraction - write(fates_log(),fmt0) 'ED_val_grass_spread = ',ED_val_grass_spread - write(fates_log(),fmt0) 'ED_val_comp_excln = ', ED_val_comp_excln + write(fates_log(),fmt0) 'ED_val_comp_excln = ',ED_val_comp_excln write(fates_log(),fmt0) 'ED_val_stress_mort = ',ED_val_stress_mort - write(fates_log(),fmt0) 'ED_val_maxspread = ',ED_val_maxspread - write(fates_log(),fmt0) 'ED_val_minspread = ',ED_val_minspread write(fates_log(),fmt0) 'ED_val_init_litter = ',ED_val_init_litter write(fates_log(),fmt0) 'ED_val_nignitions = ',ED_val_nignitions write(fates_log(),fmt0) 'ED_val_understorey_death = ',ED_val_understorey_death @@ -424,6 +439,7 @@ subroutine FatesReportParams(is_master) write(fates_log(),fmt0) 'ED_val_phen_coldtemp = ',ED_val_phen_coldtemp write(fates_log(),fmt0) 'ED_val_cohort_fusion_tol = ',ED_val_cohort_fusion_tol write(fates_log(),fmt0) 'ED_val_patch_fusion_tol = ',ED_val_patch_fusion_tol + write(fates_log(),fmt0) 'ED_val_canopy_closure_thresh = ',ED_val_canopy_closure_thresh write(fates_log(),fmt0) 'hydr_psi0 = ',hydr_psi0 write(fates_log(),fmt0) 'hydr_psicap = ',hydr_psicap write(fates_log(),fmt0) 'logging_dbhmin = ',logging_dbhmin diff --git a/main/EDPftvarcon.F90 b/main/EDPftvarcon.F90 index d4b9fdf18a..465a28fc59 100644 --- a/main/EDPftvarcon.F90 +++ b/main/EDPftvarcon.F90 @@ -26,10 +26,11 @@ module EDPftvarcon !ED specific variables. type, public :: EDPftvarcon_type real(r8), allocatable :: pft_used (:) ! Switch to turn on and off PFTs - real(r8), allocatable :: max_dbh (:) ! maximum dbh at which height growth ceases... + real(r8), allocatable :: freezetol (:) ! minimum temperature tolerance (NOT CURRENTY USED) real(r8), allocatable :: wood_density (:) ! wood density g cm^-3 ... real(r8), allocatable :: hgt_min (:) ! sapling height m + real(r8), allocatable :: dbh_repro_threshold(:) ! diameter at which mature plants shift allocation real(r8), allocatable :: dleaf (:) ! leaf characteristic dimension length (m) real(r8), allocatable :: z0mr (:) ! ratio of roughness length of vegetation to height (-) real(r8), allocatable :: displar (:) ! ratio of displacement height to canopy top height (-) @@ -45,6 +46,7 @@ module EDPftvarcon real(r8), allocatable :: root_long (:) ! root longevity (yrs) real(r8), allocatable :: clone_alloc (:) ! fraction of carbon balance allocated to clonal reproduction. real(r8), allocatable :: seed_alloc (:) ! fraction of carbon balance allocated to seeds. + real(r8), allocatable :: c2b (:) ! Carbon to biomass multiplier [kg/kgC] real(r8), allocatable :: woody(:) real(r8), allocatable :: stress_decid(:) real(r8), allocatable :: season_decid(:) @@ -98,7 +100,9 @@ module EDPftvarcon ! Equation 16 Thonicke et al 2010 ! Allometry Parameters - ! -------------------------------------------------------------------------------------------- + ! -------------------------------------------------------------------------------------------- + real(r8), allocatable :: allom_dbh_maxheight(:) ! dbh at which height growth ceases + real(r8), allocatable :: allom_hmode(:) ! height allometry function type real(r8), allocatable :: allom_lmode(:) ! maximum leaf allometry function type real(r8), allocatable :: allom_fmode(:) ! maximum root allometry function type @@ -117,13 +121,22 @@ module EDPftvarcon real(r8), allocatable :: allom_d2bl2(:) ! Parameter 2 for d2bl allometry (slope) real(r8), allocatable :: allom_d2bl3(:) ! Parameter 3 for d2bl allometry (optional) real(r8), allocatable :: allom_sai_scaler(:) ! - real(r8), allocatable :: allom_d2bl_slascaler(:) ! real(r8), allocatable :: allom_blca_expnt_diff(:) ! Any difference in the exponent between the leaf ! biomass and crown area scaling + real(r8), allocatable :: allom_d2ca_coefficient_max(:) ! upper (savanna) value for crown area to dbh coefficient + real(r8), allocatable :: allom_d2ca_coefficient_min(:) ! lower (closed-canopy forest) value for crown area to dbh coefficient real(r8), allocatable :: allom_agb1(:) ! Parameter 1 for agb allometry real(r8), allocatable :: allom_agb2(:) ! Parameter 2 for agb allometry real(r8), allocatable :: allom_agb3(:) ! Parameter 3 for agb allometry real(r8), allocatable :: allom_agb4(:) ! Parameter 3 for agb allometry + + ! Prescribed Physiology Mode Parameters + real(r8), allocatable :: prescribed_npp_canopy(:) ! this is only for the special prescribed_physiology_mode + real(r8), allocatable :: prescribed_npp_understory(:) ! this is only for the special prescribed_physiology_mode + real(r8), allocatable :: prescribed_mortality_canopy(:) ! this is only for the special prescribed_physiology_mode + real(r8), allocatable :: prescribed_mortality_understory(:) ! this is only for the special prescribed_physiology_mode + real(r8), allocatable :: prescribed_recruitment(:) ! this is only for the special prescribed_physiology_mode + ! Plant Hydraulic Parameters ! --------------------------------------------------------------------------------------------- @@ -242,7 +255,7 @@ subroutine Register_PFT(this, fates_params) call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) - name = 'fates_max_dbh' + name = 'fates_dbh_repro_threshold' call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) @@ -302,6 +315,10 @@ subroutine Register_PFT(this, fates_params) call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) + name = 'fates_c2b' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) + name = 'fates_woody' call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) @@ -394,10 +411,34 @@ subroutine Register_PFT(this, fates_params) call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) + name = 'fates_prescribed_npp_canopy' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) + + name = 'fates_prescribed_npp_understory' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) + + name = 'fates_prescribed_mortality_canopy' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) + + name = 'fates_prescribed_mortality_understory' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) + + name = 'fates_prescribed_recruitment' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) + name = 'fates_alpha_SH' call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) + name = 'fates_allom_dbh_maxheight' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) + name = 'fates_allom_hmode' call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) @@ -462,7 +503,11 @@ subroutine Register_PFT(this, fates_params) call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) - name = 'fates_allom_d2bl_slascaler' + name = 'fates_allom_d2ca_coefficient_max' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) + + name = 'fates_allom_d2ca_coefficient_min' call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) @@ -610,9 +655,9 @@ subroutine Receive_PFT(this, fates_params) call fates_params%RetreiveParameterAllocate(name=name, & data=this%pft_used) - name = 'fates_max_dbh' + name = 'fates_dbh_repro_threshold' call fates_params%RetreiveParameterAllocate(name=name, & - data=this%max_dbh) + data=this%dbh_repro_threshold) name = 'fates_freezetol' call fates_params%RetreiveParameterAllocate(name=name, & @@ -670,6 +715,10 @@ subroutine Receive_PFT(this, fates_params) call fates_params%RetreiveParameterAllocate(name=name, & data=this%seed_alloc) + name = 'fates_c2b' + call fates_params%RetreiveParameterAllocate(name=name, & + data=this%c2b) + name = 'fates_woody' call fates_params%RetreiveParameterAllocate(name=name, & data=this%woody) @@ -758,10 +807,34 @@ subroutine Receive_PFT(this, fates_params) call fates_params%RetreiveParameterAllocate(name=name, & data=this%grperc) + name = 'fates_prescribed_npp_canopy' + call fates_params%RetreiveParameterAllocate(name=name, & + data=this%prescribed_npp_canopy) + + name = 'fates_prescribed_npp_understory' + call fates_params%RetreiveParameterAllocate(name=name, & + data=this%prescribed_npp_understory) + + name = 'fates_prescribed_mortality_canopy' + call fates_params%RetreiveParameterAllocate(name=name, & + data=this%prescribed_mortality_canopy) + + name = 'fates_prescribed_mortality_understory' + call fates_params%RetreiveParameterAllocate(name=name, & + data=this%prescribed_mortality_understory) + + name = 'fates_prescribed_recruitment' + call fates_params%RetreiveParameterAllocate(name=name, & + data=this%prescribed_recruitment) + name = 'fates_alpha_SH' call fates_params%RetreiveParameterAllocate(name=name, & data=this%fire_alpha_SH) + name = 'fates_allom_dbh_maxheight' + call fates_params%RetreiveParameterAllocate(name=name, & + data=this%allom_dbh_maxheight) + name = 'fates_allom_hmode' call fates_params%RetreiveParameterAllocate(name=name, & data=this%allom_hmode) @@ -830,9 +903,13 @@ subroutine Receive_PFT(this, fates_params) call fates_params%RetreiveParameterAllocate(name=name, & data=this%allom_blca_expnt_diff) - name = 'fates_allom_d2bl_slascaler' + name = 'fates_allom_d2ca_coefficient_max' call fates_params%RetreiveParameterAllocate(name=name, & - data=this%allom_d2bl_slascaler) + data=this%allom_d2ca_coefficient_max) + + name = 'fates_allom_d2ca_coefficient_min' + call fates_params%RetreiveParameterAllocate(name=name, & + data=this%allom_d2ca_coefficient_min) name = 'fates_allom_sai_scaler' call fates_params%RetreiveParameterAllocate(name=name, & @@ -954,7 +1031,6 @@ subroutine Receive_PFT(this, fates_params) call fates_params%RetreiveParameterAllocate(name=name, & data=this%displar) - end subroutine Receive_PFT !----------------------------------------------------------------------- @@ -1313,7 +1389,8 @@ subroutine FatesReportPFTParams(is_master) write(fates_log(),*) '----------- FATES PFT Parameters -----------------' write(fates_log(),fmt0) 'pft_used = ',EDPftvarcon_inst%pft_used - write(fates_log(),fmt0) 'max_dbh = ',EDPftvarcon_inst%max_dbh + write(fates_log(),fmt0) 'dbh max height = ',EDPftvarcon_inst%allom_dbh_maxheight + write(fates_log(),fmt0) 'dbh mature = ',EDPftvarcon_inst%dbh_repro_threshold write(fates_log(),fmt0) 'freezetol = ',EDPftvarcon_inst%freezetol write(fates_log(),fmt0) 'wood_density = ',EDPftvarcon_inst%wood_density write(fates_log(),fmt0) 'hgt_min = ',EDPftvarcon_inst%hgt_min @@ -1331,6 +1408,7 @@ subroutine FatesReportPFTParams(is_master) write(fates_log(),fmt0) 'root_long = ',EDPftvarcon_inst%root_long write(fates_log(),fmt0) 'clone_alloc = ',EDPftvarcon_inst%clone_alloc write(fates_log(),fmt0) 'seed_alloc = ',EDPftvarcon_inst%seed_alloc + write(fates_log(),fmt0) 'C2B = ',EDPftvarcon_inst%c2b write(fates_log(),fmt0) 'woody = ',EDPftvarcon_inst%woody write(fates_log(),fmt0) 'stress_decid = ',EDPftvarcon_inst%stress_decid write(fates_log(),fmt0) 'season_decid = ',EDPftvarcon_inst%season_decid @@ -1392,8 +1470,9 @@ subroutine FatesReportPFTParams(is_master) write(fates_log(),fmt0) 'allom_d2bl2 = ',EDPftvarcon_inst%allom_d2bl2 write(fates_log(),fmt0) 'allom_d2bl3 = ',EDPftvarcon_inst%allom_d2bl3 write(fates_log(),fmt0) 'allom_sai_scaler = ',EDPftvarcon_inst%allom_sai_scaler - write(fates_log(),fmt0) 'allom_d2bl_slascaler = ',EDPftvarcon_inst%allom_d2bl_slascaler write(fates_log(),fmt0) 'allom_blca_expnt_diff = ',EDPftvarcon_inst%allom_blca_expnt_diff + write(fates_log(),fmt0) 'allom_d2ca_coefficient_max = ',EDPftvarcon_inst%allom_d2ca_coefficient_max + write(fates_log(),fmt0) 'allom_d2ca_coefficient_min = ',EDPftvarcon_inst%allom_d2ca_coefficient_min write(fates_log(),fmt0) 'allom_agb1 = ',EDPftvarcon_inst%allom_agb1 write(fates_log(),fmt0) 'allom_agb2 = ',EDPftvarcon_inst%allom_agb2 write(fates_log(),fmt0) 'allom_agb3 = ',EDPftvarcon_inst%allom_agb3 diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 8d0690a658..82265c093c 100755 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -13,6 +13,7 @@ module EDTypesMod integer, parameter :: maxPatchesPerSite = 10 ! maximum number of patches to live on a site integer, parameter :: maxCohortsPerPatch = 160 ! maximum number of cohorts per patch + integer, parameter :: nclmax = 2 ! Maximum number of canopy layers integer, parameter :: ican_upper = 1 ! Nominal index for the upper canopy integer, parameter :: ican_ustory = 2 ! Nominal index for understory in two-canopy system @@ -65,8 +66,10 @@ module EDTypesMod integer , parameter :: external_recruitment = 0 ! external recruitment flag 1=yes integer , parameter :: SENES = 10 ! Window of time over which we track temp for cold sensecence (days) real(r8), parameter :: DINC_ED = 1.0_r8 ! size of LAI bins. - integer , parameter :: N_DIST_TYPES = 2 ! number of disturbance types (mortality, fire) - + integer , parameter :: N_DIST_TYPES = 3 ! Disturbance Modes 1) tree-fall, 2) fire, 3) logging + integer , parameter :: dtype_ifall = 1 ! index for naturally occuring tree-fall generated event + integer , parameter :: dtype_ifire = 2 ! index for fire generated disturbance event + integer , parameter :: dtype_ilog = 3 ! index for logging generated disturbance event ! SPITFIRE integer, parameter :: NCWD = 4 ! number of coarse woody debris pools (twig,s branch,l branch, trunk) @@ -97,25 +100,6 @@ module EDTypesMod ! special mode to cause PFTs to create seed mass of all currently-existing PFTs logical, parameter :: homogenize_seed_pfts = .false. - !the lower limit of the size classes of ED cohorts - !0-10,10-20... - integer, parameter :: nlevsclass_ed = 13 ! Number of dbh size classes for size structure analysis - ! |0-1,1-2,2-3,3-4,4-5,5-10,10-20,20-30,30-40,40-50,50-60,60-70,70-80,80-90,90-100,100+| -! real(r8), parameter, dimension(16) :: sclass_ed = (/0.0_r8,1.0_r8,2.0_r8,3.0_r8,4.0_r8,5.0_r8,10.0_r8,20.0_r8,30.0_r8,40.0_r8, & -! 50.0_r8,60.0_r8,70.0_r8,80.0_r8,90.0_r8,100.0_r8/) - - real(r8), parameter, dimension(nlevsclass_ed) :: sclass_ed = (/0.0_r8,5.0_r8,10.0_r8,15.0_r8,20.0_r8,30.0_r8,40.0_r8, & - 50.0_r8,60.0_r8,70.0_r8,80.0_r8,90.0_r8,100.0_r8/) - - integer, parameter :: nlevage_ed = 7 ! Number of patch-age classes for age structured analyses - real(r8), parameter, dimension(nlevage_ed) :: ageclass_ed = (/0.0_r8,1.0_r8,2._r8,5.0_r8,10.0_r8,20.0_r8,50.0_r8/) - - - ! integer, parameter :: nlevsclass_ed = 17 - ! real(r8), parameter, dimension(17) :: sclass_ed = (/0.1_r8, 5.0_r8,10.0_r8,15.0_r8,20.0_r8,25.0_r8, & - ! 30.0_r8,35.0_r8,40.0_r8,45.0_r8,50.0_r8,55.0_r8, & - ! 60.0_r8,65.0_r8,70.0_r8,75.0_r8,80.0_r8/) - integer, parameter :: nlevmclass_ed = 5 ! nlev "mortality" classes in ED ! Number of ways to die ! (background,hydraulic,carbon,impact,fire) @@ -124,8 +108,6 @@ module EDTypesMod (/"background","hydraulic ","carbon ","impact ","fire "/) - - !************************************ !** COHORT type structure ** !************************************ @@ -247,6 +229,13 @@ module EDTypesMod real(r8) :: imort ! mortality from impacts by others n/year real(r8) :: fmort ! fire mortality n/year + ! Logging Mortality Rate + ! Yi Xu + real(r8) :: lmort_logging ! directly logging rate %/per logging activity + real(r8) :: lmort_collateral ! collaterally damaged rate %/per logging activity + real(r8) :: lmort_infra ! mechanically damaged rate %/per logging activity + + ! NITROGEN POOLS ! ---------------------------------------------------------------------------------- ! Nitrogen pools are not prognostic in the current implementation. @@ -299,7 +288,6 @@ module EDTypesMod integer :: ncl_p ! Number of occupied canopy layers ! LEAF ORGANIZATION - real(r8) :: spread(nclmax) ! dynamic ratio of dbh to canopy area: cm/m2 real(r8) :: pft_agb_profile(maxpft,n_dbh_bins) ! binned above ground biomass, for patch fusion: KgC/m2 real(r8) :: canopy_layer_lai(nclmax) ! lai that is shading this canopy layer: m2/m2 real(r8) :: total_canopy_area ! area that is covered by vegetation : m2 @@ -368,7 +356,9 @@ module EDTypesMod real(r8) :: btran_ft(maxpft) ! btran calculated seperately for each PFT:- ! DISTURBANCE - real(r8) :: disturbance_rates(n_dist_types) ! disturbance rate from 1) mortality and 2) fire: fraction/day + real(r8) :: disturbance_rates(n_dist_types) ! disturbance rate from 1) mortality + ! 2) fire: fraction/day + ! 3) logging mortatliy real(r8) :: disturbance_rate ! larger effective disturbance rate: fraction/day ! LITTER AND COARSE WOODY DEBRIS @@ -430,6 +420,7 @@ module EDTypesMod real(r8) :: tfc_ros ! total fuel consumed - no trunks. KgC/m2/day real(r8) :: burnt_frac_litter(nfsc) ! fraction of each litter pool burned:- + ! PLANT HYDRAULICS type(ed_patch_hydr_type) , pointer :: pa_hydr ! All patch hydraulics data, see FatesHydraulicsMemMod.F90 @@ -437,15 +428,38 @@ module EDTypesMod end type ed_patch_type + + !************************************ + !** Resources management type ** + ! YX + !************************************ + type ed_resources_management_type + + real(r8) :: trunk_product_site ! Actual trunk product at site level KgC/site + + !debug variables + real(r8) :: delta_litter_stock + real(r8) :: delta_biomass_stock + real(r8) :: delta_individual + + end type ed_resources_management_type + + + !************************************ !** Site type structure ** !************************************ type ed_site_type - + ! POINTERS type (ed_patch_type), pointer :: oldest_patch => null() ! pointer to oldest patch at the site type (ed_patch_type), pointer :: youngest_patch => null() ! pointer to yngest patch at the site + + ! Resource management + type (ed_resources_management_type) :: resources_management ! resources_management at the site + + ! INDICES real(r8) :: lat ! latitude: degrees @@ -488,12 +502,6 @@ module EDTypesMod real(r8) :: nbp_integrated ! Net biosphere production accumulated over model time-steps [gC/m2] - ! DISTURBANCE - real(r8) :: disturbance_mortality ! site level disturbance rates from mortality. - real(r8) :: disturbance_fire ! site level disturbance rates from fire. - integer :: dist_type ! disturbance dist_type id. - real(r8) :: disturbance_rate ! site total dist rate - ! PHENOLOGY real(r8) :: ED_GDD_site ! ED Phenology growing degree days. integer :: status ! are leaves in this pixel on or off for cold decid @@ -525,12 +533,12 @@ module EDTypesMod ! TERMINATION, RECRUITMENT, DEMOTION, and DISTURBANCE - real(r8) :: terminated_nindivs(1:nlevsclass_ed,1:maxpft,2) ! number of individuals that were in cohorts which were terminated this timestep, on size x pft x canopy array. + real(r8), allocatable :: terminated_nindivs(:,:,:) ! number of individuals that were in cohorts which were terminated this timestep, on size x pft x canopy array. real(r8) :: termination_carbonflux(2) ! carbon flux from live to dead pools associated with termination mortality, per canopy level real(r8) :: recruitment_rate(1:maxpft) ! number of individuals that were recruited into new cohorts - real(r8) :: demotion_rate(1:nlevsclass_ed) ! rate of individuals demoted from canopy to understory per FATES timestep + real(r8), allocatable :: demotion_rate(:) ! rate of individuals demoted from canopy to understory per FATES timestep real(r8) :: demotion_carbonflux ! biomass of demoted individuals from canopy to understory [kgC/ha/day] - real(r8) :: promotion_rate(1:nlevsclass_ed) ! rate of individuals promoted from understory to canopy per FATES timestep + real(r8), allocatable :: promotion_rate(:) ! rate of individuals promoted from understory to canopy per FATES timestep real(r8) :: promotion_carbonflux ! biomass of promoted individuals from understory to canopy [kgC/ha/day] ! some diagnostic-only (i.e. not resolved by ODE solver) flux of carbon to CWD and litter pools from termination and canopy mortality @@ -539,71 +547,13 @@ module EDTypesMod real(r8) :: leaf_litter_diagnostic_input_carbonflux(1:maxpft) ! diagnostic flux to AG litter [kg C / m2 / yr] real(r8) :: root_litter_diagnostic_input_carbonflux(1:maxpft) ! diagnostic flux to BG litter [kg C / m2 / yr] + ! Canopy Spread + real(r8) :: spread ! dynamic canopy allometric term [unitless] + end type ed_site_type contains - ! ===================================================================================== - - function get_age_class_index(age) result( patch_age_class ) - - real(r8), intent(in) :: age - - integer :: patch_age_class - - patch_age_class = count(age-ageclass_ed.ge.0.0_r8) - - end function get_age_class_index - - ! ===================================================================================== - - function get_sizeage_class_index(dbh,age) result(size_by_age_class) - - ! Arguments - real(r8),intent(in) :: dbh - real(r8),intent(in) :: age - - integer :: size_class - integer :: age_class - integer :: size_by_age_class - - size_class = get_size_class_index(dbh) - - age_class = get_age_class_index(age) - - size_by_age_class = (age_class-1)*nlevsclass_ed + size_class - - end function get_sizeage_class_index - - ! ===================================================================================== - - subroutine sizetype_class_index(dbh,pft,size_class,size_by_pft_class) - - ! Arguments - real(r8),intent(in) :: dbh - integer,intent(in) :: pft - integer,intent(out) :: size_class - integer,intent(out) :: size_by_pft_class - - size_class = get_size_class_index(dbh) - - size_by_pft_class = (pft-1)*nlevsclass_ed+size_class - - return - end subroutine sizetype_class_index - - ! ===================================================================================== - - function get_size_class_index(dbh) result(cohort_size_class) - - real(r8), intent(in) :: dbh - - integer :: cohort_size_class - - cohort_size_class = count(dbh-sclass_ed.ge.0.0_r8) - - end function get_size_class_index - ! ===================================================================================== subroutine val_check_ed_vars(currentPatch,var_aliases,return_code) diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index ca91ecd640..fa7d98d539 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -14,6 +14,7 @@ module FatesHistoryInterfaceMod use FatesInterfaceMod , only : hlm_use_ed_st3 use FatesInterfaceMod , only : numpft use EDParamsMod , only : ED_val_comp_excln + use FatesInterfaceMod , only : nlevsclass, nlevage ! FIXME(bja, 2016-10) need to remove CLM dependancy use EDPftvarcon , only : EDPftvarcon_inst @@ -41,7 +42,6 @@ module FatesHistoryInterfaceMod integer, private :: ih_trimming_pa integer, private :: ih_area_plant_pa integer, private :: ih_area_treespread_pa - integer, private :: ih_canopy_spread_pa integer, private :: ih_nesterov_fire_danger_pa integer, private :: ih_spitfire_ROS_pa integer, private :: ih_effect_wspeed_pa @@ -121,6 +121,7 @@ module FatesHistoryInterfaceMod integer, private :: ih_promotion_carbonflux_si integer, private :: ih_canopy_mortality_carbonflux_si integer, private :: ih_understory_mortality_carbonflux_si + integer, private :: ih_canopy_spread_si ! Indices to (site x scpf) variables integer, private :: ih_nplant_si_scpf @@ -159,6 +160,9 @@ module FatesHistoryInterfaceMod integer, private :: ih_m5_si_scpf integer, private :: ih_m6_si_scpf + !LOGGING , make sure to add ih_m7_si_scpf and hio_m7_si_scpf + integer, private :: ih_m7_si_scpf + integer, private :: ih_ar_si_scpf integer, private :: ih_ar_grow_si_scpf integer, private :: ih_ar_maint_si_scpf @@ -246,6 +250,7 @@ module FatesHistoryInterfaceMod integer, private :: ih_ncl_si_age integer, private :: ih_npatches_si_age integer, private :: ih_zstar_si_age + integer, private :: ih_biomass_si_age ! Indices to hydraulics variables @@ -866,20 +871,17 @@ subroutine flush_hvars(this,nc,upfreq_in) class(fates_history_interface_type) :: this integer,intent(in) :: nc integer,intent(in) :: upfreq_in - integer :: ivar - type(fates_history_variable_type),pointer :: hvar integer :: lb1,ub1,lb2,ub2 do ivar=1,ubound(this%hvars,1) - associate( hvar => this%hvars(ivar) ) - if (hvar%upfreq == upfreq_in) then ! Only flush variables with update on dynamics step - call hvar%Flush(nc, this%dim_bounds, this%dim_kinds) - end if - end associate + if (this%hvars(ivar)%upfreq == upfreq_in) then ! Only flush variables with update on dynamics step + call this%hvars(ivar)%flush(nc, this%dim_bounds, this%dim_kinds) + + end if end do - end subroutine flush_hvars +end subroutine flush_hvars ! ===================================================================================== @@ -912,7 +914,6 @@ subroutine set_history_var(this, vname, units, long, use_default, avgflag, vtype ! not used ! locals - type(fates_history_variable_type), pointer :: hvar integer :: ub1, lb1, ub2, lb2 ! Bounds for allocating the var integer :: ityp @@ -1097,17 +1098,15 @@ subroutine update_history_dyn(this,nc,nsites,sites) use EDtypesMod , only : ed_patch_type use EDtypesMod , only : AREA use EDtypesMod , only : AREA_INV - use EDtypesMod , only : nlevsclass_ed - use EDtypesMod , only : nlevage_ed use EDtypesMod , only : nfsc use EDtypesMod , only : ncwd use EDtypesMod , only : ican_upper use EDtypesMod , only : ican_ustory - use EDTypesMod , only : get_sizeage_class_index + use FatesSizeAgeTypeIndicesMod, only : get_sizeage_class_index use EDTypesMod , only : nlevleaf ! Arguments - class(fates_history_interface_type) :: this + class(fates_history_interface_type) :: this integer , intent(in) :: nc ! clump index integer , intent(in) :: nsites type(ed_site_type) , intent(inout), target :: sites(nsites) @@ -1132,7 +1131,6 @@ subroutine update_history_dyn(this,nc,nsites,sites) real(r8) :: patch_scaling_scalar ! ratio of canopy to patch area for counteracting patch scaling real(r8) :: dbh ! diameter ("at breast height") - type(fates_history_variable_type),pointer :: hvar type(ed_patch_type),pointer :: cpatch type(ed_cohort_type),pointer :: ccohort @@ -1143,7 +1141,7 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_trimming_pa => this%hvars(ih_trimming_pa)%r81d, & hio_area_plant_pa => this%hvars(ih_area_plant_pa)%r81d, & hio_area_treespread_pa => this%hvars(ih_area_treespread_pa)%r81d, & - hio_canopy_spread_pa => this%hvars(ih_canopy_spread_pa)%r81d, & + hio_canopy_spread_si => this%hvars(ih_canopy_spread_si)%r81d, & hio_biomass_si_pft => this%hvars(ih_biomass_si_pft)%r82d, & hio_leafbiomass_si_pft => this%hvars(ih_leafbiomass_si_pft)%r82d, & hio_storebiomass_si_pft => this%hvars(ih_storebiomass_si_pft)%r82d, & @@ -1210,6 +1208,9 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_m4_si_scpf => this%hvars(ih_m4_si_scpf)%r82d, & hio_m5_si_scpf => this%hvars(ih_m5_si_scpf)%r82d, & hio_m6_si_scpf => this%hvars(ih_m6_si_scpf)%r82d, & + + hio_m7_si_scpf => this%hvars(ih_m7_si_scpf)%r82d, & + hio_ba_si_scls => this%hvars(ih_ba_si_scls)%r82d, & hio_nplant_canopy_si_scls => this%hvars(ih_nplant_canopy_si_scls)%r82d, & hio_nplant_understory_si_scls => this%hvars(ih_nplant_understory_si_scls)%r82d, & @@ -1261,6 +1262,7 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_ncl_si_age => this%hvars(ih_ncl_si_age)%r82d, & hio_npatches_si_age => this%hvars(ih_npatches_si_age)%r82d, & hio_zstar_si_age => this%hvars(ih_zstar_si_age)%r82d, & + hio_biomass_si_age => this%hvars(ih_biomass_si_age)%r82d, & hio_litter_moisture_si_fuel => this%hvars(ih_litter_moisture_si_fuel)%r82d, & hio_cwd_ag_si_cwdsc => this%hvars(ih_cwd_ag_si_cwdsc)%r82d, & hio_cwd_bg_si_cwdsc => this%hvars(ih_cwd_bg_si_cwdsc)%r82d, & @@ -1305,6 +1307,8 @@ subroutine update_history_dyn(this,nc,nsites,sites) ! The seed bank is a site level variable hio_seed_bank_si(io_si) = sum(sites(s)%seed_bank) * g_per_kg + hio_canopy_spread_si(io_si) = sites(s)%spread + ipa = 0 cpatch => sites(s)%oldest_patch do while(associated(cpatch)) @@ -1328,7 +1332,7 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_zstar_si_age(io_si,cpatch%age_class) = hio_zstar_si_age(io_si,cpatch%age_class) & + cpatch%zstar * cpatch%area * AREA_INV endif - + ccohort => cpatch%shortest do while(associated(ccohort)) @@ -1391,6 +1395,11 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_biomass_si_pft(io_si, ft) = hio_biomass_si_pft(io_si, ft) + & (ccohort%n * AREA_INV) * ccohort%b * g_per_kg + ! update total biomass per age bin + hio_biomass_si_age(io_si,cpatch%age_class) = hio_biomass_si_age(io_si,cpatch%age_class) & + + ccohort%b * ccohort%n * AREA_INV + + ! Site by Size-Class x PFT (SCPF) ! ------------------------------------------------------------------------ @@ -1453,6 +1462,12 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_m4_si_scpf(io_si,scpf) = hio_m4_si_scpf(io_si,scpf) + ccohort%imort*ccohort%n hio_m5_si_scpf(io_si,scpf) = hio_m5_si_scpf(io_si,scpf) + ccohort%fmort*ccohort%n + + !Y.X. + hio_m7_si_scpf(io_si,scpf) = hio_m7_si_scpf(io_si,scpf) + & + (ccohort%lmort_logging+ccohort%lmort_collateral+ccohort%lmort_infra) * ccohort%n + + ! basal area [m2/ha] hio_ba_si_scpf(io_si,scpf) = hio_ba_si_scpf(io_si,scpf) + & 0.25_r8*3.14159_r8*((dbh/100.0_r8)**2.0_r8)*ccohort%n @@ -1486,8 +1501,14 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_bleaf_canopy_si_scpf(io_si,scpf) = hio_bleaf_canopy_si_scpf(io_si,scpf) + & ccohort%bl * ccohort%n hio_canopy_biomass_pa(io_pa) = hio_canopy_biomass_pa(io_pa) + n_density * ccohort%b * g_per_kg - hio_mortality_canopy_si_scpf(io_si,scpf) = hio_mortality_canopy_si_scpf(io_si,scpf)+ & - (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%fmort) * ccohort%n + + !hio_mortality_canopy_si_scpf(io_si,scpf) = hio_mortality_canopy_si_scpf(io_si,scpf)+ & + ! (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%imort + ccohort%fmort) * ccohort%n + + hio_mortality_canopy_si_scpf(io_si,scpf) = hio_mortality_canopy_si_scpf(io_si,scpf)+ & + (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%imort + ccohort%fmort+ & + ccohort%lmort_logging + ccohort%lmort_collateral + ccohort%lmort_infra) * ccohort%n + hio_nplant_canopy_si_scpf(io_si,scpf) = hio_nplant_canopy_si_scpf(io_si,scpf) + ccohort%n hio_nplant_canopy_si_scls(io_si,scls) = hio_nplant_canopy_si_scls(io_si,scls) + ccohort%n hio_trimming_canopy_si_scls(io_si,scls) = hio_trimming_canopy_si_scls(io_si,scls) + & @@ -1503,13 +1524,18 @@ subroutine update_history_dyn(this,nc,nsites,sites) ccohort%ddbhdt*ccohort%n hio_ddbh_canopy_si_scls(io_si,scls) = hio_ddbh_canopy_si_scls(io_si,scls) + & ccohort%ddbhdt*ccohort%n + ! sum of all mortality hio_mortality_canopy_si_scls(io_si,scls) = hio_mortality_canopy_si_scls(io_si,scls) + & - (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%fmort) * ccohort%n + (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%fmort + & + ccohort%lmort_logging + ccohort%lmort_collateral + ccohort%lmort_infra) * ccohort%n + hio_canopy_mortality_carbonflux_si(io_si) = hio_canopy_mortality_carbonflux_si(io_si) + & (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%fmort) * & - ccohort%b * ccohort%n * g_per_kg * days_per_sec * years_per_day * ha_per_m2 - ! + ccohort%b * ccohort%n * g_per_kg * days_per_sec * years_per_day * ha_per_m2 + & + (ccohort%lmort_logging + ccohort%lmort_collateral + ccohort%lmort_infra)* ccohort%b * & + ccohort%n * g_per_kg * ha_per_m2 + hio_leaf_md_canopy_si_scls(io_si,scls) = hio_leaf_md_canopy_si_scls(io_si,scls) + & ccohort%leaf_md * ccohort%n hio_root_md_canopy_si_scls(io_si,scls) = hio_root_md_canopy_si_scls(io_si,scls) + & @@ -1552,8 +1578,13 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_bleaf_understory_si_scpf(io_si,scpf) = hio_bleaf_understory_si_scpf(io_si,scpf) + & ccohort%bl * ccohort%n hio_understory_biomass_pa(io_pa) = hio_understory_biomass_pa(io_pa) + n_density * ccohort%b * g_per_kg + !hio_mortality_understory_si_scpf(io_si,scpf) = hio_mortality_understory_si_scpf(io_si,scpf)+ & + ! (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%imort + ccohort%fmort) * ccohort%n + hio_mortality_understory_si_scpf(io_si,scpf) = hio_mortality_understory_si_scpf(io_si,scpf)+ & - (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%fmort) * ccohort%n + (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%imort + ccohort%fmort + & + ccohort%lmort_logging + ccohort%lmort_collateral + ccohort%lmort_infra) * ccohort%n + hio_nplant_understory_si_scpf(io_si,scpf) = hio_nplant_understory_si_scpf(io_si,scpf) + ccohort%n hio_nplant_understory_si_scls(io_si,scls) = hio_nplant_understory_si_scls(io_si,scls) + ccohort%n hio_trimming_understory_si_scls(io_si,scls) = hio_trimming_understory_si_scls(io_si,scls) + & @@ -1564,17 +1595,24 @@ subroutine update_history_dyn(this,nc,nsites,sites) n_perm2*ccohort%gpp_acc_hold hio_ar_understory_si_scpf(io_si,scpf) = hio_ar_understory_si_scpf(io_si,scpf) + & n_perm2*ccohort%resp_acc_hold + ! growth increment hio_ddbh_understory_si_scpf(io_si,scpf) = hio_ddbh_understory_si_scpf(io_si,scpf) + & ccohort%ddbhdt*ccohort%n hio_ddbh_understory_si_scls(io_si,scls) = hio_ddbh_understory_si_scls(io_si,scls) + & ccohort%ddbhdt*ccohort%n + ! sum of all mortality hio_mortality_understory_si_scls(io_si,scls) = hio_mortality_understory_si_scls(io_si,scls) + & - (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%fmort) * ccohort%n + (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%imort + ccohort%fmort +& + ccohort%lmort_logging + ccohort%lmort_collateral + ccohort%lmort_infra) * ccohort%n + hio_understory_mortality_carbonflux_si(io_si) = hio_understory_mortality_carbonflux_si(io_si) + & - (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%fmort) * & - ccohort%b * ccohort%n * g_per_kg * days_per_sec * years_per_day * ha_per_m2 + (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%imort + ccohort%fmort) * & + ccohort%b * ccohort%n * g_per_kg * days_per_sec * years_per_day * ha_per_m2 + & + (ccohort%lmort_logging + ccohort%lmort_collateral + ccohort%lmort_infra) * ccohort%b * & + ccohort%n * g_per_kg * ha_per_m2 + ! hio_leaf_md_understory_si_scls(io_si,scls) = hio_leaf_md_understory_si_scls(io_si,scls) + & ccohort%leaf_md * ccohort%n @@ -1609,7 +1647,8 @@ subroutine update_history_dyn(this,nc,nsites,sites) ccohort%canopy_layer_yesterday * ccohort%n endif ! - ! consider imort as understory mortality even if it happens in cohorts that may have been promoted as part of the patch creation... + ! consider imort as understory mortality even if it happens in + ! cohorts that may have been promoted as part of the patch creation... hio_mortality_understory_si_scpf(io_si,scpf) = hio_mortality_understory_si_scpf(io_si,scpf)+ & (ccohort%imort) * ccohort%n hio_mortality_understory_si_scls(io_si,scls) = hio_mortality_understory_si_scls(io_si,scls) + & @@ -1686,8 +1725,6 @@ subroutine update_history_dyn(this,nc,nsites,sites) g_per_kg * patch_scaling_scalar * years_per_day * days_per_sec - hio_canopy_spread_pa(io_pa) = cpatch%spread(1) - do i_cwd = 1, ncwd hio_cwd_ag_si_cwdsc(io_si, i_cwd) = hio_cwd_ag_si_cwdsc(io_si, i_cwd) + & cpatch%CWD_AG(i_cwd)*cpatch%area * AREA_INV * g_per_kg @@ -1708,7 +1745,7 @@ subroutine update_history_dyn(this,nc,nsites,sites) end do !patch loop ! divide so-far-just-summed but to-be-averaged patch-age-class variables by patch-age-class area to get mean values - do ipa2 = 1, nlevage_ed + do ipa2 = 1, nlevage if (hio_area_si_age(io_si, ipa2) .gt. tiny) then hio_lai_si_age(io_si, ipa2) = hio_lai_si_age(io_si, ipa2) / (hio_area_si_age(io_si, ipa2)*AREA) hio_ncl_si_age(io_si, ipa2) = hio_ncl_si_age(io_si, ipa2) / (hio_area_si_age(io_si, ipa2)*AREA) @@ -1721,8 +1758,8 @@ subroutine update_history_dyn(this,nc,nsites,sites) ! pass the cohort termination mortality as a flux to the history, and then reset the termination mortality buffer ! note there are various ways of reporting the total mortality, so pass to these as well do i_pft = 1, numpft - do i_scls = 1,nlevsclass_ed - i_scpf = (i_pft-1)*nlevsclass_ed + i_scls + do i_scls = 1,nlevsclass + i_scpf = (i_pft-1)*nlevsclass + i_scls hio_m6_si_scpf(io_si,i_scpf) = (sites(s)%terminated_nindivs(i_scls,i_pft,1) + & sites(s)%terminated_nindivs(i_scls,i_pft,2)) * days_per_year hio_mortality_canopy_si_scls(io_si,i_scls) = hio_mortality_canopy_si_scls(io_si,i_scls) + & @@ -1745,20 +1782,32 @@ subroutine update_history_dyn(this,nc,nsites,sites) ! summarize all of the mortality fluxes by PFT do i_pft = 1, numpft - do i_scls = 1,nlevsclass_ed - i_scpf = (i_pft-1)*nlevsclass_ed + i_scls + do i_scls = 1,nlevsclass + i_scpf = (i_pft-1)*nlevsclass + i_scls + !hio_mortality_si_pft(io_si,i_pft) = hio_mortality_si_pft(io_si,i_pft) + & + ! hio_m1_si_scpf(io_si,i_scpf) + & + ! hio_m2_si_scpf(io_si,i_scpf) + & + ! hio_m3_si_scpf(io_si,i_scpf) + & + ! hio_m4_si_scpf(io_si,i_scpf) + & + ! hio_m5_si_scpf(io_si,i_scpf) + & + ! hio_m6_si_scpf(io_si,i_scpf) + hio_mortality_si_pft(io_si,i_pft) = hio_mortality_si_pft(io_si,i_pft) + & hio_m1_si_scpf(io_si,i_scpf) + & hio_m2_si_scpf(io_si,i_scpf) + & hio_m3_si_scpf(io_si,i_scpf) + & hio_m4_si_scpf(io_si,i_scpf) + & hio_m5_si_scpf(io_si,i_scpf) + & - hio_m6_si_scpf(io_si,i_scpf) + hio_m6_si_scpf(io_si,i_scpf) + & + hio_m7_si_scpf(io_si,i_scpf) + + + end do end do ! pass demotion rates and associated carbon fluxes to history - do i_scls = 1,nlevsclass_ed + do i_scls = 1,nlevsclass hio_demotion_rate_si_scls(io_si,i_scls) = sites(s)%demotion_rate(i_scls) * days_per_year hio_promotion_rate_si_scls(io_si,i_scls) = sites(s)%promotion_rate(i_scls) * days_per_year end do @@ -1812,9 +1861,8 @@ subroutine update_history_prod(this,nc,nsites,sites,dt_tstep) ed_cohort_type, & ed_patch_type, & AREA, & - AREA_INV, & - nlevage_ed, & - nlevsclass_ed + AREA_INV + use EDTypesMod , only : nclmax, nlevleaf ! ! Arguments @@ -1836,11 +1884,10 @@ subroutine update_history_prod(this,nc,nsites,sites,dt_tstep) integer :: ft ! functional type index real(r8) :: n_density ! individual of cohort per m2. real(r8) :: n_perm2 ! individuals per m2 for the whole column - real(r8) :: patch_area_by_age(nlevage_ed) ! patch area in each bin for normalizing purposes + real(r8) :: patch_area_by_age(nlevage) ! patch area in each bin for normalizing purposes real(r8), parameter :: tiny = 1.e-5_r8 ! some small number integer :: ipa2 ! patch incrementer integer :: cnlfpft_indx, cnlf_indx, ipft, ican, ileaf ! more iterators and indices - type(fates_history_variable_type),pointer :: hvar type(ed_patch_type),pointer :: cpatch type(ed_cohort_type),pointer :: ccohort real(r8) :: per_dt_tstep ! Time step in frequency units (/s) @@ -2134,7 +2181,7 @@ subroutine update_history_prod(this,nc,nsites,sites,dt_tstep) cpatch => cpatch%younger end do !patch loop - do ipa2 = 1, nlevage_ed + do ipa2 = 1, nlevage if (patch_area_by_age(ipa2) .gt. tiny) then hio_gpp_si_age(io_si, ipa2) = hio_gpp_si_age(io_si, ipa2) / (patch_area_by_age(ipa2)) hio_npp_si_age(io_si, ipa2) = hio_npp_si_age(io_si, ipa2) / (patch_area_by_age(ipa2)) @@ -2162,14 +2209,10 @@ subroutine update_history_hydraulics(this,nc,nsites,sites,dt_tstep) use EDtypesMod , only : ed_site_type, & ed_cohort_type, & ed_patch_type, & - AREA, & - nlevage_ed, & - sclass_ed, & - nlevsclass_ed + AREA use FatesHydraulicsMemMod, only : ed_cohort_hydr_type use FatesHydraulicsMemMod, only : nlevsoi_hyd - use EDTypesMod , only : nlevsclass_ed use EDTypesMod , only : maxpft @@ -2191,14 +2234,13 @@ subroutine update_history_hydraulics(this,nc,nsites,sites,dt_tstep) real(r8) :: n_density ! individual of cohort per m2. real(r8) :: n_perm2 ! individuals per m2 for the whole column real(r8), parameter :: tiny = 1.e-5_r8 ! some small number - real(r8) :: ncohort_scpf(nlevsclass_ed*maxpft) ! Bins to count up cohorts counts used in weighting + real(r8) :: ncohort_scpf(nlevsclass*maxpft) ! Bins to count up cohorts counts used in weighting ! should be "hio_nplant_si_scpf" real(r8) :: number_fraction real(r8) :: number_fraction_rate integer :: ipa2 ! patch incrementer integer :: iscpf ! index of the scpf group - type(fates_history_variable_type),pointer :: hvar type(ed_patch_type),pointer :: cpatch type(ed_cohort_type),pointer :: ccohort type(ed_cohort_hydr_type), pointer :: ccohort_hydr @@ -2381,7 +2423,7 @@ subroutine update_history_hydraulics(this,nc,nsites,sites,dt_tstep) end do !patch loop if(hlm_use_ed_st3.eq.ifalse) then - do scpf=1,nlevsclass_ed*numpft + do scpf=1,nlevsclass*numpft if( abs(hio_nplant_si_scpf(io_si, scpf)-ncohort_scpf(scpf)) > 1.0E-8_r8 ) then write(fates_log(),*) 'nplant check on hio_nplant_si_scpf fails during hydraulics history updates' call endrun(msg=errMsg(sourcefile, __LINE__)) @@ -2506,8 +2548,8 @@ subroutine define_history_vars(this, initialize_variables) call this%set_history_var(vname='CANOPY_SPREAD', units='0-1', & long='Scaling factor between tree basal area and canopy area', & use_default='active', & - avgflag='A', vtype=patch_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, & - ivar=ivar, initialize=initialize_variables, index = ih_canopy_spread_pa) + avgflag='A', vtype=site_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, & + ivar=ivar, initialize=initialize_variables, index = ih_canopy_spread_si) call this%set_history_var(vname='PFTbiomass', units='gC/m2', & long='total PFT level biomass', use_default='active', & @@ -2576,6 +2618,12 @@ subroutine define_history_vars(this, initialize_variables) avgflag='A', vtype=site_age_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, & ivar=ivar, initialize=initialize_variables, index = ih_zstar_si_age ) + call this%set_history_var(vname='BIOMASS_BY_AGE', units='m', & + long='Total Biomass within a given patch age bin (kg C)', & + use_default='inactive', & + avgflag='A', vtype=site_age_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, & + ivar=ivar, initialize=initialize_variables, index = ih_biomass_si_age ) + ! Fire Variables call this%set_history_var(vname='FIRE_NESTEROV_INDEX', units='none', & @@ -3136,6 +3184,14 @@ subroutine define_history_vars(this, initialize_variables) avgflag='A', vtype=site_size_pft_r8, hlms='CLM:ALM', flushval=0.0_r8, & upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_m6_si_scpf ) + + !Logging + call this%set_history_var(vname='M7_SCPF', units = 'N/ha/event', & + long='logging mortalities by pft/size',use_default='inactive', & + avgflag='A', vtype=site_size_pft_r8, hlms='CLM:ALM', flushval=0.0_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_m7_si_scpf ) + + call this%set_history_var(vname='MORTALITY_CANOPY_SCPF', units = 'N/ha/yr', & long='total mortality of canopy plants by pft/size', use_default='inactive', & avgflag='A', vtype=site_size_pft_r8, hlms='CLM:ALM', flushval=0.0_r8, & diff --git a/main/FatesHistoryVariableType.F90 b/main/FatesHistoryVariableType.F90 index e76531062c..86b9a4f875 100644 --- a/main/FatesHistoryVariableType.F90 +++ b/main/FatesHistoryVariableType.F90 @@ -89,7 +89,7 @@ subroutine Init(this, vname, units, long, use_default, & call dim_kinds(dk_index)%set_active() call this%GetBounds(0, dim_bounds, dim_kinds, lb1, ub1, lb2, ub2) - + ! NOTE(rgk, 2016-09) currently, all array spaces are flushed each ! time the update is called. The flush here on the initialization ! may be redundant, but will prevent issues in the future if we @@ -167,7 +167,7 @@ subroutine Init(this, vname, units, long, use_default, & end subroutine Init ! ===================================================================================== - + subroutine GetBounds(this, thread, dim_bounds, dim_kinds, lb1, ub1, lb2, ub2) use FatesIODimensionsMod, only : fates_io_dimension_type @@ -176,7 +176,7 @@ subroutine GetBounds(this, thread, dim_bounds, dim_kinds, lb1, ub1, lb2, ub2) class(fates_history_variable_type), intent(inout) :: this integer, intent(in) :: thread - class(fates_io_dimension_type), intent(in) :: dim_bounds(:) + type(fates_io_dimension_type), intent(in) :: dim_bounds(:) type(fates_io_variable_kind_type), intent(in) :: dim_kinds(:) integer, intent(out) :: lb1 integer, intent(out) :: ub1 @@ -205,14 +205,17 @@ subroutine GetBounds(this, thread, dim_bounds, dim_kinds, lb1, ub1, lb2, ub2) ub2 = dim_bounds(d_index)%upper_bound end if else + d_index = dim_kinds(this%dim_kinds_index)%dim1_index lb1 = dim_bounds(d_index)%clump_lower_bound(thread) ub1 = dim_bounds(d_index)%clump_upper_bound(thread) + if(ndims>1)then d_index = dim_kinds(this%dim_kinds_index)%dim2_index lb2 = dim_bounds(d_index)%clump_lower_bound(thread) ub2 = dim_bounds(d_index)%clump_upper_bound(thread) end if + end if end subroutine GetBounds diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index ba4597fd8d..54a7abbb3a 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -19,7 +19,7 @@ module FatesInterfaceMod use EDTypesMod , only : nlevleaf use EDTypesMod , only : maxpft use FatesConstantsMod , only : r8 => fates_r8 - use FatesConstantsMod , only : itrue + use FatesConstantsMod , only : itrue,ifalse use FatesGlobals , only : fates_global_verbose use FatesGlobals , only : fates_log use FatesGlobals , only : endrun => fates_endrun @@ -121,6 +121,10 @@ module FatesInterfaceMod integer, protected :: hlm_use_spitfire ! This flag signals whether or not to use SPITFIRE ! 1 = TRUE, 0 = FALSE + + integer, protected :: hlm_use_logging ! This flag signals whether or not to use + ! the logging module + integer, protected :: hlm_use_planthydro ! This flag signals whether or not to use ! plant hydraulics (bchristo/xu methods) ! 1 = TRUE, 0 = FALSE @@ -235,7 +239,9 @@ module FatesInterfaceMod ! ! ------------------------------------------------------------------------------------- - integer, protected :: numpft ! The total number of PFTs defined in the simulation + integer, protected :: numpft ! The total number of PFTs defined in the simulation + integer, protected :: nlevsclass ! The total number of cohort size class bins output to history + integer, protected :: nlevage ! The total number of patch age bins output to history ! ------------------------------------------------------------------------------------- @@ -841,7 +847,7 @@ subroutine set_fates_global_elements(use_fates) ! This subroutine is called directly from the HLM, and is the first FATES routine ! that is called. ! - ! This subroutine MUST BE CALLED AFTER the FATES parameter file has been read in, + ! This subroutine MUST BE CALLED AFTER the FATES PFT parameter file has been read in, ! and the EDPftvarcon_inst structure has been made. ! This subroutine must ALSO BE CALLED BEFORE the history file dimensions ! are set. @@ -852,13 +858,19 @@ subroutine set_fates_global_elements(use_fates) ! ! -------------------------------------------------------------------------------- - + use EDParamsMod, only : ED_val_history_sizeclass_bin_edges, ED_val_history_ageclass_bin_edges + use CLMFatesParamInterfaceMod , only : FatesReadParameters implicit none logical,intent(in) :: use_fates ! Is fates turned on? + integer :: i + if (use_fates) then + ! first read the non-PFT parameters + call FatesReadParameters() + ! Identify the number of PFTs by evaluating a pft array ! Using wood density as that is not expected to be deprecated any time soon @@ -888,6 +900,33 @@ subroutine set_fates_global_elements(use_fates) fates_maxElementsPerSite = maxPatchesPerSite * fates_maxElementsPerPatch + ! Identify number of size and age class bins for history output + ! assume these arrays are 1-indexed + nlevsclass = size(ED_val_history_sizeclass_bin_edges,dim=1) + nlevage = size(ED_val_history_ageclass_bin_edges,dim=1) + + ! do some checks on the size and age bin arrays to make sure they make sense: + ! make sure that both start at zero, and that both are monotonically increasing + if ( ED_val_history_sizeclass_bin_edges(1) .ne. 0._r8 ) then + write(fates_log(), *) 'size class bins specified in parameter file must start at zero' + call endrun(msg=errMsg(sourcefile, __LINE__)) + endif + if ( ED_val_history_ageclass_bin_edges(1) .ne. 0._r8 ) then + write(fates_log(), *) 'age class bins specified in parameter file must start at zero' + call endrun(msg=errMsg(sourcefile, __LINE__)) + endif + do i = 2,nlevsclass + if ( (ED_val_history_sizeclass_bin_edges(i) - ED_val_history_sizeclass_bin_edges(i-1)) .le. 0._r8) then + write(fates_log(), *) 'size class bins specified in parameter file must be monotonically increasing' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + end do + do i = 2,nlevage + if ( (ED_val_history_ageclass_bin_edges(i) - ED_val_history_ageclass_bin_edges(i-1)) .le. 0._r8) then + write(fates_log(), *) 'age class bins specified in parameter file must be monotonically increasing' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + end do ! Set Various Mapping Arrays used in history output as well ! These will not be used if use_ed or use_fates is false @@ -913,15 +952,12 @@ end subroutine set_fates_global_elements subroutine fates_history_maps - use EDTypesMod, only : nlevsclass_ed use EDTypesMod, only : NFSC use EDTypesMod, only : NCWD - use EDTypesMod, only : nlevage_ed - use EDTypesMod, only : nlevsclass_ed use EDTypesMod, only : nclmax use EDTypesMod, only : nlevleaf - use EDTypesMod, only : sclass_ed - use EDTypesMod, only : ageclass_ed + use EDParamsMod, only : ED_val_history_sizeclass_bin_edges + use EDParamsMod, only : ED_val_history_ageclass_bin_edges ! ------------------------------------------------------------------------------------------ ! This subroutine allocates and populates the variables @@ -940,13 +976,13 @@ subroutine fates_history_maps integer :: ileaf integer :: iage - allocate( fates_hdim_levsclass(1:nlevsclass_ed )) - allocate( fates_hdim_pfmap_levscpf(1:nlevsclass_ed*numpft)) - allocate( fates_hdim_scmap_levscpf(1:nlevsclass_ed*numpft)) + allocate( fates_hdim_levsclass(1:nlevsclass )) + allocate( fates_hdim_pfmap_levscpf(1:nlevsclass*numpft)) + allocate( fates_hdim_scmap_levscpf(1:nlevsclass*numpft)) allocate( fates_hdim_levpft(1:numpft )) allocate( fates_hdim_levfuel(1:NFSC )) allocate( fates_hdim_levcwdsc(1:NCWD )) - allocate( fates_hdim_levage(1:nlevage_ed )) + allocate( fates_hdim_levage(1:nlevage )) allocate( fates_hdim_levcan(nclmax)) allocate( fates_hdim_canmap_levcnlf(nlevleaf*nclmax)) @@ -954,15 +990,12 @@ subroutine fates_history_maps allocate( fates_hdim_canmap_levcnlfpf(nlevleaf*nclmax*numpft)) allocate( fates_hdim_lfmap_levcnlfpf(nlevleaf*nclmax*numpft)) allocate( fates_hdim_pftmap_levcnlfpf(nlevleaf*nclmax*numpft)) - allocate( fates_hdim_scmap_levscag(nlevsclass_ed * nlevage_ed )) - allocate( fates_hdim_agmap_levscag(nlevsclass_ed * nlevage_ed )) + allocate( fates_hdim_scmap_levscag(nlevsclass * nlevage )) + allocate( fates_hdim_agmap_levscag(nlevsclass * nlevage )) ! Fill the IO array of plant size classes - ! For some reason the history files did not like - ! a hard allocation of sclass_ed - fates_hdim_levsclass(:) = sclass_ed(:) - - fates_hdim_levage(:) = ageclass_ed(:) + fates_hdim_levsclass(:) = ED_val_history_sizeclass_bin_edges(:) + fates_hdim_levage(:) = ED_val_history_ageclass_bin_edges(:) ! make pft array do ipft=1,numpft @@ -987,7 +1020,7 @@ subroutine fates_history_maps ! Fill the IO arrays that match pft and size class to their combined array i=0 do ipft=1,numpft - do isc=1,nlevsclass_ed + do isc=1,nlevsclass i=i+1 fates_hdim_pfmap_levscpf(i) = ipft fates_hdim_scmap_levscpf(i) = isc @@ -1004,8 +1037,8 @@ subroutine fates_history_maps end do i=0 - do iage=1,nlevage_ed - do isc=1,nlevsclass_ed + do iage=1,nlevage + do isc=1,nlevsclass i=i+1 fates_hdim_scmap_levscag(i) = isc fates_hdim_agmap_levscag(i) = iage @@ -1120,6 +1153,7 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) hlm_use_vertsoilc = unset_int hlm_use_spitfire = unset_int hlm_use_planthydro = unset_int + hlm_use_logging = unset_int hlm_use_ed_st3 = unset_int hlm_use_ed_prescribed_phys = unset_int hlm_use_inventory_init = unset_int @@ -1161,6 +1195,13 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) call endrun(msg=errMsg(sourcefile, __LINE__)) end if + if ( .not.((hlm_use_logging .eq.1).or.(hlm_use_logging.eq.0)) ) then + if (fates_global_verbose()) then + write(fates_log(), *) 'The FATES namelist use_logging flag must be 0 or 1, exiting' + end if + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + if ( .not.((hlm_use_ed_st3.eq.1).or.(hlm_use_ed_st3.eq.0)) ) then if (fates_global_verbose()) then write(fates_log(), *) 'The FATES namelist stand structure flag must be 0 or 1, exiting' @@ -1392,6 +1433,12 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) write(fates_log(),*) 'Transfering hlm_use_planthydro= ',ival,' to FATES' end if + case('use_logging') + hlm_use_logging = ival + if (fates_global_verbose()) then + write(fates_log(),*) 'Transfering hlm_use_logging= ',ival,' to FATES' + end if + case('use_ed_st3') hlm_use_ed_st3 = ival if (fates_global_verbose()) then diff --git a/main/FatesInventoryInitMod.F90 b/main/FatesInventoryInitMod.F90 index 52a75e480c..d01d9c3d0d 100644 --- a/main/FatesInventoryInitMod.F90 +++ b/main/FatesInventoryInitMod.F90 @@ -77,7 +77,6 @@ subroutine initialize_sites_by_inventory(nsites,sites,bc_in) use EDTypesMod, only : nclmax use EDTypesMod, only : maxpft use EDTypesMod, only : ncwd - use EDParamsMod, only : ED_val_maxspread use EDPatchDynamicsMod, only : create_patch use EDPatchDynamicsMod, only : fuse_patches use EDCohortDynamicsMod, only : fuse_cohorts @@ -104,7 +103,6 @@ subroutine initialize_sites_by_inventory(nsites,sites,bc_in) character(len=line_strlen) :: header_str ! large string for whole lines real(r8) :: age_init ! dummy value for creating a patch real(r8) :: area_init ! dummy value for creating a patch - real(r8) :: spread_init(nclmax) ! dummy value for creating a patch real(r8) :: cwd_ag_init(ncwd) ! dummy value for creating a patch real(r8) :: cwd_bg_init(ncwd) ! dummy value for creating a patch real(r8) :: leaf_litter_init(maxpft) ! dummy value for creating a patch @@ -247,13 +245,12 @@ subroutine initialize_sites_by_inventory(nsites,sites,bc_in) age_init = 0.0_r8 area_init = 0.0_r8 - spread_init(:) = ED_val_maxspread cwd_ag_init(:) = 0.0_r8 cwd_bg_init(:) = 0.0_r8 leaf_litter_init(:) = 0.0_r8 root_litter_init(:) = 0.0_r8 - call create_patch(sites(s), newpatch, age_init, area_init, spread_init, & + call create_patch(sites(s), newpatch, age_init, area_init, & cwd_ag_init, cwd_bg_init, & leaf_litter_init, root_litter_init ) @@ -640,7 +637,7 @@ subroutine set_inventory_edpatch_type1(newpatch,pss_file_unit,ipa,ios,patch_name ! fsn (kg/m2) Fast Soil Nitrogen ! -------------------------------------------------------------------------------------------- - use EDTypesMod, only: get_age_class_index + use FatesSizeAgeTypeIndicesMod, only: get_age_class_index use EDtypesMod, only: AREA use EDTypesMod, only: ncwd use SFParamsMod , only : SF_val_CWD_frac diff --git a/main/FatesParametersInterface.F90 b/main/FatesParametersInterface.F90 index bb0cc1f795..f562b26f71 100644 --- a/main/FatesParametersInterface.F90 +++ b/main/FatesParametersInterface.F90 @@ -29,6 +29,8 @@ module FatesParametersInterface character(len=*), parameter, public :: dimension_name_allpfts = 'fates_allpfts' character(len=*), parameter, public :: dimension_name_variants = 'fates_variants' character(len=*), parameter, public :: dimension_name_hydr_organs = 'fates_hydr_organs' + character(len=*), parameter, public :: dimension_name_history_size_bins = 'fates_history_size_bins' + character(len=*), parameter, public :: dimension_name_history_age_bins = 'fates_history_age_bins' ! Dimensions in the host namespace: character(len=*), parameter, public :: dimension_name_host_allpfts = 'allpfts' diff --git a/main/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index 347bfc9052..44ac2ea58f 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -69,6 +69,7 @@ module FatesRestartInterfaceMod integer, private :: ir_fates_to_bgc_this_ts_si integer, private :: ir_fates_to_bgc_last_ts_si integer, private :: ir_seedrainflux_si + integer, private :: ir_trunk_product_si integer, private :: ir_ncohort_pa integer, private :: ir_balive_co integer, private :: ir_bdead_co @@ -99,6 +100,13 @@ module FatesRestartInterfaceMod integer, private :: ir_cmort_co integer, private :: ir_imort_co integer, private :: ir_fmort_co + + !Logging + integer, private :: ir_lmort_logging_co + integer, private :: ir_lmort_collateral_co + integer, private :: ir_lmort_infra_co + + integer, private :: ir_ddbhdt_co integer, private :: ir_dbalivedt_co integer, private :: ir_dbdeaddt_co @@ -114,7 +122,7 @@ module FatesRestartInterfaceMod integer, private :: ir_leaf_litter_in_paft integer, private :: ir_root_litter_in_paft integer, private :: ir_seed_bank_sift - integer, private :: ir_spread_pacl + integer, private :: ir_spread_si integer, private :: ir_livegrass_pa integer, private :: ir_age_pa integer, private :: ir_area_pa @@ -584,6 +592,12 @@ subroutine define_restart_vars(this, initialize_variables) units='kgC/m2/year', flushval = flushzero, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_seedrainflux_si ) + call this%set_restart_var(vname='fates_trunk_product_site', vtype=site_r8, & + long_name='Accumulate trunk product flux at site', & + units='kgC/m2', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_trunk_product_si ) + + ! ----------------------------------------------------------------------------------- ! Variables stored within cohort vectors ! Note: Some of these are multi-dimensional variables in the patch/site dimension @@ -735,6 +749,23 @@ subroutine define_restart_vars(this, initialize_variables) units='/year', flushval = flushzero, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_fmort_co ) + + call this%set_restart_var(vname='fates_lmort_logging', vtype=cohort_r8, & + long_name='ed cohort - directly logging mortality rate', & + units='%/event', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_lmort_logging_co ) + + call this%set_restart_var(vname='fates_lmort_collateral', vtype=cohort_r8, & + long_name='ed cohort - collateral mortality rate', & + units='%/event', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_lmort_collateral_co ) + + call this%set_restart_var(vname='fates_lmort_in', vtype=cohort_r8, & + long_name='ed cohort - mechanical mortality rate', & + units='%/event', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_lmort_infra_co ) + + call this%set_restart_var(vname='fates_ddbhdt', vtype=cohort_r8, & long_name='ed cohort - differential: ddbh/dt', & units='cm/year', flushval = flushzero, & @@ -812,10 +843,10 @@ subroutine define_restart_vars(this, initialize_variables) units='kgC/m2/year', flushval = flushzero, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_seed_bank_sift ) - call this%set_restart_var(vname='fates_spread', vtype=cohort_r8, & + call this%set_restart_var(vname='fates_spread', vtype=site_r8, & long_name='dynamic ratio of dbh to canopy area, by patch x canopy-layer', & units='cm/m2', flushval = flushzero, & - hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_spread_pacl ) + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_spread_si ) call this%set_restart_var(vname='fates_livegrass', vtype=cohort_r8, & long_name='total AGB from grass, by patch', & @@ -961,7 +992,6 @@ subroutine set_restart_vectors(this,nc,nsites,sites) integer :: io_idx_co ! cohort index integer :: io_idx_pa_pft ! each pft within each patch (pa_pft) integer :: io_idx_pa_cwd ! each cwd class within each patch (pa_cwd) - integer :: io_idx_pa_cl ! each canopy layer class within each patch (pa_cl) integer :: io_idx_pa_sunz ! index for the combined dimensions for radiation integer :: io_idx_si_wmem ! each water memory class within each site @@ -1001,6 +1031,7 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_fates_to_bgc_this_ts_si => this%rvars(ir_fates_to_bgc_this_ts_si)%r81d, & rio_fates_to_bgc_last_ts_si => this%rvars(ir_fates_to_bgc_last_ts_si)%r81d, & rio_seedrainflux_si => this%rvars(ir_seedrainflux_si)%r81d, & + rio_trunk_product_si => this%rvars(ir_trunk_product_si)%r81d, & rio_ncohort_pa => this%rvars(ir_ncohort_pa)%int1d, & rio_balive_co => this%rvars(ir_balive_co)%r81d, & rio_bdead_co => this%rvars(ir_bdead_co)%r81d, & @@ -1031,6 +1062,14 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_cmort_co => this%rvars(ir_cmort_co)%r81d, & rio_imort_co => this%rvars(ir_imort_co)%r81d, & rio_fmort_co => this%rvars(ir_fmort_co)%r81d, & + + + + rio_lmort_logging_co => this%rvars(ir_lmort_logging_co)%r81d, & + rio_lmort_collateral_co => this%rvars(ir_lmort_collateral_co)%r81d, & + rio_lmort_infra_co => this%rvars(ir_lmort_infra_co)%r81d, & + + rio_ddbhdt_co => this%rvars(ir_ddbhdt_co)%r81d, & rio_dbalivedt_co => this%rvars(ir_dbalivedt_co)%r81d, & rio_dbdeaddt_co => this%rvars(ir_dbdeaddt_co)%r81d, & @@ -1046,7 +1085,7 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_leaf_litter_in_paft => this%rvars(ir_leaf_litter_in_paft)%r81d, & rio_root_litter_in_paft => this%rvars(ir_root_litter_in_paft)%r81d, & rio_seed_bank_sift => this%rvars(ir_seed_bank_sift)%r81d, & - rio_spread_pacl => this%rvars(ir_spread_pacl)%r81d, & + rio_spread_si => this%rvars(ir_spread_si)%r81d, & rio_livegrass_pa => this%rvars(ir_livegrass_pa)%r81d, & rio_age_pa => this%rvars(ir_age_pa)%r81d, & rio_area_pa => this%rvars(ir_area_pa)%r81d, & @@ -1078,7 +1117,6 @@ subroutine set_restart_vectors(this,nc,nsites,sites) io_idx_co = io_idx_co_1st io_idx_pa_pft = io_idx_co_1st io_idx_pa_cwd = io_idx_co_1st - io_idx_pa_cl = io_idx_co_1st io_idx_si_wmem = io_idx_co_1st io_idx_pa_sunz = io_idx_co_1st @@ -1086,6 +1124,9 @@ subroutine set_restart_vectors(this,nc,nsites,sites) do i = 1,numpft rio_seed_bank_sift(io_idx_co_1st+i-1) = sites(s)%seed_bank(i) end do + + ! canopy spread term + rio_spread_si(io_idx_si) = sites(s)%spread cpatch => sites(s)%oldest_patch @@ -1143,6 +1184,14 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_cmort_co(io_idx_co) = ccohort%cmort rio_imort_co(io_idx_co) = ccohort%imort rio_fmort_co(io_idx_co) = ccohort%fmort + + !Logging + rio_lmort_logging_co(io_idx_co) = ccohort%lmort_logging + rio_lmort_collateral_co(io_idx_co) = ccohort%lmort_collateral + rio_lmort_infra_co(io_idx_co) = ccohort%lmort_infra + + + rio_ddbhdt_co(io_idx_co) = ccohort%ddbhdt rio_dbalivedt_co(io_idx_co) = ccohort%dbalivedt rio_dbdeaddt_co(io_idx_co) = ccohort%dbdeaddt @@ -1200,11 +1249,6 @@ subroutine set_restart_vectors(this,nc,nsites,sites) io_idx_pa_cwd = io_idx_pa_cwd + 1 end do - do i = 1,nclmax ! nclmax currently 2 - rio_spread_pacl(io_idx_pa_cl) = cpatch%spread(i) - io_idx_pa_cl = io_idx_pa_cl + 1 - end do - if ( DEBUG ) write(fates_log(),*) 'CLTV io_idx_pa_sunz 1 ',io_idx_pa_sunz if ( DEBUG ) write(fates_log(),*) 'CLTV 1186 ',nlevleaf,numpft,nclmax @@ -1232,7 +1276,6 @@ subroutine set_restart_vectors(this,nc,nsites,sites) ! reset counters so that they are all advanced evenly. io_idx_pa_pft = io_idx_co_1st io_idx_pa_cwd = io_idx_co_1st - io_idx_pa_cl = io_idx_co_1st io_idx_co = io_idx_co_1st io_idx_pa_sunz = io_idx_co_1st @@ -1270,8 +1313,11 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_fates_to_bgc_this_ts_si(io_idx_si) = sites(s)%fates_to_bgc_this_ts rio_fates_to_bgc_last_ts_si(io_idx_si) = sites(s)%fates_to_bgc_last_ts rio_seedrainflux_si(io_idx_si) = sites(s)%tot_seed_rain_flux - + + ! Accumulated trunk product + rio_trunk_product_si(io_idx_si) = sites(s)%resources_management%trunk_product_site ! set numpatches for this column + rio_npatch_si(io_idx_si) = patchespersite do i = 1,numWaterMem ! numWaterMem currently 10 @@ -1312,7 +1358,7 @@ subroutine create_patchcohort_structure(this, nc, nsites, sites, bc_in) use EDGrowthFunctionsMod, only : Dbh use EDCohortDynamicsMod, only : create_cohort use EDInitMod, only : zero_site - use EDParamsMod, only : ED_val_maxspread + use EDInitMod, only : init_site_vars use EDPatchDynamicsMod, only : create_patch use EDPftvarcon, only : EDPftvarcon_inst @@ -1329,7 +1375,6 @@ subroutine create_patchcohort_structure(this, nc, nsites, sites, bc_in) type(ed_cohort_type), allocatable :: temp_cohort real(r8) :: cwd_ag_local(ncwd) real(r8) :: cwd_bg_local(ncwd) - real(r8) :: spread_local(nclmax) real(r8) :: leaf_litter_local(maxpft) real(r8) :: root_litter_local(maxpft) real(r8) :: patch_age @@ -1348,7 +1393,6 @@ subroutine create_patchcohort_structure(this, nc, nsites, sites, bc_in) cwd_bg_local(:) = 0.0_r8 leaf_litter_local(:) = 0.0_r8 root_litter_local(:) = 0.0_r8 - spread_local(:) = ED_val_maxspread patch_age = 0.0_r8 ! ---------------------------------------------------------------------------------- @@ -1365,6 +1409,7 @@ subroutine create_patchcohort_structure(this, nc, nsites, sites, bc_in) io_idx_si = this%restart_map(nc)%site_index(s) io_idx_co_1st = this%restart_map(nc)%cohort1_index(s) + call init_site_vars( sites(s) ) call zero_site( sites(s) ) ! @@ -1396,7 +1441,7 @@ subroutine create_patchcohort_structure(this, nc, nsites, sites, bc_in) ! make new patch call create_patch(sites(s), newp, patch_age, area, & - spread_local, cwd_ag_local, cwd_bg_local, & + cwd_ag_local, cwd_bg_local, & leaf_litter_local, root_litter_local) newp%siteptr => sites(s) @@ -1537,7 +1582,6 @@ subroutine get_restart_vectors(this, nc, nsites, sites) integer :: io_idx_co ! cohort index integer :: io_idx_pa_pft ! each pft within each patch (pa_pft) integer :: io_idx_pa_cwd ! each cwd class within each patch (pa_cwd) - integer :: io_idx_pa_cl ! each canopy layer class within each patch (pa_cl) integer :: io_idx_pa_sunz ! index for the combined dimensions for radiation integer :: io_idx_si_wmem ! each water memory class within each site @@ -1571,6 +1615,7 @@ subroutine get_restart_vectors(this, nc, nsites, sites) rio_fates_to_bgc_this_ts_si => this%rvars(ir_fates_to_bgc_this_ts_si)%r81d, & rio_fates_to_bgc_last_ts_si => this%rvars(ir_fates_to_bgc_last_ts_si)%r81d, & rio_seedrainflux_si => this%rvars(ir_seedrainflux_si)%r81d, & + rio_trunk_product_si => this%rvars(ir_trunk_product_si)%r81d, & rio_ncohort_pa => this%rvars(ir_ncohort_pa)%int1d, & rio_balive_co => this%rvars(ir_balive_co)%r81d, & rio_bdead_co => this%rvars(ir_bdead_co)%r81d, & @@ -1601,6 +1646,13 @@ subroutine get_restart_vectors(this, nc, nsites, sites) rio_cmort_co => this%rvars(ir_cmort_co)%r81d, & rio_imort_co => this%rvars(ir_imort_co)%r81d, & rio_fmort_co => this%rvars(ir_fmort_co)%r81d, & + + rio_lmort_logging_co => this%rvars(ir_lmort_logging_co)%r81d, & + rio_lmort_collateral_co => this%rvars(ir_lmort_collateral_co)%r81d, & + rio_lmort_infra_co => this%rvars(ir_lmort_infra_co)%r81d, & + + + rio_ddbhdt_co => this%rvars(ir_ddbhdt_co)%r81d, & rio_dbalivedt_co => this%rvars(ir_dbalivedt_co)%r81d, & rio_dbdeaddt_co => this%rvars(ir_dbdeaddt_co)%r81d, & @@ -1616,7 +1668,7 @@ subroutine get_restart_vectors(this, nc, nsites, sites) rio_leaf_litter_in_paft => this%rvars(ir_leaf_litter_in_paft)%r81d, & rio_root_litter_in_paft => this%rvars(ir_root_litter_in_paft)%r81d, & rio_seed_bank_sift => this%rvars(ir_seed_bank_sift)%r81d, & - rio_spread_pacl => this%rvars(ir_spread_pacl)%r81d, & + rio_spread_si => this%rvars(ir_spread_si)%r81d, & rio_livegrass_pa => this%rvars(ir_livegrass_pa)%r81d, & rio_age_pa => this%rvars(ir_age_pa)%r81d, & rio_area_pa => this%rvars(ir_area_pa)%r81d, & @@ -1637,7 +1689,6 @@ subroutine get_restart_vectors(this, nc, nsites, sites) io_idx_co = io_idx_co_1st io_idx_pa_pft = io_idx_co_1st io_idx_pa_cwd = io_idx_co_1st - io_idx_pa_cl = io_idx_co_1st io_idx_pa_sunz = io_idx_co_1st io_idx_si_wmem = io_idx_co_1st @@ -1645,6 +1696,8 @@ subroutine get_restart_vectors(this, nc, nsites, sites) do i = 1,numpft sites(s)%seed_bank(i) = rio_seed_bank_sift(io_idx_co_1st+i-1) enddo + + sites(s)%spread = rio_spread_si(io_idx_si) ! Perform a check on the number of patches per site patchespersite = 0 @@ -1698,6 +1751,14 @@ subroutine get_restart_vectors(this, nc, nsites, sites) ccohort%cmort = rio_cmort_co(io_idx_co) ccohort%imort = rio_imort_co(io_idx_co) ccohort%fmort = rio_fmort_co(io_idx_co) + + !Logging + ccohort%lmort_logging = rio_lmort_logging_co(io_idx_co) + ccohort%lmort_collateral = rio_lmort_collateral_co(io_idx_co) + ccohort%lmort_infra = rio_lmort_infra_co(io_idx_co) + + + ccohort%ddbhdt = rio_ddbhdt_co(io_idx_co) ccohort%dbalivedt = rio_dbalivedt_co(io_idx_co) ccohort%dbdeaddt = rio_dbdeaddt_co(io_idx_co) @@ -1725,7 +1786,6 @@ subroutine get_restart_vectors(this, nc, nsites, sites) cpatch%root_litter(:) = 0.0_r8 cpatch%leaf_litter_in(:) = 0.0_r8 cpatch%root_litter_in(:) = 0.0_r8 - cpatch%spread(:) = 0.0_r8 ! ! deal with patch level fields here @@ -1760,11 +1820,6 @@ subroutine get_restart_vectors(this, nc, nsites, sites) io_idx_pa_cwd = io_idx_pa_cwd + 1 enddo - do i = 1,nclmax ! nclmax currently 2 - cpatch%spread(i) = rio_spread_pacl(io_idx_pa_cl) - io_idx_pa_cl = io_idx_pa_cl + 1 - enddo - if ( DEBUG ) write(fates_log(),*) 'CVTL io_idx_pa_sunz 1 ',io_idx_pa_sunz do k = 1,nlevleaf ! nlevleaf currently 40 @@ -1790,7 +1845,6 @@ subroutine get_restart_vectors(this, nc, nsites, sites) ! and max the number of allowed cohorts per patch io_idx_pa_pft = io_idx_co_1st io_idx_pa_cwd = io_idx_co_1st - io_idx_pa_cl = io_idx_co_1st io_idx_co = io_idx_co_1st io_idx_pa_sunz = io_idx_co_1st @@ -1838,7 +1892,8 @@ subroutine get_restart_vectors(this, nc, nsites, sites) sites(s)%fates_to_bgc_this_ts = rio_fates_to_bgc_this_ts_si(io_idx_si) sites(s)%fates_to_bgc_last_ts = rio_fates_to_bgc_last_ts_si(io_idx_si) sites(s)%tot_seed_rain_flux = rio_seedrainflux_si(io_idx_si) - + sites(s)%resources_management%trunk_product_site = rio_trunk_product_si(io_idx_si) + end do if ( DEBUG ) then diff --git a/main/FatesRestartVariableType.F90 b/main/FatesRestartVariableType.F90 index 40648fb4c6..48152adf7f 100644 --- a/main/FatesRestartVariableType.F90 +++ b/main/FatesRestartVariableType.F90 @@ -118,7 +118,7 @@ subroutine GetBounds(this, thread, dim_bounds, dim_kinds, lb1, ub1, lb2, ub2) class(fates_restart_variable_type), intent(inout) :: this integer, intent(in) :: thread - class(fates_io_dimension_type), intent(in) :: dim_bounds(:) + type(fates_io_dimension_type), intent(in) :: dim_bounds(:) type(fates_io_variable_kind_type), intent(in) :: dim_kinds(:) integer, intent(out) :: lb1 integer, intent(out) :: ub1 diff --git a/main/FatesSizeAgeTypeIndicesMod.F90 b/main/FatesSizeAgeTypeIndicesMod.F90 new file mode 100644 index 0000000000..bac132ba12 --- /dev/null +++ b/main/FatesSizeAgeTypeIndicesMod.F90 @@ -0,0 +1,73 @@ +module FatesSizeAgeTypeIndicesMod + + use FatesConstantsMod, only : r8 => fates_r8 + use FatesInterfaceMod, only : nlevsclass + use EDParamsMod, only : ED_val_history_sizeclass_bin_edges + use EDParamsMod, only : ED_val_history_ageclass_bin_edges + + implicit none + +contains + + ! ===================================================================================== + + function get_age_class_index(age) result( patch_age_class ) + + real(r8), intent(in) :: age + + integer :: patch_age_class + + patch_age_class = count(age-ED_val_history_ageclass_bin_edges.ge.0.0_r8) + + end function get_age_class_index + + ! ===================================================================================== + + function get_sizeage_class_index(dbh,age) result(size_by_age_class) + + ! Arguments + real(r8),intent(in) :: dbh + real(r8),intent(in) :: age + + integer :: size_class + integer :: age_class + integer :: size_by_age_class + + size_class = get_size_class_index(dbh) + + age_class = get_age_class_index(age) + + size_by_age_class = (age_class-1)*nlevsclass + size_class + + end function get_sizeage_class_index + + ! ===================================================================================== + + subroutine sizetype_class_index(dbh,pft,size_class,size_by_pft_class) + + ! Arguments + real(r8),intent(in) :: dbh + integer,intent(in) :: pft + integer,intent(out) :: size_class + integer,intent(out) :: size_by_pft_class + + size_class = get_size_class_index(dbh) + + size_by_pft_class = (pft-1)*nlevsclass+size_class + + return + end subroutine sizetype_class_index + + ! ===================================================================================== + + function get_size_class_index(dbh) result(cohort_size_class) + + real(r8), intent(in) :: dbh + + integer :: cohort_size_class + + cohort_size_class = count(dbh-ED_val_history_sizeclass_bin_edges.ge.0.0_r8) + + end function get_size_class_index + +end module FatesSizeAgeTypeIndicesMod