diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index e98549286e..088ca1eabf 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -22,6 +22,7 @@ module EDCohortDynamicsMod use EDTypesMod , only : min_npm2, min_nppatch use EDTypesMod , only : min_n_safemath use EDTypesMod , only : nlevleaf + use EDTypesMod , only : ican_upper use FatesInterfaceMod , only : hlm_use_planthydro use FatesInterfaceMod , only : hlm_parteh_mode use FatesPlantHydraulicsMod, only : FuseCohortHydraulics @@ -633,21 +634,25 @@ subroutine terminate_cohorts( currentSite, currentPatch, level ) endif ! if (.not.currentCohort%isnew .and. level == 2) then if (terminate == 1) then + ! preserve a record of the to-be-terminated cohort for mortality accounting - if (currentCohort%canopy_layer .eq. 1) then - levcan = 1 - else - levcan = 2 - endif + levcan = currentCohort%canopy_layer - currentSite%terminated_nindivs(currentCohort%size_class,currentCohort%pft,levcan) = & - currentSite%terminated_nindivs(currentCohort%size_class,currentCohort%pft,levcan) + currentCohort%n - ! - currentSite%termination_carbonflux(levcan) = currentSite%termination_carbonflux(levcan) + & - currentCohort%n * (struct_c+sapw_c+leaf_c+fnrt_c+store_c+repro_c) - - if( hlm_use_planthydro == itrue ) then + if( hlm_use_planthydro == itrue ) & call AccumulateMortalityWaterStorage(currentSite,currentCohort,currentCohort%n) + + if(levcan==ican_upper) then + currentSite%term_nindivs_canopy(currentCohort%size_class,currentCohort%pft) = & + currentSite%term_nindivs_canopy(currentCohort%size_class,currentCohort%pft) + currentCohort%n + + currentSite%term_carbonflux_canopy = currentSite%term_carbonflux_canopy + & + currentCohort%n * (struct_c+sapw_c+leaf_c+fnrt_c+store_c+repro_c) + else + currentSite%term_nindivs_ustory(currentCohort%size_class,currentCohort%pft) = & + currentSite%term_nindivs_ustory(currentCohort%size_class,currentCohort%pft) + currentCohort%n + + currentSite%term_carbonflux_ustory = currentSite%term_carbonflux_ustory + & + currentCohort%n * (struct_c+sapw_c+leaf_c+fnrt_c+store_c+repro_c) end if !put the litter from the terminated cohorts straight into the fragmenting pools @@ -965,7 +970,6 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) currentCohort%cmort = (currentCohort%n*currentCohort%cmort + nextc%n*nextc%cmort)/newn currentCohort%hmort = (currentCohort%n*currentCohort%hmort + nextc%n*nextc%hmort)/newn currentCohort%bmort = (currentCohort%n*currentCohort%bmort + nextc%n*nextc%bmort)/newn - currentCohort%fmort = (currentCohort%n*currentCohort%fmort + nextc%n*nextc%fmort)/newn currentCohort%frmort = (currentCohort%n*currentCohort%frmort + nextc%n*nextc%frmort)/newn ! logging mortality, Yi Xu @@ -1352,7 +1356,6 @@ subroutine copy_cohort( currentCohort,copyc ) ! Mortality diagnostics n%cmort = o%cmort n%bmort = o%bmort - n%fmort = o%fmort n%hmort = o%hmort n%frmort = o%frmort diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 98546557e3..c7f059ca9d 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -17,6 +17,7 @@ module EDPatchDynamicsMod use EDTypesMod , only : dtype_ifall use EDTypesMod , only : dtype_ilog use EDTypesMod , only : dtype_ifire + use EDTypesMod , only : ican_upper use FatesInterfaceMod , only : hlm_use_planthydro use FatesInterfaceMod , only : hlm_numSWb use FatesInterfaceMod , only : bc_in_type @@ -139,7 +140,6 @@ subroutine disturbance_rates( site_in, bc_in) currentCohort%bmort = bmort currentCohort%hmort = hmort currentCohort%frmort = frmort - currentCohort%fmort = 0.0_r8 ! Fire mortality is initialized as zero, but may be changed call LoggingMortality_frac(currentCohort%pft, currentCohort%dbh, & lmort_direct,lmort_collateral,lmort_infra ) @@ -216,7 +216,6 @@ subroutine disturbance_rates( site_in, bc_in) currentCohort => currentPatch%shortest do while(associated(currentCohort)) if(currentCohort%canopy_layer == 1)then - 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) @@ -267,7 +266,6 @@ subroutine disturbance_rates( site_in, bc_in) currentCohort%lmort_direct = 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 @@ -323,12 +321,14 @@ 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 + integer :: levcan ! canopy level real(r8) :: leaf_c ! leaf carbon [kg] - real(r8) :: fnrt_c ! fineroot carbon [kg] - real(r8) :: sapw_c ! sapwood carbon [kg] - real(r8) :: store_c ! storage carbon [kg] - real(r8) :: struct_c ! structure carbon [kg] - real(r8) :: total_c ! total carbon of plant [kg] + real(r8) :: fnrt_c ! fineroot carbon [kg] + real(r8) :: sapw_c ! sapwood carbon [kg] + real(r8) :: store_c ! storage carbon [kg] + real(r8) :: struct_c ! structure carbon [kg] + real(r8) :: total_c ! total carbon of plant [kg] + !--------------------------------------------------------------------- storesmallcohort => null() ! storage of the smallest cohort for insertion routine @@ -351,12 +351,17 @@ subroutine spawn_patches( currentSite, bc_in) call endrun(msg=errMsg(sourcefile, __LINE__)) end if - site_areadis = site_areadis + currentPatch%area * currentPatch%disturbance_rate + ! Only create new patches that have non-negligible amount of land + if((currentPatch%area*currentPatch%disturbance_rate) > nearzero ) then + site_areadis = site_areadis + currentPatch%area * currentPatch%disturbance_rate + end if + currentPatch => currentPatch%older enddo ! end loop over patches. sum area disturbed for all patches. - if (site_areadis > 0.0_r8) then + if (site_areadis > nearzero) then + cwd_ag_local = 0.0_r8 cwd_bg_local = 0.0_r8 leaf_litter_local = 0.0_r8 @@ -378,6 +383,8 @@ subroutine spawn_patches( currentSite, bc_in) ! This is the amount of patch area that is disturbed, and donated by the donor patch_site_areadis = currentPatch%area * currentPatch%disturbance_rate + if (patch_site_areadis > nearzero) then + call average_patch_properties(currentPatch, new_patch, patch_site_areadis) if (currentPatch%disturbance_rates(dtype_ilog) > currentPatch%disturbance_rates(dtype_ifall) .and. & @@ -442,7 +449,6 @@ subroutine spawn_patches( currentSite, bc_in) 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%frmort = nan nc%lmort_direct = nan nc%lmort_collateral = nan @@ -485,7 +491,6 @@ 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%cmort = currentCohort%cmort nc%hmort = currentCohort%hmort nc%bmort = currentCohort%bmort @@ -510,7 +515,6 @@ subroutine spawn_patches( currentSite, bc_in) ! Those remaining in the existing currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area) - nc%fmort = 0.0_r8 nc%cmort = currentCohort%cmort nc%hmort = currentCohort%hmort nc%bmort = currentCohort%bmort @@ -533,11 +537,40 @@ subroutine spawn_patches( currentSite, bc_in) ! loss of individuals from source patch due to area shrinking currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area) + levcan = currentCohort%canopy_layer + + if(levcan==ican_upper) then + + ! before changing number densities, track total rate of trees that died + ! due to fire, as well as from each fire mortality term + currentSite%fmort_rate_canopy(currentCohort%size_class, currentCohort%pft) = & + currentSite%fmort_rate_canopy(currentCohort%size_class, currentCohort%pft) + & + nc%n * currentCohort%fire_mort / hlm_freq_day + + currentSite%fmort_carbonflux_canopy = currentSite%fmort_carbonflux_canopy + & + (nc%n * currentCohort%fire_mort) * & + total_c * g_per_kg * days_per_sec * ha_per_m2 + + else + currentSite%fmort_rate_ustory(currentCohort%size_class, currentCohort%pft) = & + currentSite%fmort_rate_ustory(currentCohort%size_class, currentCohort%pft) + & + nc%n * currentCohort%fire_mort / hlm_freq_day + + currentSite%fmort_carbonflux_ustory = currentSite%fmort_carbonflux_ustory + & + (nc%n * currentCohort%fire_mort) * & + total_c * g_per_kg * days_per_sec * ha_per_m2 + end if + + currentSite%fmort_rate_cambial(currentCohort%size_class, currentCohort%pft) = & + currentSite%fmort_rate_cambial(currentCohort%size_class, currentCohort%pft) + & + nc%n * currentCohort%cambial_mort / hlm_freq_day + currentSite%fmort_rate_crown(currentCohort%size_class, currentCohort%pft) = & + currentSite%fmort_rate_crown(currentCohort%size_class, currentCohort%pft) + & + nc%n * currentCohort%crownfire_mort / hlm_freq_day + ! 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%cmort = currentCohort%cmort nc%hmort = currentCohort%hmort nc%bmort = currentCohort%bmort @@ -546,7 +579,8 @@ subroutine spawn_patches( currentSite, bc_in) nc%lmort_direct = currentCohort%lmort_direct 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 @@ -567,7 +601,6 @@ subroutine spawn_patches( currentSite, bc_in) nc%hmort = nan nc%bmort = nan nc%frmort = nan - nc%fmort = nan nc%lmort_direct = nan nc%lmort_collateral = nan nc%lmort_infra = nan @@ -609,7 +642,6 @@ subroutine spawn_patches( currentSite, bc_in) currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area) - nc%fmort = 0.0_r8 nc%cmort = currentCohort%cmort nc%hmort = currentCohort%hmort nc%bmort = currentCohort%bmort @@ -629,7 +661,6 @@ subroutine spawn_patches( currentSite, bc_in) 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%cmort = currentCohort%cmort nc%hmort = currentCohort%hmort nc%bmort = currentCohort%bmort @@ -682,10 +713,6 @@ subroutine spawn_patches( currentSite, bc_in) enddo ! currentCohort call sort_cohorts(currentPatch) - !zero disturbance accumulators - currentPatch%disturbance_rate = 0._r8 - currentPatch%disturbance_rates = 0._r8 - !update area of donor patch currentPatch%area = currentPatch%area - patch_site_areadis @@ -698,6 +725,12 @@ subroutine spawn_patches( currentSite, bc_in) call terminate_cohorts(currentSite, currentPatch, 2) call sort_cohorts(currentPatch) + end if ! if (patch_site_areadis > nearzero) then + + !zero disturbance rate trackers + currentPatch%disturbance_rate = 0._r8 + currentPatch%disturbance_rates = 0._r8 + currentPatch => currentPatch%younger enddo ! currentPatch patch loop. @@ -723,9 +756,11 @@ subroutine spawn_patches( currentSite, bc_in) endif !end new_patch area + call check_patch_area(currentSite) call set_patchno(currentSite) + return end subroutine spawn_patches ! ============================================================================ diff --git a/biogeophys/FatesPlantRespPhotosynthMod.F90 b/biogeophys/FatesPlantRespPhotosynthMod.F90 index 2fae4eceb2..4b181ed7b3 100644 --- a/biogeophys/FatesPlantRespPhotosynthMod.F90 +++ b/biogeophys/FatesPlantRespPhotosynthMod.F90 @@ -938,7 +938,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in if ( parsun_lsl <= 0._r8 ) then ! night time - anet_av_out = 0._r8 + anet_av_out = -lmr psn_out = 0._r8 rstoma_out = min(rsmax0, 1._r8/bbb * cf) diff --git a/fire/SFMainMod.F90 b/fire/SFMainMod.F90 index 2a7a1ca01c..505437634d 100644 --- a/fire/SFMainMod.F90 +++ b/fire/SFMainMod.F90 @@ -158,14 +158,15 @@ end subroutine fire_danger_index subroutine charecteristics_of_fuel ( currentSite ) !***************************************************************** - use SFParamsMod, only : SF_val_alpha_FMC, SF_val_SAV, SF_val_FBD + use SFParamsMod, only : SF_val_drying_ratio, SF_val_SAV, SF_val_FBD type(ed_site_type), intent(in), target :: currentSite type(ed_patch_type), pointer :: currentPatch type(ed_cohort_type), pointer :: currentCohort - real(r8) timeav_swc + real(r8) timeav_swc + real(r8) alpha_FMC(nfsc) ! Relative fuel moisture adjusted per drying ratio real(r8) fuel_moisture(nfsc) ! Scaled moisture content of small litter fuels. real(r8) MEF(nfsc) ! Moisture extinction factor of fuels integer n @@ -235,13 +236,15 @@ subroutine charecteristics_of_fuel ( currentSite ) ! Equation 6 in Thonicke et al. 2010. across leaves,twig, small branch, and large branch ! dead leaves and twigs included in 1hr pool per Thonicke (2010) ! Calculate fuel moisture for trunks to hold value for fuel consumption - fuel_moisture(dl_sf:tr_sf) = exp(-1.0_r8 * SF_val_alpha_FMC(dl_sf:tr_sf) * currentSite%acc_NI) + alpha_FMC(dl_sf:tr_sf) = SF_val_SAV(dl_sf:tr_sf)/SF_val_drying_ratio + + fuel_moisture(dl_sf:tr_sf) = exp(-1.0_r8 * alpha_FMC(dl_sf:tr_sf) * currentSite%acc_NI) if(write_SF == itrue)then if ( hlm_masterproc == itrue ) write(fates_log(),*) 'ff3 ',currentPatch%fuel_frac if ( hlm_masterproc == itrue ) write(fates_log(),*) 'fm ',fuel_moisture if ( hlm_masterproc == itrue ) write(fates_log(),*) 'csa ',currentSite%acc_NI - if ( hlm_masterproc == itrue ) write(fates_log(),*) 'sfv ',SF_val_alpha_FMC + if ( hlm_masterproc == itrue ) write(fates_log(),*) 'sfv ',alpha_FMC endif ! FIX(RF,032414): needs refactoring. ! average water content !is this the correct metric? @@ -408,8 +411,8 @@ subroutine rate_of_spread ( currentSite ) use SFParamsMod, only : SF_val_miner_total, & SF_val_part_dens, & SF_val_miner_damp, & - SF_val_fuel_energy, & - SF_val_wind_max + SF_val_fuel_energy + use FatesInterfaceMod, only : hlm_current_day, hlm_current_month type(ed_site_type), intent(in), target :: currentSite @@ -426,7 +429,6 @@ subroutine rate_of_spread ( currentSite ) real(r8) beta_ratio ! ratio of beta/beta_op real(r8) a_beta ! dummy variable for product of a* beta_ratio for react_v_opt equation real(r8) a,b,c,e ! function of fuel sav - real(r8) wind_elev_fire !wind speed (m/min) at elevevation relevant for fire logical,parameter :: debug_windspeed = .false. !for debugging @@ -493,21 +495,9 @@ subroutine rate_of_spread ( currentSite ) ! Equation A5 in Thonicke et al. 2010 ! phi_wind (unitless) - ! convert wind_elev_fire from m/min to ft/min for Rothermel ROS eqn - ! wind max per Lasslop et al 2014 to linearly reduce ROS for high wind speeds - !OLD! phi_wind = c * ((3.281_r8*currentPatch%effect_wspeed)**b)*(beta_ratio**(-e)) - if (currentPatch%effect_wspeed .le. SF_val_wind_max) then - wind_elev_fire = currentPatch%effect_wspeed - phi_wind = c * ((3.281_r8*wind_elev_fire)**b)*(beta_ratio**(-e)) - if (debug_windspeed) write(fates_log(),*) 'SF wind LESS max ', currentPatch%effect_wspeed - if (debug_windspeed) write(fates_log(),*) 'month and day', hlm_current_month, hlm_current_day - else - !max condition 225 ft/min (FIREMIP Rabin table A10 JSBACH-Spitfire) convert to 68.577 m/min - wind_elev_fire = max(0.0_r8,(68.577_r8-0.5_r8*currentPatch%effect_wspeed)) - phi_wind = c * ((3.281_r8*wind_elev_fire)**b)*(beta_ratio**(-e)) - if (debug_windspeed) write(fates_log(),*) 'SF wind GREATER max ', currentPatch%effect_wspeed - if (debug_windspeed) write(fates_log(),*) 'month and day', hlm_current_month, hlm_current_day - endif + ! convert current_wspeed (wind at elev relevant to fire) from m/min to ft/min for Rothermel ROS eqn + phi_wind = c * ((3.281_r8*currentPatch%effect_wspeed)**b)*(beta_ratio**(-e)) + ! ---propagating flux---- ! Equation A2 in Thonicke et al.2010 and Eq. 42 Rothermal 1972 diff --git a/fire/SFParamsMod.F90 b/fire/SFParamsMod.F90 index a8556577fa..ecf4751f81 100644 --- a/fire/SFParamsMod.F90 +++ b/fire/SFParamsMod.F90 @@ -23,8 +23,7 @@ module SFParamsMod real(r8),protected :: SF_val_miner_damp real(r8),protected :: SF_val_max_durat real(r8),protected :: SF_val_durat_slope - real(r8),protected :: SF_val_wind_max ! Maximum wind speed expected by fire model (m/min) - real(r8),protected :: SF_val_alpha_FMC(NFSC) + real(r8),protected :: SF_val_drying_ratio real(r8),protected :: SF_val_CWD_frac(NCWD) real(r8),protected :: SF_val_max_decomp(NFSC) real(r8),protected :: SF_val_SAV(NFSC) @@ -45,7 +44,7 @@ module SFParamsMod character(len=param_string_length),parameter :: SF_name_miner_damp = "fates_miner_damp" character(len=param_string_length),parameter :: SF_name_max_durat = "fates_max_durat" character(len=param_string_length),parameter :: SF_name_durat_slope = "fates_durat_slope" - character(len=param_string_length),parameter :: SF_name_alpha_FMC = "fates_alpha_FMC" + character(len=param_string_length),parameter :: SF_name_drying_ratio = "fates_drying_ratio" character(len=param_string_length),parameter :: SF_name_CWD_frac = "fates_CWD_frac" character(len=param_string_length),parameter :: SF_name_max_decomp = "fates_max_decomp" character(len=param_string_length),parameter :: SF_name_SAV = "fates_SAV" @@ -56,7 +55,6 @@ module SFParamsMod character(len=param_string_length),parameter :: SF_name_low_moisture_Slope = "fates_low_moisture_Slope" character(len=param_string_length),parameter :: SF_name_mid_moisture_Coeff = "fates_mid_moisture_Coeff" character(len=param_string_length),parameter :: SF_name_mid_moisture_Slope = "fates_mid_moisture_Slope" - character(len=param_string_length),parameter :: SF_name_wind_max = "fates_fire_wind_max" public :: SpitFireRegisterParams public :: SpitFireReceiveParams @@ -90,11 +88,10 @@ subroutine SpitFireParamsInit() SF_val_miner_damp = nan SF_val_max_durat = nan SF_val_durat_slope = nan - SF_val_wind_max = nan + SF_val_drying_ratio = nan SF_val_CWD_frac(:) = nan - SF_val_alpha_FMC(:) = nan SF_val_max_decomp(:) = nan SF_val_SAV(:) = nan @@ -150,9 +147,6 @@ subroutine SpitFireRegisterScalars(fates_params) character(len=param_string_length), parameter :: dim_names_scalar(1) = (/dimension_name_scalar/) - call fates_params%RegisterParameter(name=SF_name_wind_max, dimension_shape=dimension_shape_scalar, & - dimension_names=dim_names_scalar) - call fates_params%RegisterParameter(name=SF_name_fdi_a, dimension_shape=dimension_shape_scalar, & dimension_names=dim_names_scalar) @@ -180,6 +174,9 @@ subroutine SpitFireRegisterScalars(fates_params) call fates_params%RegisterParameter(name=SF_name_durat_slope, dimension_shape=dimension_shape_scalar, & dimension_names=dim_names_scalar) + call fates_params%RegisterParameter(name=SF_name_drying_ratio, dimension_shape=dimension_shape_scalar, & + dimension_names=dim_names_scalar) + end subroutine SpitFireRegisterScalars !----------------------------------------------------------------------- @@ -191,9 +188,6 @@ subroutine SpitFireReceiveScalars(fates_params) class(fates_parameters_type), intent(inout) :: fates_params - call fates_params%RetreiveParameter(name=SF_name_wind_max, & - data=SF_val_wind_max) - call fates_params%RetreiveParameter(name=SF_name_fdi_a, & data=SF_val_fdi_a) @@ -221,6 +215,9 @@ subroutine SpitFireReceiveScalars(fates_params) call fates_params%RetreiveParameter(name=SF_name_durat_slope, & data=SF_val_durat_slope) + call fates_params%RetreiveParameter(name=SF_name_drying_ratio, & + data=SF_val_drying_ratio) + end subroutine SpitFireReceiveScalars !----------------------------------------------------------------------- @@ -288,9 +285,6 @@ subroutine SpitFireRegisterNFSC(fates_params) call fates_params%RegisterParameter(name=SF_name_mid_moisture_Slope, dimension_shape=dimension_shape_1d, & dimension_names=dim_names) - call fates_params%RegisterParameter(name=SF_name_alpha_FMC, dimension_shape=dimension_shape_1d, & - dimension_names=dim_names) - call fates_params%RegisterParameter(name=SF_name_max_decomp, dimension_shape=dimension_shape_1d, & dimension_names=dim_names) @@ -330,9 +324,6 @@ subroutine SpitFireReceiveNFSC(fates_params) call fates_params%RetreiveParameter(name=SF_name_mid_moisture_Slope, & data=SF_val_mid_moisture_Slope) - call fates_params%RetreiveParameter(name=SF_name_alpha_FMC, & - data=SF_val_alpha_FMC) - call fates_params%RetreiveParameter(name=SF_name_max_decomp, & data=SF_val_max_decomp) diff --git a/main/EDInitMod.F90 b/main/EDInitMod.F90 index 696ba70a72..50aa6d77c2 100644 --- a/main/EDInitMod.F90 +++ b/main/EDInitMod.F90 @@ -70,11 +70,17 @@ subroutine init_site_vars( site_in ) ! !LOCAL VARIABLES: !---------------------------------------------------------------------- ! - allocate(site_in%terminated_nindivs(1:nlevsclass,1:numpft,2)) + allocate(site_in%term_nindivs_canopy(1:nlevsclass,1:numpft)) + allocate(site_in%term_nindivs_ustory(1:nlevsclass,1:numpft)) allocate(site_in%demotion_rate(1:nlevsclass)) allocate(site_in%promotion_rate(1:nlevsclass)) allocate(site_in%imort_rate(1:nlevsclass,1:numpft)) + allocate(site_in%fmort_rate_canopy(1:nlevsclass,1:numpft)) + allocate(site_in%fmort_rate_ustory(1:nlevsclass,1:numpft)) + allocate(site_in%fmort_rate_cambial(1:nlevsclass,1:numpft)) + allocate(site_in%fmort_rate_crown(1:nlevsclass,1:numpft)) allocate(site_in%growthflux_fusion(1:nlevsclass,1:numpft)) + ! end subroutine init_site_vars @@ -123,11 +129,19 @@ subroutine zero_site( site_in ) site_in%fates_to_bgc_last_ts = 0.0_r8 ! termination and recruitment info - site_in%terminated_nindivs(:,:,:) = 0._r8 - site_in%termination_carbonflux(:) = 0._r8 + site_in%term_nindivs_canopy(:,:) = 0._r8 + site_in%term_nindivs_ustory(:,:) = 0._r8 + site_in%term_carbonflux_canopy = 0._r8 + site_in%term_carbonflux_ustory = 0._r8 site_in%recruitment_rate(:) = 0._r8 site_in%imort_rate(:,:) = 0._r8 site_in%imort_carbonflux = 0._r8 + site_in%fmort_rate_canopy(:,:) = 0._r8 + site_in%fmort_rate_ustory(:,:) = 0._r8 + site_in%fmort_carbonflux_canopy = 0._r8 + site_in%fmort_carbonflux_ustory = 0._r8 + site_in%fmort_rate_cambial(:,:) = 0._r8 + site_in%fmort_rate_crown(:,:) = 0._r8 ! fusoin-induced growth flux of individuals site_in%growthflux_fusion(:,:) = 0._r8 diff --git a/main/EDMainMod.F90 b/main/EDMainMod.F90 index 12aba1c7d3..5e9830a1dd 100644 --- a/main/EDMainMod.F90 +++ b/main/EDMainMod.F90 @@ -538,6 +538,12 @@ subroutine ed_total_balance_check (currentSite, call_index ) 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 + real(r8) :: leaf_c + real(r8) :: fnrt_c + real(r8) :: sapw_c + real(r8) :: store_c + real(r8) :: struct_c + real(r8) :: repro_c ! nb. There is no time associated with these variables ! because this routine can be called between any two @@ -579,7 +585,7 @@ subroutine ed_total_balance_check (currentSite, call_index ) ! burned_litter * new_patch%area !kG/site/day ! ----------------------------------------------------------------------------------- - if ( error_frac > 10e-6 ) then + if ( error_frac > 10e-6_r8 ) 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 @@ -593,13 +599,17 @@ subroutine ed_total_balance_check (currentSite, call_index ) write(fates_log(),*) 'seeds',seed_stock write(fates_log(),*) 'previous total',currentSite%old_stock - if(debug)then + 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 + 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 + currentPatch => currentSite%oldest_patch do while(associated(currentPatch)) write(fates_log(),*) '---------------------------------------' write(fates_log(),*) currentPatch%area , sum(currentPatch%cwd_ag), sum(currentPatch%cwd_bg) @@ -607,19 +617,22 @@ subroutine ed_total_balance_check (currentSite, call_index ) write(fates_log(),*)'---' currentCohort => currentPatch%tallest do while(associated(currentCohort)) - write(fates_log(),*) 'structure: ',currentCohort%prt%GetState(struct_organ,all_carbon_elements) - write(fates_log(),*) 'storage: ',currentCohort%prt%GetState(store_organ,all_carbon_elements) + write(fates_log(),*) 'pft: ',currentCohort%pft + write(fates_log(),*) 'dbh: ',currentCohort%dbh + leaf_c = currentCohort%prt%GetState(leaf_organ,all_carbon_elements) + struct_c = currentCohort%prt%GetState(struct_organ,all_carbon_elements) + store_c = currentCohort%prt%GetState(store_organ,all_carbon_elements) + fnrt_c = currentCohort%prt%GetState(fnrt_organ,all_carbon_elements) + repro_c = currentCohort%prt%GetState(repro_organ,all_carbon_elements) + sapw_c = currentCohort%prt%GetState(sapw_organ,all_carbon_elements) + + write(fates_log(),*) 'lc: ',leaf_c,' dc: ',struct_c,' stc: ',store_c + write(fates_log(),*) 'fc: ',fnrt_c,' rc: ',repro_c,' sac: ',sapw_c write(fates_log(),*) 'N plant: ',currentCohort%n - currentCohort => currentCohort%shorter; - enddo !end cohort loop + 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 + enddo !end patch loop write(fates_log(),*) 'aborting on date:',hlm_current_year,hlm_current_month,hlm_current_day call endrun(msg=errMsg(sourcefile, __LINE__)) end if @@ -673,7 +686,6 @@ subroutine bypass_dynamics(currentSite) currentCohort%bmort = 0.0_r8 currentCohort%hmort = 0.0_r8 currentCohort%cmort = 0.0_r8 - currentCohort%fmort = 0.0_r8 currentCohort%frmort = 0.0_r8 currentCohort%dndt = 0.0_r8 diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index f0d54dc90e..8f9457ed10 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -20,7 +20,9 @@ module EDTypesMod 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 + integer, parameter :: ican_ustory = 2 ! Nominal index for diagnostics that refer + ! to understory layers (all layers that + ! are not the top canopy layer) integer, parameter :: nlevleaf = 40 ! number of leaf layers in canopy layer integer, parameter :: maxpft = 15 ! maximum number of PFTs allowed @@ -144,11 +146,11 @@ module EDTypesMod integer , parameter :: N_HITE_BINS = 60 ! no. of hite bins used to distribute LAI ! COHORT TERMINATION - real(r8), parameter :: min_npm2 = 1.0E-8_r8 ! minimum cohort number density per m2 before termination - real(r8), parameter :: min_patch_area = 0.001_r8 ! smallest allowable patch area before termination - real(r8), parameter :: min_nppatch = 1.0E-11_r8 ! minimum number of cohorts per patch (min_npm2*min_patch_area) - real(r8), parameter :: min_n_safemath = 1.0E-15_r8 ! in some cases, we want to immediately remove super small - ! number densities of cohorts to prevent FPEs + real(r8), parameter :: min_npm2 = 1.0E-7_r8 ! minimum cohort number density per m2 before termination + real(r8), parameter :: min_patch_area = 0.01_r8 ! smallest allowable patch area before termination + real(r8), parameter :: min_nppatch = min_npm2*min_patch_area ! minimum number of cohorts per patch (min_npm2*min_patch_area) + real(r8), parameter :: min_n_safemath = 1.0E-12_r8 ! in some cases, we want to immediately remove super small + ! number densities of cohorts to prevent FPEs character*4 yearchar @@ -261,7 +263,6 @@ module EDTypesMod real(r8) :: bmort ! background mortality rate n/year real(r8) :: cmort ! carbon starvation mortality rate n/year real(r8) :: hmort ! hydraulic failure mortality rate n/year - real(r8) :: fmort ! fire mortality n/year real(r8) :: frmort ! freezing mortality n/year ! Logging Mortality Rate @@ -286,8 +287,8 @@ module EDTypesMod ! FIRE real(r8) :: fraction_crown_burned ! proportion of crown affected by fire:- - real(r8) :: cambial_mort ! probability that trees dies due to cambial char:- - real(r8) :: crownfire_mort ! probability of tree post-fire mortality due to crown scorch:- + real(r8) :: cambial_mort ! probability that trees dies due to cambial char (conditional on the tree being subjected to the fire) + real(r8) :: crownfire_mort ! probability of tree post-fire mortality due to crown scorch (conditional on the tree being subjected to the fire) real(r8) :: fire_mort ! post-fire mortality from cambial and crown damage assuming two are independent:- ! Hydraulics @@ -603,19 +604,41 @@ module EDTypesMod ! PLANT HYDRAULICS type(ed_site_hydr_type), pointer :: si_hydr - + + ! DIAGNOSTICS ! TERMINATION, RECRUITMENT, DEMOTION, and DISTURBANCE - - 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), 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), 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] + + real(r8), allocatable :: term_nindivs_canopy(:,:) ! number of canopy individuals that were in cohorts which + ! were terminated this timestep, on size x pft + real(r8), allocatable :: term_nindivs_ustory(:,:) ! number of understory individuals that were in cohorts which + ! were terminated this timestep, on size x pft + real(r8) :: term_carbonflux_canopy ! carbon flux from live to dead pools associated + ! with termination mortality, per canopy level + real(r8) :: term_carbonflux_ustory ! 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), 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), 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] real(r8), allocatable :: imort_rate(:,:) ! rate of individuals killed due to impact mortality per year. on size x pft array real(r8) :: imort_carbonflux ! biomass of individuals killed due to impact mortality per year. [kgC/ha/day] - real(r8), allocatable :: growthflux_fusion(:,:) ! rate of individuals moving into a given size class bin due to fusion in a given day. on size x pft array + + real(r8), allocatable :: fmort_rate_canopy(:,:) ! rate of canopy individuals killed due to fire mortality per year. + ! on size x pft array (1:nlevsclass,1:numpft) + real(r8), allocatable :: fmort_rate_ustory(:,:) ! rate of understory individuals killed due to fire mortality per year. + ! on size x pft array (1:nlevsclass,1:numpft) + real(r8) :: fmort_carbonflux_canopy ! biomass of canopy indivs killed due to fire per year. [gC/m2/sec] + real(r8) :: fmort_carbonflux_ustory ! biomass of understory indivs killed due to fire per year [gC/m2/sec] + real(r8), allocatable :: fmort_rate_cambial(:,:) ! rate of individuals killed due to fire mortality + ! from cambial damage per year. on size x pft array + real(r8), allocatable :: fmort_rate_crown(:,:) ! rate of individuals killed due to fire mortality + ! from crown damage per year. on size x pft array + + real(r8), allocatable :: growthflux_fusion(:,:) ! rate of individuals moving into a given size class bin + ! due to fusion in a given day. on size x pft array + ! some diagnostic-only (i.e. not resolved by ODE solver) flux of carbon to CWD and litter pools from termination and canopy mortality real(r8) :: CWD_AG_diagnostic_input_carbonflux(1:ncwd) ! diagnostic flux to AG CWD [kg C / m2 / yr] @@ -808,7 +831,6 @@ subroutine dump_cohort(ccohort) write(fates_log(),*) 'co%c_area = ', ccohort%c_area write(fates_log(),*) 'co%cmort = ', ccohort%cmort write(fates_log(),*) 'co%bmort = ', ccohort%bmort - write(fates_log(),*) 'co%fmort = ', ccohort%fmort write(fates_log(),*) 'co%hmort = ', ccohort%hmort write(fates_log(),*) 'co%frmort = ', ccohort%frmort write(fates_log(),*) 'co%isnew = ', ccohort%isnew diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index 760d180a51..9913aeca3f 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -8,7 +8,8 @@ module FatesHistoryInterfaceMod use FatesConstantsMod , only : calloc_abs_error use FatesGlobals , only : fates_log use FatesGlobals , only : endrun => fates_endrun - + use EDTypesMod , only : nclmax + use EDTypesMod , only : ican_upper use FatesIODimensionsMod , only : fates_io_dimension_type use FatesIOVariableKindMod , only : fates_io_variable_kind_type use FatesHistoryVariableType , only : fates_history_variable_type @@ -198,6 +199,8 @@ module FatesHistoryInterfaceMod integer, private :: ih_m6_si_scpf integer, private :: ih_m7_si_scpf integer, private :: ih_m8_si_scpf + integer, private :: ih_crownfiremort_si_scpf + integer, private :: ih_cambialfiremort_si_scpf integer, private :: ih_ar_si_scpf @@ -1431,6 +1434,8 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_m6_si_scpf => this%hvars(ih_m6_si_scpf)%r82d, & hio_m7_si_scpf => this%hvars(ih_m7_si_scpf)%r82d, & hio_m8_si_scpf => this%hvars(ih_m8_si_scpf)%r82d, & + hio_crownfiremort_si_scpf => this%hvars(ih_crownfiremort_si_scpf)%r82d, & + hio_cambialfiremort_si_scpf => this%hvars(ih_cambialfiremort_si_scpf)%r82d, & hio_m1_si_scls => this%hvars(ih_m1_si_scls)%r82d, & hio_m2_si_scls => this%hvars(ih_m2_si_scls)%r82d, & @@ -1779,7 +1784,6 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_m1_si_scpf(io_si,scpf) = hio_m1_si_scpf(io_si,scpf) + ccohort%bmort*ccohort%n hio_m2_si_scpf(io_si,scpf) = hio_m2_si_scpf(io_si,scpf) + ccohort%hmort*ccohort%n hio_m3_si_scpf(io_si,scpf) = hio_m3_si_scpf(io_si,scpf) + ccohort%cmort*ccohort%n - hio_m5_si_scpf(io_si,scpf) = hio_m5_si_scpf(io_si,scpf) + ccohort%fmort*ccohort%n hio_m7_si_scpf(io_si,scpf) = hio_m7_si_scpf(io_si,scpf) + & (ccohort%lmort_direct+ccohort%lmort_collateral+ccohort%lmort_infra) * ccohort%n hio_m8_si_scpf(io_si,scpf) = hio_m8_si_scpf(io_si,scpf) + ccohort%frmort*ccohort%n @@ -1787,7 +1791,6 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_m1_si_scls(io_si,scls) = hio_m1_si_scls(io_si,scls) + ccohort%bmort*ccohort%n hio_m2_si_scls(io_si,scls) = hio_m2_si_scls(io_si,scls) + ccohort%hmort*ccohort%n hio_m3_si_scls(io_si,scls) = hio_m3_si_scls(io_si,scls) + ccohort%cmort*ccohort%n - hio_m5_si_scls(io_si,scls) = hio_m5_si_scls(io_si,scls) + ccohort%fmort*ccohort%n hio_m7_si_scls(io_si,scls) = hio_m7_si_scls(io_si,scls) + & (ccohort%lmort_direct+ccohort%lmort_collateral+ccohort%lmort_infra) * ccohort%n hio_m8_si_scls(io_si,scls) = hio_m8_si_scls(io_si,scls) + & @@ -1845,7 +1848,7 @@ subroutine update_history_dyn(this,nc,nsites,sites) if (ccohort%canopy_layer .eq. 1) then hio_nplant_canopy_si_scag(io_si,iscag) = hio_nplant_canopy_si_scag(io_si,iscag) + ccohort%n hio_mortality_canopy_si_scag(io_si,iscag) = hio_mortality_canopy_si_scag(io_si,iscag) + & - (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%fmort + ccohort%frmort) * ccohort%n + (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%frmort) * ccohort%n hio_ddbh_canopy_si_scag(io_si,iscag) = hio_ddbh_canopy_si_scag(io_si,iscag) + & ccohort%ddbhdt*ccohort%n hio_bstor_canopy_si_scpf(io_si,scpf) = hio_bstor_canopy_si_scpf(io_si,scpf) + & @@ -1856,10 +1859,10 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_canopy_biomass_pa(io_pa) = hio_canopy_biomass_pa(io_pa) + n_density * total_c * 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%frmort) * ccohort%n + ! (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%frmort) * 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%fmort + ccohort%frmort) * ccohort%n + & + (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%frmort) * ccohort%n + & (ccohort%lmort_direct + ccohort%lmort_collateral + ccohort%lmort_infra) * & ccohort%n * sec_per_day * days_per_year @@ -1885,12 +1888,12 @@ subroutine update_history_dyn(this,nc,nsites,sites) ! 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%frmort ) * ccohort%n + & + (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%frmort ) * ccohort%n + & (ccohort%lmort_direct + ccohort%lmort_collateral + ccohort%lmort_infra) * & ccohort%n * sec_per_day * days_per_year hio_canopy_mortality_carbonflux_si(io_si) = hio_canopy_mortality_carbonflux_si(io_si) + & - (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%fmort + ccohort%frmort) * & + (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%frmort) * & total_c * ccohort%n * g_per_kg * days_per_sec * years_per_day * ha_per_m2 + & (ccohort%lmort_direct + ccohort%lmort_collateral + ccohort%lmort_infra) * total_c * & ccohort%n * g_per_kg * ha_per_m2 @@ -1932,7 +1935,7 @@ subroutine update_history_dyn(this,nc,nsites,sites) else hio_nplant_understory_si_scag(io_si,iscag) = hio_nplant_understory_si_scag(io_si,iscag) + ccohort%n hio_mortality_understory_si_scag(io_si,iscag) = hio_mortality_understory_si_scag(io_si,iscag) + & - (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%fmort + ccohort%frmort) * ccohort%n + (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%frmort) * ccohort%n hio_ddbh_understory_si_scag(io_si,iscag) = hio_ddbh_understory_si_scag(io_si,iscag) + & ccohort%ddbhdt*ccohort%n hio_bstor_understory_si_scpf(io_si,scpf) = hio_bstor_understory_si_scpf(io_si,scpf) + & @@ -1941,9 +1944,11 @@ subroutine update_history_dyn(this,nc,nsites,sites) leaf_c * ccohort%n hio_understory_biomass_pa(io_pa) = hio_understory_biomass_pa(io_pa) + & n_density * total_c * 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%frmort) * 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%frmort ) * ccohort%n + & + (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%frmort ) * ccohort%n + & (ccohort%lmort_direct + ccohort%lmort_collateral + ccohort%lmort_infra) * & ccohort%n * sec_per_day * days_per_year @@ -1970,12 +1975,12 @@ subroutine update_history_dyn(this,nc,nsites,sites) ! 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%frmort ) * ccohort%n + & + (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%frmort ) * ccohort%n + & (ccohort%lmort_direct + ccohort%lmort_collateral + ccohort%lmort_infra) * & ccohort%n * sec_per_day * days_per_year hio_understory_mortality_carbonflux_si(io_si) = hio_understory_mortality_carbonflux_si(io_si) + & - (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%fmort + ccohort%frmort) * & + (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%frmort) * & total_c * ccohort%n * g_per_kg * days_per_sec * years_per_day * ha_per_m2 + & (ccohort%lmort_direct + ccohort%lmort_collateral + ccohort%lmort_infra) * total_c * & ccohort%n * g_per_kg * ha_per_m2 @@ -2137,21 +2142,28 @@ subroutine update_history_dyn(this,nc,nsites,sites) i_scpf = (i_pft-1)*nlevsclass + i_scls ! ! termination mortality. sum of canopy and understory indices - 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_m6_si_scpf(io_si,i_scpf) = (sites(s)%term_nindivs_canopy(i_scls,i_pft) + & + sites(s)%term_nindivs_ustory(i_scls,i_pft)) * days_per_year + hio_m6_si_scls(io_si,i_scls) = hio_m6_si_scls(io_si,i_scls) + & - (sites(s)%terminated_nindivs(i_scls,i_pft,1) + & - sites(s)%terminated_nindivs(i_scls,i_pft,2)) * days_per_year + (sites(s)%term_nindivs_canopy(i_scls,i_pft) + & + sites(s)%term_nindivs_ustory(i_scls,i_pft)) * days_per_year + + ! ! add termination mortality to canopy and understory mortality hio_mortality_canopy_si_scls(io_si,i_scls) = hio_mortality_canopy_si_scls(io_si,i_scls) + & - sites(s)%terminated_nindivs(i_scls,i_pft,1) * days_per_year + sites(s)%term_nindivs_canopy(i_scls,i_pft) * days_per_year + hio_mortality_understory_si_scls(io_si,i_scls) = hio_mortality_understory_si_scls(io_si,i_scls) + & - sites(s)%terminated_nindivs(i_scls,i_pft,2) * days_per_year + sites(s)%term_nindivs_ustory(i_scls,i_pft) * days_per_year + hio_mortality_canopy_si_scpf(io_si,i_scpf) = hio_mortality_canopy_si_scpf(io_si,i_scpf) + & - sites(s)%terminated_nindivs(i_scls,i_pft,1) * days_per_year + sites(s)%term_nindivs_canopy(i_scls,i_pft) * days_per_year + hio_mortality_understory_si_scpf(io_si,i_scpf) = hio_mortality_understory_si_scpf(io_si,i_scpf) + & - sites(s)%terminated_nindivs(i_scls,i_pft,2) * days_per_year + sites(s)%term_nindivs_ustory(i_scls,i_pft) * days_per_year + ! ! imort on its own hio_m4_si_scpf(io_si,i_scpf) = sites(s)%imort_rate(i_scls, i_pft) @@ -2160,7 +2172,7 @@ subroutine update_history_dyn(this,nc,nsites,sites) ! add imort to other mortality terms. consider imort as understory mortality even if it happens in ! cohorts that may have been promoted as part of the patch creation, and use the pre-calculated site-level ! values to avoid biasing the results by the dramatically-reduced number densities in cohorts that are subject to imort - hio_mortality_understory_si_scpf(io_si,i_scpf) = hio_mortality_understory_si_scpf(io_si,i_scpf)+ & + hio_mortality_understory_si_scpf(io_si,i_scpf) = hio_mortality_understory_si_scpf(io_si,i_scpf) + & sites(s)%imort_rate(i_scls, i_pft) hio_mortality_understory_si_scls(io_si,i_scls) = hio_mortality_understory_si_scls(io_si,i_scls) + & sites(s)%imort_rate(i_scls, i_pft) @@ -2168,10 +2180,50 @@ subroutine update_history_dyn(this,nc,nsites,sites) iscag = i_scls ! since imort is by definition something that only happens in newly disturbed patches, treat as such hio_mortality_understory_si_scag(io_si,iscag) = hio_mortality_understory_si_scag(io_si,iscag) + & sites(s)%imort_rate(i_scls, i_pft) + + ! fire mortality from the site-level diagnostic rates + hio_m5_si_scpf(io_si,i_scpf) = sites(s)%fmort_rate_canopy(i_scls, i_pft) + & + sites(s)%fmort_rate_ustory(i_scls, i_pft) + hio_m5_si_scls(io_si,i_scls) = hio_m5_si_scls(io_si,i_scls) + & + sites(s)%fmort_rate_canopy(i_scls, i_pft) + sites(s)%fmort_rate_ustory(i_scls, i_pft) + ! + hio_crownfiremort_si_scpf(io_si,i_scpf) = sites(s)%fmort_rate_crown(i_scls, i_pft) + hio_cambialfiremort_si_scpf(io_si,i_scpf) = sites(s)%fmort_rate_cambial(i_scls, i_pft) + ! + ! fire components of overall canopy and understory mortality + hio_mortality_canopy_si_scpf(io_si,i_scpf) = hio_mortality_canopy_si_scpf(io_si,i_scpf) + & + sites(s)%fmort_rate_canopy(i_scls, i_pft) + hio_mortality_canopy_si_scls(io_si,i_scls) = hio_mortality_canopy_si_scls(io_si,i_scls) + & + sites(s)%fmort_rate_canopy(i_scls, i_pft) + + ! the fire mortality rates for each layer are total dead, since the usable + ! output will then normalize by the counts, we are allowed to sum over layers + hio_mortality_understory_si_scpf(io_si,i_scpf) = hio_mortality_understory_si_scpf(io_si,i_scpf) + & + sites(s)%fmort_rate_ustory(i_scls, i_pft) + + hio_mortality_understory_si_scls(io_si,i_scls) = hio_mortality_understory_si_scls(io_si,i_scls) + & + sites(s)%fmort_rate_ustory(i_scls, i_pft) + + ! + ! carbon flux associated with mortality of trees dying by fire + hio_canopy_mortality_carbonflux_si(io_si) = hio_canopy_mortality_carbonflux_si(io_si) + & + sites(s)%fmort_carbonflux_canopy + + hio_understory_mortality_carbonflux_si(io_si) = hio_understory_mortality_carbonflux_si(io_si) + & + sites(s)%fmort_carbonflux_ustory + ! + ! for scag variables, also treat as happening in the newly-disurbed patch + + hio_mortality_canopy_si_scag(io_si,iscag) = hio_mortality_canopy_si_scag(io_si,iscag) + & + sites(s)%fmort_rate_canopy(i_scls, i_pft) + hio_mortality_understory_si_scag(io_si,iscag) = hio_mortality_understory_si_scag(io_si,iscag) + & + sites(s)%fmort_rate_ustory(i_scls, i_pft) + ! while in this loop, pass the fusion-induced growth rate flux to history hio_growthflux_fusion_si_scpf(io_si,i_scpf) = hio_growthflux_fusion_si_scpf(io_si,i_scpf) + & sites(s)%growthflux_fusion(i_scls, i_pft) * days_per_year + end do end do ! @@ -2179,10 +2231,16 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_understory_mortality_carbonflux_si(io_si) = hio_understory_mortality_carbonflux_si(io_si) + & sites(s)%imort_carbonflux ! - sites(s)%terminated_nindivs(:,:,:) = 0._r8 + sites(s)%term_nindivs_canopy(:,:) = 0._r8 + sites(s)%term_nindivs_ustory(:,:) = 0._r8 sites(s)%imort_carbonflux = 0._r8 sites(s)%imort_rate(:,:) = 0._r8 - ! + sites(s)%fmort_rate_canopy(:,:) = 0._r8 + sites(s)%fmort_rate_ustory(:,:) = 0._r8 + sites(s)%fmort_carbonflux_canopy = 0._r8 + sites(s)%fmort_carbonflux_ustory = 0._r8 + sites(s)%fmort_rate_cambial(:,:) = 0._r8 + sites(s)%fmort_rate_crown(:,:) = 0._r8 sites(s)%growthflux_fusion(:,:) = 0._r8 ! pass the recruitment rate as a flux to the history, and then reset the recruitment buffer @@ -2222,11 +2280,14 @@ subroutine update_history_dyn(this,nc,nsites,sites) ! mortality-associated carbon fluxes hio_canopy_mortality_carbonflux_si(io_si) = hio_canopy_mortality_carbonflux_si(io_si) + & - sites(s)%termination_carbonflux(ican_upper) * g_per_kg * days_per_sec * ha_per_m2 + sites(s)%term_carbonflux_canopy * g_per_kg * days_per_sec * ha_per_m2 + hio_understory_mortality_carbonflux_si(io_si) = hio_understory_mortality_carbonflux_si(io_si) + & - sites(s)%termination_carbonflux(ican_ustory) * g_per_kg * days_per_sec * ha_per_m2 + sites(s)%term_carbonflux_ustory * g_per_kg * days_per_sec * ha_per_m2 + ! and zero the site-level termination carbon flux variable - sites(s)%termination_carbonflux(:) = 0._r8 + sites(s)%term_carbonflux_canopy = 0._r8 + sites(s)%term_carbonflux_ustory = 0._r8 ! ! add the site-level disturbance-associated cwd and litter input fluxes to thir respective flux fields do i_cwd = 1, ncwd @@ -2551,8 +2612,8 @@ subroutine update_history_prod(this,nc,nsites,sites,dt_tstep) ! summarize radiation profiles through the canopy do ipft=1,numpft - do ican=1,nclmax - do ileaf=1,nlevleaf + do ican=1,nclmax ! cpatch%ncl_p ? + do ileaf=1,nlevleaf ! cpatch%ncan(ican,ipft) ? ! calculate where we are on multiplexed dimensions cnlfpft_indx = ileaf + (ican-1) * nlevleaf + (ipft-1) * nlevleaf * nclmax cnlf_indx = ileaf + (ican-1) * nlevleaf @@ -3882,6 +3943,16 @@ 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_m5_si_scpf ) + call this%set_history_var(vname='CROWNFIREMORT_SCPF', units = 'N/ha/yr', & + long='crown fire mortality 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_crownfiremort_si_scpf ) + + call this%set_history_var(vname='CAMBIALFIREMORT_SCPF', units = 'N/ha/yr', & + long='cambial fire mortality 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_cambialfiremort_si_scpf ) + call this%set_history_var(vname='M6_SCPF', units = 'N/ha/yr', & long='termination mortality by pft/size',use_default='inactive', & avgflag='A', vtype=site_size_pft_r8, hlms='CLM:ALM', flushval=0.0_r8, & diff --git a/main/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index c108be8757..fb9b316f45 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -1,14 +1,14 @@ module FatesRestartInterfaceMod - use FatesConstantsMod , only : r8 => fates_r8 - use FatesConstantsMod , only : fates_avg_flag_length - use FatesConstantsMod , only : fates_short_string_length - use FatesConstantsMod , only : fates_long_string_length - use FatesConstantsMod , only : itrue - use FatesConstantsMod , only : ifalse - use FatesGlobals , only : fates_log - use FatesGlobals , only : endrun => fates_endrun + use FatesConstantsMod, only : r8 => fates_r8 + use FatesConstantsMod, only : fates_avg_flag_length + use FatesConstantsMod, only : fates_short_string_length + use FatesConstantsMod, only : fates_long_string_length + use FatesConstantsMod, only : itrue + use FatesConstantsMod, only : ifalse + use FatesGlobals, only : fates_log + use FatesGlobals, only : endrun => fates_endrun use FatesIODimensionsMod, only : fates_io_dimension_type use FatesIOVariableKindMod, only : fates_io_variable_kind_type use FatesRestartVariableMod, only : fates_restart_variable_type @@ -26,6 +26,9 @@ module FatesRestartInterfaceMod use EDCohortDynamicsMod, only : zero_cohort use EDCohortDynamicsMod, only : InitPRTCohort use FatesPlantHydraulicsMod, only : InitHydrCohort + use FatesInterfaceMod, only : nlevsclass + use PRTGenericMod, only : prt_global + ! CIME GLOBALS use shr_log_mod , only : errMsg => shr_log_errMsg @@ -105,7 +108,6 @@ module FatesRestartInterfaceMod integer, private :: ir_bmort_co integer, private :: ir_hmort_co integer, private :: ir_cmort_co - integer, private :: ir_fmort_co integer, private :: ir_frmort_co !Logging @@ -132,13 +134,35 @@ module FatesRestartInterfaceMod integer, private :: ir_root_litter_paft integer, private :: ir_leaf_litter_in_paft integer, private :: ir_root_litter_in_paft - integer, private :: ir_seed_bank_sift - integer, private :: ir_spread_si + integer, private :: ir_livegrass_pa integer, private :: ir_age_pa integer, private :: ir_area_pa + ! Site level integer, private :: ir_watermem_siwm + integer, private :: ir_seed_bank_sift + integer, private :: ir_spread_si + integer, private :: ir_recrate_sift + integer, private :: ir_fmortrate_cano_siscpf + integer, private :: ir_fmortrate_usto_siscpf + integer, private :: ir_imortrate_siscpf + integer, private :: ir_fmortrate_crown_siscpf + integer, private :: ir_fmortrate_cambi_siscpf + integer, private :: ir_termnindiv_cano_siscpf + integer, private :: ir_termnindiv_usto_siscpf + integer, private :: ir_growflx_fusion_siscpf + integer, private :: ir_demorate_sisc + integer, private :: ir_promrate_sisc + integer, private :: ir_termcflux_cano_si + integer, private :: ir_termcflux_usto_si + integer, private :: ir_democflux_si + integer, private :: ir_promcflux_si + integer, private :: ir_imortcflux_si + integer, private :: ir_fmortcflux_cano_si + integer, private :: ir_fmortcflux_usto_si + + integer, private :: ir_prt_base ! Base index for all PRT variables @@ -732,11 +756,6 @@ subroutine define_restart_vars(this, initialize_variables) units='/year', flushval = flushzero, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_cmort_co ) - call this%set_restart_var(vname='fates_fmort', vtype=cohort_r8, & - long_name='ed cohort - fire mortality rate', & - units='/year', flushval = flushzero, & - hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_fmort_co ) - call this%set_restart_var(vname='fates_frmort', vtype=cohort_r8, & long_name='ed cohort - freezing mortality rate', & units='/year', flushval = flushzero, & @@ -942,6 +961,99 @@ subroutine define_restart_vars(this, initialize_variables) units='m3/m3', flushval = flushzero, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_watermem_siwm ) + + call this%set_restart_var(vname='fates_recrate', vtype=cohort_r8, & + long_name='fates diagnostics on recruitment', & + units='indiv/ha/day', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_recrate_sift) + + call this%set_restart_var(vname='fates_fmortrate_canopy', vtype=cohort_r8, & + long_name='fates diagnostics on fire mortality canopy', & + units='indiv/ha/year', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_fmortrate_cano_siscpf) + + call this%set_restart_var(vname='fates_fmortrate_ustory', vtype=cohort_r8, & + long_name='fates diagnostics on fire mortality ustory', & + units='indiv/ha/year', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_fmortrate_usto_siscpf) + + call this%set_restart_var(vname='fates_imortrate', vtype=cohort_r8, & + long_name='fates diagnostics on impact mortality', & + units='indiv/ha/year', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_imortrate_siscpf) + + call this%set_restart_var(vname='fates_fmortrate_crown', vtype=cohort_r8, & + long_name='fates diagnostics on crown fire mortality', & + units='indiv/ha/year', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_fmortrate_crown_siscpf) + + call this%set_restart_var(vname='fates_fmortrate_cambi', vtype=cohort_r8, & + long_name='fates diagnostics on fire cambial mortality', & + units='indiv/ha/year', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_fmortrate_cambi_siscpf) + + call this%set_restart_var(vname='fates_termn_canopy', vtype=cohort_r8, & + long_name='fates diagnostics on termin mortality canopy', & + units='indiv/ha/day', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_termnindiv_cano_siscpf) + + call this%set_restart_var(vname='fates_termn_ustory', vtype=cohort_r8, & + long_name='fates diagnostics on term mortality ustory', & + units='indiv/ha/day', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_termnindiv_usto_siscpf) + + call this%set_restart_var(vname='fates_growflx_fusion', vtype=cohort_r8, & + long_name='fates diag: rate of indivs moving via fusion', & + units='indiv/ha/day', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_growflx_fusion_siscpf) + + call this%set_restart_var(vname='fates_demorate', vtype=cohort_r8, & + long_name='fates diagnoatic rate of indivs demoted', & + units='indiv/ha/day', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_demorate_sisc) + + call this%set_restart_var(vname='fates_promrate', vtype=cohort_r8, & + long_name='fates diagnostic rate of indivs promoted', & + units='indiv/ha/da', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_promrate_sisc) + + call this%set_restart_var(vname='fates_imortcflux', vtype=site_r8, & + long_name='biomass of indivs killed due to impact mort', & + units='kgC/ha/day', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_imortcflux_si) + + call this%set_restart_var(vname='fates_fmortcflux_canopy', vtype=site_r8, & + long_name='fates diagnostic biomass of canopy fire', & + units='gC/m2/sec', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_fmortcflux_cano_si) + + call this%set_restart_var(vname='fates_fmortcflux_ustory', vtype=site_r8, & + long_name='fates diagnostic biomass of understory fire', & + units='gC/m2/sec', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_fmortcflux_usto_si) + + call this%set_restart_var(vname='fates_termcflux_canopy', vtype=site_r8, & + long_name='fates diagnostic term carbon flux canopy', & + units='', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_termcflux_cano_si ) + + call this%set_restart_var(vname='fates_termcflux_ustory', vtype=site_r8, & + long_name='fates diagnostic term carbon flux understory', & + units='', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_termcflux_usto_si ) + + call this%set_restart_var(vname='fates_democflux', vtype=site_r8, & + long_name='fates diagnostic demotion carbon flux', & + units='', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_democflux_si ) + + call this%set_restart_var(vname='fates_promcflux', vtype=site_r8, & + long_name='fates diagnostic promotion carbon flux ', & + units='', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_promcflux_si ) + + + ! Register all of the PRT states and fluxes @@ -1288,6 +1400,8 @@ subroutine set_restart_vectors(this,nc,nsites,sites) integer :: io_idx_pa_ib ! each SW band (vis/ir) per patch (pa_ib) integer :: io_idx_si_wmem ! each water memory class within each site integer :: io_idx_si_lyr_shell ! site - layer x shell index + integer :: io_idx_si_scpf ! each size-class x pft index within site + integer :: io_idx_si_sc ! each size-class index within site ! Some counters (for checking mostly) integer :: totalcohorts ! total cohort count on this thread (diagnostic) @@ -1299,6 +1413,8 @@ subroutine set_restart_vectors(this,nc,nsites,sites) integer :: ir_prt_var ! loop counter for var x position integer :: i_var ! loop counter for PRT variables integer :: i_pos ! loop counter for discrete PRT positions + integer :: i_scls ! loop counter for size-class + integer :: i_pft ! loop counter for pft type(fates_restart_variable_type) :: rvar type(ed_patch_type),pointer :: cpatch @@ -1349,7 +1465,6 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_bmort_co => this%rvars(ir_bmort_co)%r81d, & rio_hmort_co => this%rvars(ir_hmort_co)%r81d, & rio_cmort_co => this%rvars(ir_cmort_co)%r81d, & - rio_fmort_co => this%rvars(ir_fmort_co)%r81d, & rio_frmort_co => this%rvars(ir_frmort_co)%r81d, & rio_lmort_direct_co => this%rvars(ir_lmort_direct_co)%r81d, & rio_lmort_collateral_co => this%rvars(ir_lmort_collateral_co)%r81d, & @@ -1372,7 +1487,25 @@ subroutine set_restart_vectors(this,nc,nsites,sites) 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, & - rio_watermem_siwm => this%rvars(ir_watermem_siwm)%r81d) + rio_watermem_siwm => this%rvars(ir_watermem_siwm)%r81d, & + rio_recrate_sift => this%rvars(ir_recrate_sift)%r81d, & + rio_fmortrate_cano_siscpf => this%rvars(ir_fmortrate_cano_siscpf)%r81d, & + rio_fmortrate_usto_siscpf => this%rvars(ir_fmortrate_usto_siscpf)%r81d, & + rio_imortrate_siscpf => this%rvars(ir_imortrate_siscpf)%r81d, & + rio_fmortrate_crown_siscpf => this%rvars(ir_fmortrate_crown_siscpf)%r81d, & + rio_fmortrate_cambi_siscpf => this%rvars(ir_fmortrate_cambi_siscpf)%r81d, & + rio_termnindiv_cano_siscpf => this%rvars(ir_termnindiv_cano_siscpf)%r81d, & + rio_termnindiv_usto_siscpf => this%rvars(ir_termnindiv_usto_siscpf)%r81d, & + rio_growflx_fusion_siscpf => this%rvars(ir_growflx_fusion_siscpf)%r81d, & + rio_demorate_sisc => this%rvars(ir_demorate_sisc)%r81d, & + rio_promrate_sisc => this%rvars(ir_promrate_sisc)%r81d, & + rio_termcflux_cano_si => this%rvars(ir_termcflux_cano_si)%r81d, & + rio_termcflux_usto_si => this%rvars(ir_termcflux_usto_si)%r81d, & + rio_democflux_si => this%rvars(ir_democflux_si)%r81d, & + rio_promcflux_si => this%rvars(ir_promcflux_si)%r81d, & + rio_imortcflux_si => this%rvars(ir_imortcflux_si)%r81d, & + rio_fmortcflux_cano_si => this%rvars(ir_fmortcflux_cano_si)%r81d, & + rio_fmortcflux_usto_si => this%rvars(ir_fmortcflux_usto_si)%r81d) totalCohorts = 0 @@ -1399,12 +1532,17 @@ subroutine set_restart_vectors(this,nc,nsites,sites) io_idx_pa_ib = io_idx_co_1st io_idx_si_wmem = io_idx_co_1st + ! Hydraulics counters lyr = hydraulic layer, shell = rhizosphere shell io_idx_si_lyr_shell = io_idx_co_1st + io_idx_si_scpf = io_idx_co_1st + io_idx_si_sc = io_idx_co_1st + ! write seed_bank info(site-level, but PFT-resolved) - do i = 1,numpft - rio_seed_bank_sift(io_idx_co_1st+i-1) = sites(s)%seed_bank(i) + do i_pft = 1,numpft + rio_seed_bank_sift(io_idx_co_1st+i_pft-1) = sites(s)%seed_bank(i_pft) + rio_recrate_sift(io_idx_co_1st+i_pft-1) = sites(s)%recruitment_rate(i_pft) end do ! canopy spread term @@ -1511,7 +1649,6 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_bmort_co(io_idx_co) = ccohort%bmort rio_hmort_co(io_idx_co) = ccohort%hmort rio_cmort_co(io_idx_co) = ccohort%cmort - rio_fmort_co(io_idx_co) = ccohort%fmort rio_frmort_co(io_idx_co) = ccohort%frmort !Logging @@ -1607,6 +1744,41 @@ subroutine set_restart_vectors(this,nc,nsites,sites) enddo ! cpatch do while + + ! Fill the site level diagnostics arrays + do i_scls = 1, nlevsclass + + do i_pft = 1, numpft + + rio_fmortrate_cano_siscpf(io_idx_si_scpf) = sites(s)%fmort_rate_canopy(i_scls, i_pft) + rio_fmortrate_usto_siscpf(io_idx_si_scpf) = sites(s)%fmort_rate_ustory(i_scls, i_pft) + rio_imortrate_siscpf(io_idx_si_scpf) = sites(s)%imort_rate(i_scls, i_pft) + rio_fmortrate_crown_siscpf(io_idx_si_scpf) = sites(s)%fmort_rate_crown(i_scls, i_pft) + rio_fmortrate_cambi_siscpf(io_idx_si_scpf) = sites(s)%fmort_rate_cambial(i_scls, i_pft) + rio_termnindiv_cano_siscpf(io_idx_si_scpf) = sites(s)%term_nindivs_canopy(i_scls,i_pft) + rio_termnindiv_usto_siscpf(io_idx_si_scpf) = sites(s)%term_nindivs_ustory(i_scls,i_pft) + rio_growflx_fusion_siscpf(io_idx_si_scpf) = sites(s)%growthflux_fusion(i_scls, i_pft) + + io_idx_si_scpf = io_idx_si_scpf + 1 + end do + + rio_demorate_sisc(io_idx_si_sc) = sites(s)%demotion_rate(i_scls) + rio_promrate_sisc(io_idx_si_sc) = sites(s)%promotion_rate(i_scls) + + io_idx_si_sc = io_idx_si_sc + 1 + end do + + + rio_termcflux_cano_si(io_idx_si) = sites(s)%term_carbonflux_canopy + rio_termcflux_usto_si(io_idx_si) = sites(s)%term_carbonflux_ustory + rio_democflux_si(io_idx_si) = sites(s)%demotion_carbonflux + rio_promcflux_si(io_idx_si) = sites(s)%promotion_carbonflux + rio_imortcflux_si(io_idx_si) = sites(s)%imort_carbonflux + rio_fmortcflux_cano_si(io_idx_si) = sites(s)%fmort_carbonflux_canopy + rio_fmortcflux_usto_si(io_idx_si) = sites(s)%fmort_carbonflux_ustory + + + rio_old_stock_si(io_idx_si) = sites(s)%old_stock rio_cd_status_si(io_idx_si) = sites(s)%status rio_dd_status_si(io_idx_si) = sites(s)%dstatus @@ -1932,6 +2104,8 @@ subroutine get_restart_vectors(this, nc, nsites, sites) integer :: io_idx_pa_ib ! each SW radiation band per patch (pa_ib) integer :: io_idx_si_wmem ! each water memory class within each site integer :: io_idx_si_lyr_shell ! site - layer x shell index + integer :: io_idx_si_scpf ! each size-class x pft index within site + integer :: io_idx_si_sc ! each size-class index within site ! Some counters (for checking mostly) integer :: totalcohorts ! total cohort count on this thread (diagnostic) @@ -1940,7 +2114,8 @@ subroutine get_restart_vectors(this, nc, nsites, sites) integer :: ir_prt_var ! loop counter for var x position integer :: i_var ! loop counter for PRT variables integer :: i_pos ! loop counter for discrete PRT positions - + integer :: i_pft ! loop counter for pft + integer :: i_scls ! loop counter for size-class associate( rio_npatch_si => this%rvars(ir_npatch_si)%int1d, & rio_old_stock_si => this%rvars(ir_oldstock_si)%r81d, & @@ -1986,7 +2161,6 @@ subroutine get_restart_vectors(this, nc, nsites, sites) rio_bmort_co => this%rvars(ir_bmort_co)%r81d, & rio_hmort_co => this%rvars(ir_hmort_co)%r81d, & rio_cmort_co => this%rvars(ir_cmort_co)%r81d, & - rio_fmort_co => this%rvars(ir_fmort_co)%r81d, & rio_frmort_co => this%rvars(ir_frmort_co)%r81d, & rio_lmort_direct_co => this%rvars(ir_lmort_direct_co)%r81d, & rio_lmort_collateral_co => this%rvars(ir_lmort_collateral_co)%r81d, & @@ -2009,8 +2183,27 @@ subroutine get_restart_vectors(this, nc, nsites, sites) 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, & - rio_watermem_siwm => this%rvars(ir_watermem_siwm)%r81d) + rio_watermem_siwm => this%rvars(ir_watermem_siwm)%r81d, & + rio_recrate_sift => this%rvars(ir_recrate_sift)%r81d, & + rio_fmortrate_cano_siscpf => this%rvars(ir_fmortrate_cano_siscpf)%r81d, & + rio_fmortrate_usto_siscpf => this%rvars(ir_fmortrate_usto_siscpf)%r81d, & + rio_imortrate_siscpf => this%rvars(ir_imortrate_siscpf)%r81d, & + rio_fmortrate_crown_siscpf => this%rvars(ir_fmortrate_crown_siscpf)%r81d, & + rio_fmortrate_cambi_siscpf => this%rvars(ir_fmortrate_cambi_siscpf)%r81d, & + rio_termnindiv_cano_siscpf => this%rvars(ir_termnindiv_cano_siscpf)%r81d, & + rio_termnindiv_usto_siscpf => this%rvars(ir_termnindiv_usto_siscpf)%r81d, & + rio_growflx_fusion_siscpf => this%rvars(ir_growflx_fusion_siscpf)%r81d, & + rio_demorate_sisc => this%rvars(ir_demorate_sisc)%r81d, & + rio_promrate_sisc => this%rvars(ir_promrate_sisc)%r81d, & + rio_termcflux_cano_si => this%rvars(ir_termcflux_cano_si)%r81d, & + rio_termcflux_usto_si => this%rvars(ir_termcflux_usto_si)%r81d, & + rio_democflux_si => this%rvars(ir_democflux_si)%r81d, & + rio_promcflux_si => this%rvars(ir_promcflux_si)%r81d, & + rio_imortcflux_si => this%rvars(ir_imortcflux_si)%r81d, & + rio_fmortcflux_cano_si => this%rvars(ir_fmortcflux_cano_si)%r81d, & + rio_fmortcflux_usto_si => this%rvars(ir_fmortcflux_usto_si)%r81d) + totalcohorts = 0 do s = 1,nsites @@ -2027,10 +2220,13 @@ subroutine get_restart_vectors(this, nc, nsites, sites) ! Hydraulics counters lyr = hydraulic layer, shell = rhizosphere shell io_idx_si_lyr_shell = io_idx_co_1st + io_idx_si_scpf = io_idx_co_1st + io_idx_si_sc = io_idx_co_1st ! read seed_bank info(site-level, but PFT-resolved) - do i = 1,numpft - sites(s)%seed_bank(i) = rio_seed_bank_sift(io_idx_co_1st+i-1) + do i_pft = 1,numpft + sites(s)%seed_bank(i_pft) = rio_seed_bank_sift(io_idx_co_1st+i_pft-1) + sites(s)%recruitment_rate(i_pft) = rio_recrate_sift(io_idx_co_1st+i_pft-1) enddo sites(s)%spread = rio_spread_si(io_idx_si) @@ -2106,7 +2302,6 @@ subroutine get_restart_vectors(this, nc, nsites, sites) ccohort%bmort = rio_bmort_co(io_idx_co) ccohort%hmort = rio_hmort_co(io_idx_co) ccohort%cmort = rio_cmort_co(io_idx_co) - ccohort%fmort = rio_fmort_co(io_idx_co) ccohort%frmort = rio_frmort_co(io_idx_co) !Logging @@ -2268,6 +2463,39 @@ subroutine get_restart_vectors(this, nc, nsites, sites) end if + + ! Fill the site level diagnostics arrays + do i_scls = 1,nlevsclass + + do i_pft = 1, numpft + + sites(s)%fmort_rate_canopy(i_scls, i_pft) = rio_fmortrate_cano_siscpf(io_idx_si_scpf) + sites(s)%fmort_rate_ustory(i_scls, i_pft) = rio_fmortrate_usto_siscpf(io_idx_si_scpf) + sites(s)%imort_rate(i_scls, i_pft) = rio_imortrate_siscpf(io_idx_si_scpf) + sites(s)%fmort_rate_crown(i_scls, i_pft) = rio_fmortrate_crown_siscpf(io_idx_si_scpf) + sites(s)%fmort_rate_cambial(i_scls, i_pft) = rio_fmortrate_cambi_siscpf(io_idx_si_scpf) + sites(s)%term_nindivs_canopy(i_scls,i_pft) = rio_termnindiv_cano_siscpf(io_idx_si_scpf) + sites(s)%term_nindivs_ustory(i_scls,i_pft) = rio_termnindiv_usto_siscpf(io_idx_si_scpf) + sites(s)%growthflux_fusion(i_scls, i_pft) = rio_growflx_fusion_siscpf(io_idx_si_scpf) + + io_idx_si_scpf = io_idx_si_scpf + 1 + end do + + sites(s)%demotion_rate(i_scls) = rio_demorate_sisc(io_idx_si_sc) + sites(s)%promotion_rate(i_scls) = rio_promrate_sisc(io_idx_si_sc) + + io_idx_si_sc = io_idx_si_sc + 1 + end do + + + sites(s)%term_carbonflux_canopy = rio_termcflux_cano_si(io_idx_si) + sites(s)%term_carbonflux_ustory = rio_termcflux_usto_si(io_idx_si) + sites(s)%demotion_carbonflux = rio_democflux_si(io_idx_si) + sites(s)%promotion_carbonflux = rio_promcflux_si(io_idx_si) + sites(s)%imort_carbonflux = rio_imortcflux_si(io_idx_si) + sites(s)%fmort_carbonflux_canopy = rio_fmortcflux_cano_si(io_idx_si) + sites(s)%fmort_carbonflux_ustory = rio_fmortcflux_usto_si(io_idx_si) + sites(s)%old_stock = rio_old_stock_si(io_idx_si) sites(s)%status = rio_cd_status_si(io_idx_si) sites(s)%dstatus = rio_dd_status_si(io_idx_si) diff --git a/parameter_files/fates_params_14pfts.cdl b/parameter_files/fates_params_14pfts.cdl index 80c469e6a0..be2fcdda0b 100644 --- a/parameter_files/fates_params_14pfts.cdl +++ b/parameter_files/fates_params_14pfts.cdl @@ -527,9 +527,6 @@ variables: float fates_z0mr(fates_pft) ; fates_z0mr:units = "unitless" ; fates_z0mr:long_name = "Ratio of momentum roughness length to canopy top height" ; - float fates_alpha_FMC(fates_litterclass) ; - fates_alpha_FMC:units = "NA" ; - fates_alpha_FMC:long_name = "spitfire parameter related to fuel moisture content, Equation 6 Thonicke et al 2010" ; float fates_FBD(fates_litterclass) ; fates_FBD:units = "NA" ; fates_FBD:long_name = "spitfire parameter related to fuel bulk density, see SFMain.F90" ; @@ -560,6 +557,9 @@ variables: float fates_CWD_frac(fates_NCWD) ; fates_CWD_frac:units = "fraction" ; fates_CWD_frac:long_name = "fraction of woody (bdead+bsw) biomass destined for CWD pool" ; + float fates_drying_ratio ; + fates_drying_ratio:units = "NA" ; + fates_drying_ratio:long_name = "spitfire parameter, fire drying ratio for fuel moisture, alpha_FMC EQ 6 Thonicke et al 2010" ; float fates_durat_slope ; fates_durat_slope:units = "NA" ; fates_durat_slope:long_name = "spitfire parameter, fire max duration slope, Equation 14 Thonicke et al 2010" ; @@ -572,9 +572,6 @@ variables: float fates_fdi_b ; fates_fdi_b:units = "NA" ; fates_fdi_b:long_name = "spitfire parameter (unknown) " ; - float fates_fire_wind_max ; - fates_fire_wind_max:units = "m/min" ; - fates_fire_wind_max:long_name = "maximum wind speed expected by the fire model" ; float fates_fuel_energy ; fates_fuel_energy:units = "kJ/kg" ; fates_fuel_energy:long_name = "pitfire parameter, heat content of fuel" ; @@ -1165,8 +1162,6 @@ data: fates_z0mr = 0.055, 0.055, 0.055, 0.055, 0.055, 0.055, 0.055, 0.055, 0.055, 0.055, 0.055, 0.055, 0.055, 0.055 ; - fates_alpha_FMC = 0.0050769, 0.001, 0.0002754, 7.54e-05, 1.54e-05, 999 ; - fates_FBD = 4, 15.4, 16.8, 19.6, 999, 4 ; fates_low_moisture_Coeff = 1.15, 1.12, 1.09, 0.98, 0.8, 1.15 ; @@ -1187,6 +1182,8 @@ data: fates_CWD_frac = 0.045, 0.075, 0.21, 0.67 ; + fates_drying_ratio = 13000 ; + fates_durat_slope = -11.06 ; fates_fdi_a = 17.62 ; @@ -1195,8 +1192,6 @@ data: fates_fdi_b = 243.12 ; - fates_fire_wind_max = 45.718 ; - fates_fuel_energy = 18000 ; fates_max_durat = 240 ; diff --git a/parameter_files/fates_params_coastal_veg.cdl b/parameter_files/fates_params_coastal_veg.cdl index dddcc9d159..60c887fb83 100644 --- a/parameter_files/fates_params_coastal_veg.cdl +++ b/parameter_files/fates_params_coastal_veg.cdl @@ -23,9 +23,6 @@ variables: float fates_SAV(fates_litterclass) ; fates_SAV:units = "NA" ; fates_SAV:long_name = "spitfire parameter related to surface area to volume ratio, see SFMain.F90" ; - float fates_alpha_FMC(fates_litterclass) ; - fates_alpha_FMC:units = "NA" ; - fates_alpha_FMC:long_name = "spitfire parameter related to fuel moisture content, Equation 6 Thonicke et al 2010" ; float fates_history_height_bin_edges(fates_history_height_bins) ; fates_history_height_bin_edges:units = "m" ; fates_history_height_bin_edges:long_name = "Lower edges for height bins used in height-resolved history output" ; @@ -515,6 +512,9 @@ variables: float fates_phen_ncolddayslim(fates_scalar) ; fates_phen_ncolddayslim:units = "days" ; fates_phen_ncolddayslim:long_name = "day threshold exceedance for temperature leaf-drop" ; + float fates_drying_ratio ; + fates_drying_ratio:units = "NA" ; + fates_drying_ratio:long_name = "spitfire parameter, fire drying ratio for fuel moisture, alpha_FMC EQ 6 Thonicke et al 2010" ; float fates_durat_slope ; fates_durat_slope:units = "NA" ; fates_durat_slope:long_name = "spitfire parameter, fire max duration slope, Equation 14 Thonicke et al 2010" ; @@ -527,9 +527,6 @@ variables: float fates_fdi_b ; fates_fdi_b:units = "NA" ; fates_fdi_b:long_name = "spitfire parameter (unknown) " ; - float fates_fire_wind_max ; - fates_fire_wind_max:units = "m/min" ; - fates_fire_wind_max:long_name = "maximum wind speed expected by the fire model" ; float fates_fuel_energy ; fates_fuel_energy:units = "kJ/kg" ; fates_fuel_energy:long_name = "pitfire parameter, heat content of fuel" ; @@ -571,8 +568,6 @@ data: fates_SAV = 66, 13, 3.58, 0.98, 0.2, 66 ; - fates_alpha_FMC = 0.0050769, 0.001, 0.0002754, 7.54e-05, 1.54e-05, 999 ; - fates_history_height_bin_edges = 0, 0.1, 0.3, 0.5, 1, 2 ; fates_low_moisture_Coeff = 1.15, 1.12, 1.09, 0.98, 0.8, 1.15 ; @@ -939,6 +934,8 @@ data: fates_phen_ncolddayslim = 5 ; + fates_drying_ratio = 13000 ; + fates_durat_slope = -11.06 ; fates_fdi_a = 17.62 ; @@ -947,8 +944,6 @@ data: fates_fdi_b = 243.12 ; - fates_fire_wind_max = 45.718 ; - fates_fuel_energy = 18000 ; fates_max_durat = 240 ; diff --git a/parameter_files/fates_params_default.cdl b/parameter_files/fates_params_default.cdl index df152ce87b..0579b9405a 100644 --- a/parameter_files/fates_params_default.cdl +++ b/parameter_files/fates_params_default.cdl @@ -24,9 +24,6 @@ variables: float fates_SAV(fates_litterclass) ; fates_SAV:units = "NA" ; fates_SAV:long_name = "spitfire parameter related to surface area to volume ratio, see SFMain.F90" ; - float fates_alpha_FMC(fates_litterclass) ; - fates_alpha_FMC:units = "NA" ; - fates_alpha_FMC:long_name = "spitfire parameter related to fuel moisture content, Equation 6 Thonicke et al 2010" ; float fates_history_height_bin_edges(fates_history_height_bins) ; fates_history_height_bin_edges:units = "m" ; fates_history_height_bin_edges:long_name = "Lower edges for height bins used in height-resolved history output" ; @@ -560,6 +557,9 @@ variables: float fates_phen_ncolddayslim(fates_scalar) ; fates_phen_ncolddayslim:units = "days" ; fates_phen_ncolddayslim:long_name = "day threshold exceedance for temperature leaf-drop" ; + float fates_drying_ratio ; + fates_drying_ratio:units = "NA" ; + fates_drying_ratio:long_name = "spitfire parameter, fire drying ratio for fuel moisture, alpha_FMC EQ 6 Thonicke et al 2010" ; float fates_durat_slope ; fates_durat_slope:units = "NA" ; fates_durat_slope:long_name = "spitfire parameter, fire max duration slope, Equation 14 Thonicke et al 2010" ; @@ -572,9 +572,6 @@ variables: float fates_fdi_b ; fates_fdi_b:units = "NA" ; fates_fdi_b:long_name = "spitfire parameter (unknown) " ; - float fates_fire_wind_max ; - fates_fire_wind_max:units = "m/min" ; - fates_fire_wind_max:long_name = "maximum wind speed expected by the fire model" ; float fates_fuel_energy ; fates_fuel_energy:units = "kJ/kg" ; fates_fuel_energy:long_name = "pitfire parameter, heat content of fuel" ; @@ -617,8 +614,6 @@ data: fates_SAV = 66, 13, 3.58, 0.98, 0.2, 66 ; - fates_alpha_FMC = 0.0050769, 0.001, 0.0002754, 7.54e-05, 1.54e-05, 999 ; - fates_history_height_bin_edges = 0, 0.1, 0.3, 1, 3, 10 ; fates_low_moisture_Coeff = 1.15, 1.12, 1.09, 0.98, 0.8, 1.15 ; @@ -1063,6 +1058,8 @@ data: fates_phen_ncolddayslim = 5 ; + fates_drying_ratio = 13000 ; + fates_durat_slope = -11.06 ; fates_fdi_a = 17.62 ; @@ -1071,8 +1068,6 @@ data: fates_fdi_b = 243.12 ; - fates_fire_wind_max = 45.718 ; - fates_fuel_energy = 18000 ; fates_max_durat = 240 ;